Paper deep dive
Transformers converge to invariant algorithmic cores
Joshua S. Schiffman
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 7/20/2026, 10:15:02 AM
Summary
This paper introduces Algorithmic Core Extraction (ACE), a method to identify compact, necessary, and sufficient subspaces (algorithmic cores) within transformer models that implement task-specific computations. The study demonstrates that independently trained models converge to functionally equivalent but geometrically distinct cores. Key findings include the recovery of Markov chain dynamics in single-layer transformers, the emergence of rotational dynamics in modular addition tasks at 'grokking', and the identification of a shared, steerable one-dimensional axis governing subject-verb agreement across six large language models (GPT-2, LLaMA-3.1, Gemma-2, Qwen2.5). The work suggests that mechanistic interpretability should focus on these invariant computational structures rather than specific parameterizations.
Entities (12)
Relation Signals (13)
Algorithmic Core Extraction → identifies → Algorithmic Core
confidence 98% · ACE is introduced to isolate these subspaces, causally validate them, and recover the algorithms they implement
Algorithmic Core → isnecessaryandsufficientfor → Task Performance
confidence 95% · A core is sufficient if the projection preserves task performance, and necessary if its complement reduces it to near chance.
Transformer → convergesto → Algorithmic Core
confidence 92% · Transformers converge to invariant algorithmic cores
Subject-Verb Agreement → isgovernedby → Shared 1D Axis
confidence 92% · subject–verb agreement is governed by a single, steerable axis that aligns across architectures.
Qwen2.5 → exhibits → Shared 1D Axis
confidence 90% · In six language models... (Qwen2.5)
GPT-2 → exhibits → Shared 1D Axis
confidence 90% · In six language models... (GPT-2 Small/Medium/Large...)
Llama-3.1 → exhibits → Shared 1D Axis
confidence 90% · In six language models... (LLaMA-3.1...)
Gemma-2 → exhibits → Shared 1D Axis
confidence 90% · In six language models... (Gemma-2...)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Training selects for behavior, not circuitry: many weight configurations can implement the same function. Studying any single trained neural network thus risks describing accidents of one training run rather than the computation itself. This work shifts focus from what transformers happen to do to what they must do by extracting algorithmic cores, compact subspaces that are necessary and sufficient for a task and that recur across independently trained models. Here, Algorithmic Core Extraction (ACE) is introduced to isolate these subspaces, causally validate them, and recover the algorithms they implement across settings ranging from synthetic tasks to large-scale pretrained models. Markov-chain transformers embed three-dimensional cores in nearly orthogonal subspaces yet recover identical transition spectra. Modular-addition transformers form compact cyclic cores at grokking that later inflate under continued regularization, redundantly distributing the same computation across many functionally equivalent modes. This functional redundancy is found to accelerate the transition from memorization to generalization, yielding an inverse scaling law for grokking time. In six language models spanning more than two orders of magnitude in scale (GPT-2 Small/Medium/Large, LLaMA-3.1, Gemma-2, and Qwen2.5), subject-verb agreement is governed by a single, steerable axis that aligns across architectures. Flipping this axis inverts grammatical number throughout open-ended generation. Together these results suggest that beneath the apparent complexity of trained transformers lies a simpler, shared computational structure, and that targeting invariants rather than parameterizations may offer a more tractable path to mechanistic understanding and control. Code: this https URL
Tags
Links
- Source: https://arxiv.org/abs/2602.22600v2
- Canonical: https://arxiv.org/abs/2602.22600v2
Trouble viewing inline? Open PDF directly →
Full Text
75,524 characters extracted from source content.
Expand or collapse full text
Transformers Converge to Invariant Algorithmic Cores Joshua S. Schiffman Abstract Training selects for behavior, not circuitry: many weight configurations can implement the same function. Studying any single trained neural network thus risks describing accidents of one training run rather than the computation itself. This work shifts focus from what transformers happen to do to what they must do by extracting algorithmic cores, compact subspaces that are necessary and sufficient for a task and that recur across independently trained models. Here, Algorithmic Core Extraction (ACE) is introduced to isolate these subspaces, causally validate them, and recover the algorithms they implement across settings ranging from synthetic tasks to large-scale pretrained models. Markov-chain transformers embed three-dimensional cores in nearly orthogonal subspaces yet recover identical transition spectra. Modular-addition transformers form compact cyclic cores at grokking that later inflate under continued regularization, redundantly distributing the same computation across many functionally equivalent modes. This functional redundancy is found to accelerate the transition from memorization to generalization, yielding an inverse scaling law for grokking time. In six language models spanning more than two orders of magnitude in scale (GPT-2 Small/Medium/Large, LLaMA-3.1, Gemma-2, and Qwen2.5), subject–verb agreement is governed by a single, steerable axis that aligns across architectures. Flipping this axis inverts grammatical number throughout open-ended generation. Together these results suggest that beneath the apparent complexity of trained transformers lies a simpler, shared computational structure, and that targeting invariants rather than parameterizations may offer a more tractable path to mechanistic understanding and control. Code: https://github.com/joshseth/cores Mechanistic Interpretability, Transformers, Algorithmic Cores 1 Introduction A key obstacle to mechanistic interpretability (Elhage et al., 2021; Sharkey et al., 2025) is underdetermination: while training constrains model behavior – how inputs are mapped to outputs – it generally does not constrain how behavior is realized internally. This poses a fundamental challenge for interpretability: if mechanisms don’t generalize across realizations, which explanations are real? Such functional equivalence is routinely observed among independently trained artificial neural networks, and has been investigated in loss landscape geometry and model merging (Draxler et al., 2018; Garipov et al., 2018; Ainsworth et al., 2022), the nonidentifiability of mechanistic circuits (Méloux et al., 2025), representational similarity (Kornblith et al., 2019), and in the Rashomon effect (Breiman, 2001). Yet, this phenomenon is not restricted to neural networks and has been explored across scientific disciplines. In biology, it appears as degeneracy (Edelman and Gally, 2001) (e.g., in the genetic code), and in evolution as system drift (Schiffman and Ralph, 2022), where the wiring of a gene network changes but its function does not. In control theory, different realizations (Kalman, 1962, 1963) induce identical observable dynamics, and in physics, gauge symmetry indicates that many mathematical descriptions represent the same state. A natural response is to shift focus from individual realizations to equivalence classes, studying the invariants shared across them. If mechanistic explanations of language models are to generalize across random seeds (Gurnee et al., 2024), checkpoints, and architectures, they should be tethered to stable, implementation-invariant quantities rather than idiosyncratic details that vary across training runs. To explore this perspective, Algorithmic Core Extraction (ACE) is introduced to isolate algorithmic cores: low-dimensional subspaces that are necessary and sufficient for a task and shared across independent realizations. Applying ACE across three settings of escalating complexity demonstrates that functionally equivalent models can converge on compact, invariant mechanisms. In single-layer transformers (Vaswani et al., 2017), ACE recovers ground-truth Markov chain dynamics. In modular addition, it isolates the emergence of rotational dynamics at grokking. Finally, in six pretrained language models (spanning GPT-2, LLaMA-3.1, Gemma-2, and Qwen2.5) (Radford et al., 2019; Grattafiori et al., 2024; Riviere et al., 2024; Yang et al., 2025), ACE identifies a shared, one-dimensional core that causally steers subject–verb agreement during open-ended generation. This work contributes: (1) a conceptual framework for mechanistic interpretability that shifts focus from realization-specific circuitry to invariants; (2) ACE, a method for isolating compact subspaces that are causally necessary and sufficient for task performance; (3) evidence that these cores isolate interpretable mechanisms; (4) a theory linking functional equivalence, regularization, and grokking time; and (5) a steerable one-dimensional subject–verb agreement core shared by six distinct LLMs. 2 Methods The structure–function relationship is often many-to-one, but how many different structures can implement the same function? In linear system theory this can be answered with the Kalman decomposition (Kalman, 1962, 1963; Anderson et al., 1966; Kalman et al., 1969), which guarantees the existence of a minimal realization – a dynamical system that can be empirically recovered via balanced truncation (Moore, 1981). Algorithmic Core Extraction (ACE) operationalizes this principle for transformers by first extracting activation subspaces that are both highly active and relevant, then causally validating them with ablations, and finally fitting operators to identify the computations they perform (Appendix A). Extract. Fix a transformer layer with hidden dimension D. For N inputs, let ∈ℝN×DH ^N× D denote the mean-centered activations, with rows i⊤h_i , and let f:ℝD→ℝKf:R^D ^K map activations to K task-relevant outputs. Stack the Jacobians as ≔[(∂f/∂1)⊤⋯(∂f/∂N)⊤]⊤∈ℝNK×DJ [(∂ f/ _1) ·s(∂ f/ _N) ] ^NK× D. ACE finds directions that are jointly active and relevant by computing the SVD of their interaction:111When NK≫DNK\! \!D, use SVD of ⊤∈ℝD×DL ^D\!×\!D instead, where ⊤=⊤+εLL \!\!=\!H H\!+\! and ⊤=⊤ \!\!=\!J J. ⊤=⊤.HJ \;=\;U . The singular values quantify the joint activity and relevance of each direction and provide a principled criterion for rank selection. The algorithmic core is obtained by mapping the leading r columns of U back into activation space: ≔span(⊤r),C\, \,span\! (H U_r ), and QR decomposition yields an orthonormal basis ∈ℝD×rQ ^D× r and core projector ≔⊤P . Validate. A core is sufficient if the projection Ph preserves task performance, and necessary if its complement −h-Ph reduces it to near chance. Identify. A core’s computational structure is recovered by examining its coordinates =⊤z=Q h directly, or by fitting an operator A (e.g., t+1≈tz_t+1 _t by least squares) and inspecting its spectrum. 3 Algorithmic Core Necessity and Sufficiency The central goal of this manuscript is to determine whether low-dimensional subspaces, or algorithmic cores, within higher-dimensional trained transformers exist that are functionally necessary and sufficient for task performance. If so, are such cores shared across independently trained models, and do they admit simple mechanistic characterizations? Recovering algorithmic cores. The analysis begins in a fully controlled setting: three single-layer transformers (dmodel=64d_ model=64, dff=256d_ f=256, |V|=4|V|=4) trained with independent random seeds on a four-state Markov chain (Appendix B). Although each reached near Bayes-optimal test accuracy, their learned weights exhibited near-zero cosine similarity, indicating highly divergent parameterizations (Fig. 1A). To search for a shared internal representation, ACE was applied to each model’s 64-dimensional hidden state, successfully isolating a 3-dimensional algorithmic core. Ablations using all test data confirmed these cores were both necessary (removing the core drops accuracy to chance) and sufficient (retaining only the core preserves baseline accuracy) for the task (Fig. 1B; Table A1). Figure 1: Transformers trained on the same Markov task converge to a shared 3D causal core. Three one-layer transformers were trained with different random seeds on next-token prediction for a four-state Markov chain (Appendix B). (A) Learned weights differ substantially across runs measured by cosine similarity. (B) 3D core extracted from each 64D hidden state are necessary and sufficient under ablation (Table A1), and compared to optimal ∑iπimaxjTij _i _i _jT_ij and chance max() ( π) theoretical controls, with transition matrix T and stationary distribution π; points show individual accuracies and bars denote mean± . (C) Cores geometrically diverge with projector overlaps near zero and principal angles nearly orthogonal (Table A2). (D) Cores statistically align with mean cross-core CCAs reaching near unity (also see Table A2). (E) Dynamics fit in core recover the Markov chain nontrivial spectrum (Table A3). Geometric dissimilarity, statistical equivalence. To assess universality, each core recovered from the independently trained transformers was compared geometrically and statistically. Despite meeting equivalent causal criteria, cores were embedded in nearly orthogonal subspaces: projector overlap was 0.020.02–0.040.04, and principal angles ranged from 75∘75 –90∘90 (Fig. 1C; Table A2). Yet canonical correlation analysis (CCA) (Morcos et al., 2018) revealed nearly exact statistical alignment, with mean CCA correlations near 0.990.99 (Fig. 1D; Table A2). This suggests the cores encode the same information in different geometric coordinates – a signature of functionally equivalent yet structurally divergent realizations. Algorithmic cores encode Markov dynamics. To interpret what algorithm the cores implement, a linear operator was fit to next-token dynamics inside each core, and relative to “oracle” prediction, these operators achieved strong fits: Rcore2/Roracle2>0.98R^2_ core/R^2_ oracle>0.98 (Appendix B). Eigenvalues (the spectrum) of a linear operator determine its dynamics – such as oscillations and growth rates – so matching eigenvalues can indicate matching dynamics. Remarkably, the eigenvalues of each fit operator matched the non-trivial eigenvalues of the true Markov transition matrix to within a few percent (Fig. 1E; Table A3). This suggests that the recovered cores learned to efficiently encode Markov dynamics: trained transformers route inputs through a minimal, shared 3D subspace – that is necessary and sufficient for performance – and internally represents transition dynamics up to a change of coordinates. 4 Algorithmic Core Emergence and Evolution Because ACE is automated, it can recover learned computations without presupposing their form and can trace how they evolve during training. Modular addition is a natural test case: transformers trained on this task exhibit grokking (Power et al., 2022; Liu et al., 2022a), with high training accuracy preceding a delayed spike in test accuracy. Prior work showed that these models learn a Fourier “clock” algorithm (Nanda et al., 2023), but doing so required hypothesizing the mechanism a priori, designing targeted probes, and manually verifying circuits. Figure 2: Modular addition cores form at grokking and implement rotational mechanics. Three two-layer transformers were trained with different random seeds. (A) Test accuracy exhibits grokking (red, mean± , left y-axis) coincident with algorithmic core formation (gray, mean± , right y-axis). (B) After grokking, the recovered cores are necessary and sufficient under ablation. (C) Automated operator fits in core coordinates reveal the emergence of a cyclic mechanism: before grokking, eigenvalues scatter inside the unit circle, while at grokking they snap onto it, indicating discovery of a rotational mechanism. Cores crystallize at grokking. Three two-layer transformers (dmodel=128d_ model=128, dff=512d_ f=512, |V|=53|V|=53) were trained on modular addition (a+b≡cmod53a+b≡ c 53) for 2×1032× 10^3 epochs under weight decay regularization to encourage generalization. All models grokked: test accuracy remained near chance until spiking around epoch 800. Coincident with this delayed generalization, algorithmic cores crystallized – condensing into low-dimensional, ablation-defined necessary and sufficient subspaces (Fig. 2A,B; Appendix C). Blind recovery of rotational dynamics in cores. At each checkpoint, a linear operator was fit to the second-layer “shift” (add 1) dynamics in each extracted core. This revealed the emergence of a cyclic computational structure: at grokking, the operators’ eigenvalues snap onto the unit circle (Fig. 2C), indicating rotational dynamics capable of modular addition. Notably, this structure emerges directly from least-squares optimization in the core, without needing to prespecify an algorithmic form. However, while all three models converged to cyclic operators, the specific rotational modes (conjugate eigenvalue pairs) differed across runs – another instance of functional equivalence without structural identity (Chughtai et al., 2023; Zhong et al., 2023; Olah, 2025). Modular addition permits multiple valid modes and multiplicities, and models need not agree on which, nor how many, to use. Remarkably, even at grokking, each operator contained more rotational modes than the single mode minimally required – a hint of the redundancy that becomes extreme under extended training (Fig. 3). Figure 3: Extended training with weight decay inflates cores. Long-term training dynamics of transformers that grokked modular addition under different weight-decay schedules. (A) After grokking, core dimension continues to increase when weight decay is maintained (black, mean± ), but remains compact when weight decay is disabled post-grokking (purple). (B) Core inflation is driven by redundancy: the number of dimensions sufficient to preserve performance is stable, while the number whose removal reduces performance to chance increases. Lines depict means across models trained with weight decay fixed. (C) (Left) Dynamics fit in the terminal epoch reveal a saturated core operator when weight decay is maintained, in contrast to a more sparsely represented operator when weight decay is disabled. (Right) Rotational modes (conjugate eigenvalue pairs) around the unit circle increase with extended training under weight decay, whereas when weight decay is removed, mode counts remain stable. Cores inflate under extended training. Extending training to 2×1042× 10^4 epochs revealed an unexpected phenomenon: under continued weight decay, cores progressively inflated from approximately 1515 to 6060 dimensions. In contrast, disabling weight decay post-grokking kept cores more compact (Fig. 3A). This inflation is driven by a pronounced increase in redundant encoding. While the number of dimensions sufficient for task performance remained stable, the number of dimensions necessary to prevent chance-level performance expanded dramatically (Fig. 3B). Operator analysis reveals how this transformer "over-education" manifests: under continued weight decay, operators accumulated rotational modes. These approached the theoretical maximum of ⌊p/2⌋=26 p/2 =26 valid harmonic representations by the terminal epoch – far exceeding the minimally required single mode (Fig. 3C). Disabling weight decay prevented this proliferation: cores remained compact, mode counts stayed sparse, and operator structure remained stable. This suggests that weight decay may actively drive the transition from parsimonious algorithmic solutions to redundantly saturated representations. 4.1 Redundancy Drives Core Inflation and Grokking That a regularization penalty designed to simplify representations should instead inflate cores seems paradoxical. The behavior, however, emerges naturally from minimizing the weight norm within a highly redundant solution space. Furthermore, this interplay between redundancy and regularization also predicts the timing of grokking itself. Minimum norm requires maximum redundancy. After grokking, task loss is negligible and the gradient is dominated by weight decay (Varma et al., 2023), driving the network toward a minimum-norm solution. By Fourier symmetry (Chughtai et al., 2023), modular addition mod p admits ⌊p/2⌋ p/2 functionally equivalent modes, each a 2D rotation with phase θk≔2πk/p _k 2π k/p. Let ≥0 α≥ 0 denote mode amplitudes and ψ their label-contrasts, where ψk≔1−cosθk>0 _k 1- _k>0 is mode k’s contribution to the classification margin. If modes are encoded in approximately orthogonal parameter subspaces222If not orthogonal (or in superposition (Elhage et al., 2022)) with ≻0S\! \!0 mode-overlap, ‖22→⊤−1\| ψ\|_2^2→ ψ S^-1 ψ, reducing effective redundancy. then the weight norm satisfies ‖2≈‖2\|W\|^2≈\| α\|^2, while correct classification requires margin ⟨,⟩≥δ α, ψ ≥δ. Training thus implicitly solves min≥ 0‖2subject to⟨,⟩≥δ. _ α\,≥\,0\;\;\| α\|^2 to α, ψ ≥δ. By the Cauchy–Schwarz inequality, ‖2\| α\|^2 is minimized when ∥ α ψ – that is, when every mode is active – with the optimal solution ∗=(δ/‖22) α = (δ/\| ψ\|_2^2 ) ψ. Weight decay thus acts as a redistribution force: rather than simplifying the representation, it spreads weight across all valid solutions. Disabling weight decay removes this pressure, consistent with observations in Fig. 3. Functional equivalence accelerates grokking. The same redistribution pressure governs the speed of grokking. Define the grokking delay τgrok≔τgen−τmem _grok _gen- _mem as the time between memorization and generalization, and model the transition to generalization as the margin m(t)≔⟨(t),⟩m(t) α(t), ψ reaching threshold δ. After memorization, with task gradients largely vanished and weight decay (ω) dominating, the expected margin trajectory follows (Appendix E.1) m˙(t)=−ωm(t)+cω‖22. m(t)\;=\;-ω\,m(t)\;+\;c\,ω\,\| ψ\|_2^2. Crucially, functional equivalence makes the margin-driving direction additive across modes, giving ‖22∝p\| ψ\|_2^2 p.333Using ∑k=1⌊p/2⌋cosθk=∑k=1⌊p/2⌋cos2θk=−12 _k=1 p/2 _k= _k=1 p/2 2 _k=- 12 and expanding (1−cosθk)2(1- _k)^2 gives ‖22=34p\| ψ\|_2^2= 34p. Each redundant mode amplifies the mean-drift velocity toward generalization, consistent with the multiple active modes observed at grokking (Fig. 2C). When p<dmodelp<d_model the initial memorized solution has negligible margin (m(0)≈0m(0)≈ 0), and grokking occurs when m(τ)=δm(τ)=δ. Solving for the expected grokking time delay yields an expression that linearizes for high redundancy (p≫pcritp p_crit) into a simple inverse scaling law: τgrok(p)=−Ωlog(1−pcritp)≈Ωpcritp∝1ωp. _grok(p)\;=\;-\, \, \! (1- p_critp )\;≈\; \,p_critp\; \; 1ω\,p. Two empirical constants govern this expression: an optimizer constant Ω∝ω−1 ω^-1 that sets the timescale of grokking when it occurs, and an architectural constant pcritp_crit that determines whether it can occur at all. Grokking time thus shrinks with both weight decay and functional redundancy. These predictions are validated by sweeping ω and p in transformers (Fig. 4; Appendix E.2). Figure 4: Grokking time scales inversely with redundancy. (Left) Time to grok after memorization τ scales inversely with weight decay ω; consistent with prior observations (Liu et al., 2022b). The observed fit (τ∝ω−1.02τ ω^-1.02, red) matches the theoretical prediction ω−1ω^-1 (gray). (Right) Grokking time scales inversely with redundancy p. The ODE solution (green; R2>0.99R^2>0.99) captures both the inverse scaling at large p and the divergence near pcritp_crit. Fit parameters: Ω^≈2,770 ≈ 2,770, p^crit≈23 p_crit≈ 23. Points represent mean± for 12 random seeds (Appendix E.2). Summary. The algorithmic core framework – automated operator extraction from causally defined, low-dimensional core subspaces – can mechanistically characterize and trace the evolution of computations transformers learn throughout training. In modular addition, the extracted cores exhibit rotational dynamics consistent with the task’s cyclic structure, crystallize at grokking, and inflate under extended weight decay. This inflation reflects transformers converging on the optimal weighting strategy under regularization: to distribute weight across all functionally equivalent representations. This same pressure – regularization utilizing redundancy – predicts the speed of grokking, explaining the transition from memorization to generalization. The next question is whether these tools scale to larger and more complex systems. Figure 5: Subject–verb agreement is mediated by a shared 1D core across LLMs. The core framework was applied to GPT-2 Small, Medium, Large, Gemma-2, LLaMA-3.1, and Qwen2.5 to isolate a low-dimensional mechanism for number agreement. (A) Layer sweep: agreement performance (AUC) vs. normalized layer depth, averaged across LLMs (lines) with per-model measurements (markers) and shaded min–max bands. Agreement performance is the probability that the model assigns a higher plural-vs-singular verb-preference score to a plural than to a singular prompt. (B) Projecting last-token hidden states onto the core produces a nearly linear control axis for the singular–plural logit margin; per-model affine fits are shown after z scoring. (C) Removing the core degrades agreement, while flipping it inverts verb preference. Box plots summarize prompt-level agreement scores under perturbation; reported p-values combine per-model paired Wilcoxon tests using Fisher’s method. Base: We hold these truths to be self-evident: that all men are created equal, that they are endowed by their Creator with certain unalienable Rights [i.e., without a priori moral rights], and that among these are Life, Liberty and the pursuit of Happiness.” Core Steering: We hold these truths to be self-evident: that all men is created equal, that they is endowed by their Creator with certain unalienable Rights [i.e., the right to life], that among these is Life[.]” Base: As a new field of research, artificial intelligence has already delivered us with its revolutionary ways to solve complex medical issues. Its potential to address more such problems can boost the healthcare system… Core Steering: As a new field of research, artificial intelligence have already made great strides to improve our lives. AI has been instrumental in providing a more efficient… But how will we know what we has the potential to do… Figure 6: Core steering induces systematic agreement violations in open-ended generation. Prompted text from Base or Core Steering, with select violations highlighted. 5 Scaling ACE to LLMs: A Universal 1D Core The preceding experiments establish the ACE framework in highly controlled, synthetic settings. The critical question is whether the ACE framework scales beyond toy models to govern complex behaviors in production-scale models. To establish an empirical foothold on this question, ACE was applied to six pretrained language models spanning four distinct families: GPT-2 Small, Medium, and Large (117M, 345M, and 774M parameters) (Radford et al., 2019; Wolf et al., 2020); LLaMA-3.1 (8B) (Grattafiori et al., 2024); Gemma-2 (9B) (Riviere et al., 2024); and Qwen2.5 (32B) (Yang et al., 2025). These models differ in architecture, training corpus, and tokenization, and span more than two orders of magnitude in parameter count. The target task is subject–verb number agreement, a tractable linguistic computation with clear ground-truth labels (singular vs. plural subject) and a well-defined behavioral output (verb selection); admitting systematic evaluation via controlled prompts and a scalar verb-preference score (Linzen et al., 2016; Marvin and Linzen, 2018; Finlayson et al., 2021) (Appendix D). Localizing a shared 1D agreement core. To localize the agreement mechanism, candidate cores were extracted at each layer and evaluated via causal ablations. Across all six models, early layers exhibited minimal causal influence, but a highly potent core consistently emerged in the late layers (Fig. 5A). At the layer of maximal effect, this agreement core is remarkably one-dimensional – a single axis separated from all remaining directions by a large spectral gap (Table A5). Causal validation and control. Observationally, this axis behaves as a graded number coordinate: projection onto it predicts the singular–plural logit margin across models (Fig. 5B), aligning with the linear representation hypothesis (Park et al., 2023). However, because subspace projections alone can be deceptive (Belinkov, 2022; Makelov et al., 2023), claims here are strictly grounded in causal ablations. Despite its compact size, this single axis is sufficient (retaining it preserves agreement; AUC ≥0.91≥ 0.91), necessary (removing it collapses agreement below chance; AUC ≤0.25≤ 0.25), and directionally controllable. Reflecting activations through this axis inverts verb preferences, inducing strong disagreement with the subject (AUC ≤0.10≤ 0.10; Table A6). At the prompt level, for instance, core inversion on “The key next to the cabinets” drives ℙ(is)P(is) from 0.510.51 down to 0.010.01, while boosting ℙ(are)P(are) from 0.060.06 to 0.710.71 (Fig. 5C). Alignment across LLMs. Projecting last-token hidden states onto each model’s agreement core (Fig. 5B) yields a signed grammatical-number coordinate that tracks verb preference. Because cores are one-dimensional, cross-model alignment reduces to fixing a sign convention and comparing projected coordinates. Within the GPT-2 family – three models that share an architecture and training procedure – coordinates align tightly (Spearman’s ρ=0.88ρ=0.88–0.920.92; Pearson’s r=0.92 r=0.92–0.970.97). More strikingly, alignment persists across families: between Qwen2.5, LLaMA-3.1, Gemma-2, and the GPT-2 models, Spearman correlations range from 0.590.59 to 0.930.93 and Pearson correlations from 0.620.62 to 0.940.94 (Table A4). The strongest cross-family correlations (Qwen × Gemma: ρ=0.93ρ=0.93; Gemma × GPT-2 Medium: ρ=0.89ρ=0.89) approach the within-family ceiling, indicating that the agreement core encodes grammatical number in a way that is largely independent of architecture, tokenization, training corpus, and scale. Core steering inverts grammar in open-ended text. A stronger test of the core’s role is whether it governs agreement throughout autoregressive generation, where each token conditions subsequent predictions. To test this, the core-axis intervention was applied adaptively at each decoding step. Modulating the intervention strength based on each token’s sensitivity to number agreement, leaving irrelevant tokens untouched (Appendix D), induced systematic agreement violations across all six models (Fig. 6). Singular subjects recruited plural verbs, plural contexts shifted toward the singular, and errors cascaded as toggling the number variable corrupted downstream predictions. Crucially, the effect generalized well beyond the specific verbs (is/are/was/were) used to define the initial preference score. The emergence of agreement failures in entirely different word classes supports the interpretation that the core encodes a global grammatical-number variable, rather than a narrow, verb-specific heuristic. Summary. Subject–verb agreement in language models is governed by a 1D causal subspace localized to late layers. This core is necessary, sufficient, and controllable, and its coordinates align across six models from four families. 6 Discussion These results suggest that transformer computations may be governed by low-dimensional mechanisms that recur across independent training runs despite substantial variation in learned parameters. These findings have implications for how we conceptualize mechanistic interpretability. Invariance, not sparsity or circuitry. Mechanistic interpretability has largely studied implementations, such as circuits of attention heads and neurons (Elhage et al., 2021; Olah et al., 2020; Wang et al., 2022; Ameisen et al., 2025; Lindsey et al., 2025), or sparse decompositions of activations into interpretable features (Cunningham et al., 2023; Bricken et al., 2023; Dunefsky et al., 2024; Templeton et al., 2024). Such descriptions can be highly precise, but they face a conceptual challenge: they may be implementation-specific. Two models might compute the same function using entirely different circuits and coordinate systems (Méloux et al., 2025; Fel et al., 2025). The core framework shifts the explanatory target from implementation to invariant. The motivation for sparse features parallels a classical aim in linear algebra: diagonalization. But the fundamental power of diagonalization lies not in sparsity per se, but in revealing invariants – eigenvalues preserved under change of basis. Sparsity is basis-dependent; invariants are not. Likewise, circuits and sparse features describe coordinates of implementation, while cores identify the causal subspaces and dynamics preserved across implementations. Where features or circuits recur across models, perhaps identified via cross-coders (Lindsey et al., 2024), the approaches converge. Where they diverge, invariance provides a reliable criterion for distinguishing structural essence from artifact. Cores as internal world models. The observation that independent models converge to the same invariant structure raises a natural question: what anchors these shared representations? If these cores are not artifacts of the architecture or training run, they might reflect the data-generating process itself. When algorithmic cores recover ground-truth task structure – Markov transition spectra, cyclic operators for modular arithmetic – they encode not merely input–output mappings but internal representations of the generative process underlying the data (Li et al., 2022; Gurnee and Tegmark, 2023; Huh et al., 2024). This aligns with two classical ideas: the good regulator theorem (Conant and Ross Ashby, 1970) and the internal model principle (Francis and Wonham, 1976) from control theory, which hold that any system achieving optimal prediction must contain a model of its environment. When a core is isomorphic to the task-generating process, interpretability may be viewed as a form of internal-model recovery. Redundancy accelerates grokking. Once a model reaches perfect training accuracy, it enters a highly degenerate zero-loss manifold in parameter space (Bushnaq et al., 2024), populated by many functionally equivalent solutions. Weight decay then biases stochastic exploration along this manifold toward a minimum-norm, maximum-margin solution. Because the target task admits multiple functionally equivalent realizations, the corrective pressure from weight decay accumulates across valid modes rather than acting on a single narrow solution, accelerating the expected trajectory toward generalization. As the network traverses the continuous margin space, it eventually crosses the discrete classification threshold, producing a sharp jump in test accuracy even though the underlying trajectory in weight space remains smooth. Speculatively, scaling may produce capability jumps once models have enough capacity to realize large classes of functionally equivalent solutions. An analogous phenomenon appears in evolutionary genetics, where robustness creates extended webs of phenotype-preserving genotypes that facilitate the discovery of new functions via neutral drift (Wagner, 2008, 2012). A practical consequence is that cores are most compact immediately after grokking and subsequently inflate. This suggests a natural interpretability window: annealing weight decay toward zero shortly after task convergence may help preserve the most compact solution. System drift and model merging. System drift describes how a gene network can preserve its phenotype while its underlying genetic wiring diverges, effectively drifting through a neutral space. Because the set of functionally equivalent realizations is not generally convex or closed under recombination, mixing divergent solutions often produces hybrid incompatibility (Schiffman and Ralph, 2022). Transformers exhibit an analogous pattern: models trained from different initializations implement identical cores embedded in nearly orthogonal subspaces, revealing substantial representational drift despite functional equivalence. This orthogonality implies that naïve weight interpolation between geometrically divergent models moves off the solution manifold, consistent with empirical difficulties in model merging (Garipov et al., 2018; Ainsworth et al., 2022). By contrast, extracting and aligning algorithmic cores may offer a principled diagnostic for merge-compatibility and a potential coordinate system for successful recombination. Limitations and future directions. Whether cores remain low-dimensional for multi-step reasoning tasks is untested. The agreement core, however, remains one-dimensional across six models spanning four architectures and over two orders of magnitude in scale (117M to 32B parameters), suggesting core dimensionality may not depend on model scale. This is compatible with the empirical success of LoRA (Hu et al., 2022), which often achieves large behavioral changes via low-dimensional weight updates. Extracting task-specific cores from multifunctional models also requires framing precise mechanistic inquiries; this work demonstrates this for subject–verb agreement, but systematic approaches to task decomposition remain open. The extraction procedure itself admits natural extensions: nonlinear dimensionality reduction in place of the active component, learned probes in place of Jacobians, and Koopman operator approximations (Brunton et al., 2021) for tasks with nonlinear dynamics. More broadly, the relevant invariants for complex tasks are not obvious a priori; future work might discover them empirically by asking what core properties are shared across independently trained models. Finally, methods that identify causally effective subspaces may enable more targeted model control: this could support auditing and debugging, but also creates misuse risks if used to induce systematic errors or circumvent intended behaviors. Conclusion. These results point toward a view of transformer computation as organized around low-dimensional invariants: subspaces that are preserved across training runs, necessary and sufficient for task performance, and structured in ways that mirror the tasks themselves. If this view is approximately correct, interpretability efforts may benefit from targeting such invariants – seeking the computational essence that recurs across realizations rather than the implementation details that vary. The algorithmic core is one operationalization of this intuition. Whether it scales to the complexity of contemporary language models remains to be seen, but the guiding principle – focus on what is preserved, not what is particular – may prove durable. Code Availability Code to reproduce all analyses and figures is available at https://github.com/joshseth/cores. Acknowledgements I would like to thank Drs. Alison Pickover and Dan Landau for their support. References S. K. Ainsworth, J. Hayase, and S. Srinivasa (2022) Git re-basin: merging models modulo permutation symmetries. arXiv preprint arXiv:2209.04836. Cited by: §1, §6. E. Ameisen, J. Lindsey, A. Pearce, W. Gurnee, N. L. Turner, B. Chen, C. Citro, D. Abrahams, S. Carter, B. Hosmer, J. Marcus, M. Sklar, A. Templeton, T. Bricken, C. McDougall, H. Cunningham, T. Henighan, A. Jermyn, A. Jones, A. Persic, Z. Qi, T. Ben Thompson, S. Zimmerman, K. Rivoire, T. Conerly, C. Olah, and J. Batson (2025) Circuit tracing: revealing computational graphs in language models. Transformer Circuits Thread. External Links: Link Cited by: §6. B. Anderson, R. Newcomb, R. Kalman, and D. Youla (1966) Equivalence of linear time-invariant dynamical systems. Journal of the Franklin Institute 281 (5), p. 371–378. Cited by: Appendix A, §2. Y. Belinkov (2022) Probing classifiers: promises, shortcomings, and advances. Computational Linguistics 48 (1), p. 207–219. Cited by: §5. R. Bellman and K. J. Åström (1970) On structural identifiability. Mathematical Biosciences 7 (3-4), p. 329–339. Cited by: Appendix A. L. Breiman (2001) Statistical modeling: the two cultures. Statistical Science. Cited by: §1. T. Bricken, A. Templeton, J. Batson, B. Chen, A. Jermyn, T. Conerly, N. Turner, C. Anil, C. Denison, A. Askell, R. Lasenby, Y. Wu, S. Kravec, N. Schiefer, T. Maxwell, N. Joseph, Z. Hatfield-Dodds, A. Tamkin, K. Nguyen, B. McLean, J. E. Burke, T. Hume, S. Carter, T. Henighan, and C. Olah (2023) Towards monosemanticity: decomposing language models with dictionary learning. Transformer Circuits Thread. Note: https://transformer-circuits.pub/2023/monosemantic-features/index.html Cited by: §6. S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz (2021) Modern koopman theory for dynamical systems. arXiv preprint arXiv:2102.12086. Cited by: §6. L. Bushnaq, J. Mendel, S. Heimersheim, D. Braun, N. Goldowsky-Dill, K. Hänni, C. Wu, and M. Hobbhahn (2024) Using degeneracy in the loss landscape for mechanistic interpretability. arXiv preprint arXiv:2405.10927. Cited by: §6. B. Chughtai, L. Chan, and N. Nanda (2023) A toy model of universality: reverse engineering how networks learn group operations. In International Conference on Machine Learning, p. 6243–6267. Cited by: §4, §4.1. R. C. Conant and W. Ross Ashby (1970) Every good regulator of a system must be a model of that system. International journal of systems science 1 (2), p. 89–97. Cited by: §6. H. Cunningham, A. Ewart, L. Riggs, R. Huben, and L. Sharkey (2023) Sparse autoencoders find highly interpretable features in language models. arXiv preprint arXiv:2309.08600. Cited by: §6. F. Draxler, K. Veschgini, M. Salmhofer, and F. Hamprecht (2018) Essentially no barriers in neural network energy landscape. In International conference on machine learning, p. 1309–1318. Cited by: §1. J. Dunefsky, P. Chlenski, and N. Nanda (2024) Transcoders find interpretable llm feature circuits. Advances in Neural Information Processing Systems 37, p. 24375–24410. Cited by: §6. G. M. Edelman and J. A. Gally (2001) Degeneracy and complexity in biological systems. Proceedings of the national academy of sciences 98 (24), p. 13763–13768. Cited by: §1. N. Elhage, T. Hume, C. Olsson, N. Schiefer, T. Henighan, S. Kravec, Z. Hatfield-Dodds, R. Lasenby, D. Drain, C. Chen, R. Grosse, S. McCandlish, J. Kaplan, D. Amodei, M. Wattenberg, and C. Olah (2022) Toy models of superposition. Note: Transformer CircuitsOnline External Links: Link Cited by: footnote 2. N. Elhage, N. Nanda, C. Olsson, T. Henighan, N. Joseph, B. Mann, A. Askell, Y. Bai, A. Chen, T. Conerly, N. DasSarma, D. Drain, D. Ganguli, Z. Hatfield-Dodds, D. Hernandez, A. Jones, J. Kernion, L. Lovitt, K. Ndousse, D. Amodei, T. Brown, J. Clark, J. Kaplan, S. McCandlish, and C. Olah (2021) A mathematical framework for transformer circuits. Transformer Circuits Thread. Note: https://transformer-circuits.pub/2021/framework/index.html Cited by: §1, §6. T. Fel, E. S. Lubana, J. S. Prince, M. Kowal, V. Boutin, I. Papadimitriou, B. Wang, M. Wattenberg, D. Ba, and T. Konkle (2025) Archetypal sae: adaptive and stable dictionary learning for concept extraction in large vision models. arXiv preprint arXiv:2502.12892. Cited by: §6. M. Finlayson, A. Mueller, S. Gehrmann, S. M. Shieber, T. Linzen, and Y. Belinkov (2021) Causal analysis of syntactic agreement mechanisms in neural language models. In Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), p. 1828–1843. Cited by: §5. B. A. Francis and W. M. Wonham (1976) The internal model principle of control theory. Automatica 12 (5), p. 457–465. Cited by: §6. T. Garipov, P. Izmailov, D. Podoprikhin, D. P. Vetrov, and A. G. Wilson (2018) Loss surfaces, mode connectivity, and fast ensembling of DNNs. Advances in neural information processing systems 31. Cited by: §1, §6. A. Grattafiori, A. Dubey, A. Jauhri, A. Pandey, A. Kadian, A. Al-Dahle, A. Letman, A. Mathur, A. Schelten, A. Vaughan, et al. (2024) The llama 3 herd of models. External Links: 2407.21783, Link Cited by: §1, §5. W. Gurnee, T. Horsley, Z. C. Guo, T. R. Kheirkhah, Q. Sun, W. Hathaway, N. Nanda, and D. Bertsimas (2024) Universal neurons in gpt2 language models. arXiv preprint arXiv:2401.12181. Cited by: §1. W. Gurnee and M. Tegmark (2023) Language models represent space and time. arXiv preprint arXiv:2310.02207. Cited by: §6. E. J. Hu, yelong shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen (2022) LoRA: low-rank adaptation of large language models. In International Conference on Learning Representations, External Links: Link Cited by: §6. M. Huh, B. Cheung, T. Wang, and P. Isola (2024) The platonic representation hypothesis. arXiv preprint arXiv:2405.07987. Cited by: §6. R. E. Kalman, P. L. Falb, and M. A. Arbib (1969) Topics in mathematical system theory. McGraw-Hill, New York (English). External Links: ISBN 0754321069 Cited by: Appendix A, §2. R. E. Kalman (1962) Canonical structure of linear dynamical systems. Proceedings of the National Academy of Sciences 48 (4), p. 596–600. Cited by: Appendix A, Appendix A, §1, §2. R. E. Kalman (1963) Mathematical description of linear dynamical systems. Journal of the Society for Industrial and Applied Mathematics, Series A: Control 1 (2), p. 152–192. Cited by: Appendix A, §1, §2. S. Kornblith, M. Norouzi, H. Lee, and G. Hinton (2019) Similarity of neural network representations revisited. In International conference on machine learning, p. 3519–3529. Cited by: §1. K. Li, A. K. Hopkins, D. Bau, F. Viégas, H. Pfister, and M. Wattenberg (2022) Emergent world representations: exploring a sequence model trained on a synthetic task. arXiv preprint arXiv:2210.13382. Cited by: §6. J. Lindsey, W. Gurnee, E. Ameisen, B. Chen, A. Pearce, N. L. Turner, C. Citro, D. Abrahams, S. Carter, B. Hosmer, J. Marcus, M. Sklar, A. Templeton, T. Bricken, C. McDougall, H. Cunningham, T. Henighan, A. Jermyn, A. Jones, A. Persic, Z. Qi, T. B. Thompson, S. Zimmerman, K. Rivoire, T. Conerly, C. Olah, and J. Batson (2025) On the biology of a large language model. Transformer Circuits Thread. External Links: Link Cited by: §6. J. Lindsey, A. Templeton, J. Marcus, T. Conerly, J. Batson, and C. Olah (2024) Sparse crosscoders for cross-layer features and model diffing. Transformer Circuits Thread, p. 3982–3992. Cited by: §6. T. Linzen, E. Dupoux, and Y. Goldberg (2016) Assessing the ability of lstms to learn syntax-sensitive dependencies. Transactions of the Association for Computational Linguistics 4, p. 521–535. Cited by: §5. Z. Liu, O. Kitouni, N. S. Nolte, E. Michaud, M. Tegmark, and M. Williams (2022a) Towards understanding grokking: an effective theory of representation learning. Advances in Neural Information Processing Systems 35, p. 34651–34663. Cited by: §4. Z. Liu, E. J. Michaud, and M. Tegmark (2022b) Omnigrok: grokking beyond algorithmic data. arXiv preprint arXiv:2210.01117. Cited by: Figure 4, Figure 4. A. Makelov, G. Lange, and N. Nanda (2023) Is this the subspace you are looking for? an interpretability illusion for subspace activation patching. arXiv preprint arXiv:2311.17030. Cited by: §5. R. Marvin and T. Linzen (2018) Targeted syntactic evaluation of language models. In Proceedings of the 2018 conference on empirical methods in natural language processing, p. 1192–1202. Cited by: §5. M. Méloux, S. Maniu, F. Portet, and M. Peyrard (2025) Everything, everywhere, all at once: is mechanistic interpretability identifiable?. arXiv preprint arXiv:2502.20914. Cited by: §1, §6. B. Moore (1981) Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Transactions on Automatic Control 26 (1), p. 17–32. External Links: Document Cited by: Appendix A, §2. A. Morcos, M. Raghu, and S. Bengio (2018) Insights on representational similarity in neural networks with canonical correlation. Advances in neural information processing systems 31. Cited by: §3. N. Nanda, L. Chan, T. Lieberum, J. Smith, and J. Steinhardt (2023) Progress measures for grokking via mechanistic interpretability. arXiv preprint arXiv:2301.05217. Cited by: §4. C. Olah, N. Cammarata, L. Schubert, G. Goh, M. Petrov, and S. Carter (2020) Zoom in: an introduction to circuits. Distill. Note: https://distill.pub/2020/circuits/zoom-in External Links: Document Cited by: §6. C. Olah (2025) Note: Transformer Circuits External Links: Link Cited by: §4. K. Park, Y. J. Choe, and V. Veitch (2023) The linear representation hypothesis and the geometry of large language models. arXiv preprint arXiv:2311.03658. Cited by: §5. A. Power, Y. Burda, H. Edwards, I. Babuschkin, and V. Misra (2022) Grokking: generalization beyond overfitting on small algorithmic datasets. arXiv preprint arXiv:2201.02177. Cited by: §4. A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever, et al. (2019) Language models are unsupervised multitask learners. OpenAI blog 1 (8), p. 9. Cited by: §1, §5. M. Riviere, S. Pathak, P. G. Sessa, C. Hardin, S. Bhupatiraju, L. Hussenot, T. Mesnard, B. Shahriari, A. Ramé, et al. (2024) Gemma 2: improving open language models at a practical size. arXiv preprint arXiv:2408.00118. Cited by: §1, §5. J. S. Schiffman and P. L. Ralph (2022) System drift and speciation. Evolution 76 (2), p. 236–251. Cited by: §1, §6. L. Sharkey, B. Chughtai, J. Batson, J. Lindsey, J. Wu, L. Bushnaq, N. Goldowsky-Dill, S. Heimersheim, A. Ortega, J. Bloom, et al. (2025) Open problems in mechanistic interpretability. arXiv preprint arXiv:2501.16496. Cited by: §1. A. Templeton, T. Conerly, J. Marcus, J. Lindsey, T. Bricken, B. Chen, A. Pearce, C. Citro, E. Ameisen, A. Jones, H. Cunningham, N. L. Turner, C. McDougall, M. MacDiarmid, C. D. Freeman, T. R. Sumers, E. Rees, J. Batson, A. Jermyn, S. Carter, C. Olah, and T. Henighan (2024) Scaling monosemanticity: extracting interpretable features from claude 3 sonnet. Transformer Circuits Thread. External Links: Link Cited by: §6. V. Varma, R. Shah, Z. Kenton, J. Kramár, and R. Kumar (2023) Explaining grokking through circuit efficiency. arXiv preprint arXiv:2309.02390. Cited by: §4.1. A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. Advances in neural information processing systems 30. Cited by: §1. A. Wagner (2008) Robustness and evolvability: a paradox resolved. Proceedings of the Royal Society B: Biological Sciences 275 (1630), p. 91–100. Cited by: §6. A. Wagner (2012) The role of robustness in phenotypic adaptation and innovation. Proceedings of the Royal Society B: Biological Sciences 279 (1732), p. 1249–1258. Cited by: §6. K. Wang, A. Variengien, A. Conmy, B. Shlegeris, and J. Steinhardt (2022) Interpretability in the wild: a circuit for indirect object identification in gpt-2 small. arXiv preprint arXiv:2211.00593. Cited by: §6. T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz, et al. (2020) Transformers: state-of-the-art natural language processing. In Proceedings of the 2020 conference on empirical methods in natural language processing: system demonstrations, p. 38–45. Cited by: §5. A. Yang, B. Yang, B. Zhang, B. Hui, B. Zheng, B. Yu, C. Li, D. Liu, F. Huang, H. Wei, H. Lin, J. Yang, J. Tu, J. Zhang, J. Yang, J. Yang, J. Zhou, J. Lin, K. Dang, K. Lu, K. Bao, K. Yang, L. Yu, M. Li, M. Xue, P. Zhang, Q. Zhu, R. Men, R. Lin, T. Li, T. Tang, T. Xia, X. Ren, X. Ren, Y. Fan, Y. Su, Y. Zhang, Y. Wan, Y. Liu, Z. Cui, Z. Zhang, and Z. Qiu (2025) Qwen2.5 technical report. External Links: 2412.15115, Link Cited by: §1, §5. Z. Zhong, Z. Liu, M. Tegmark, and J. Andreas (2023) The clock and the pizza: two stories in mechanistic explanation of neural networks. Advances in neural information processing systems 36, p. 27223–27250. Cited by: §4. Appendix Appendix A Algorithmic Core Extraction Functional equivalence and minimal realizations. The structure–function relationship is often many-to-one (Bellman and Åström, 1970): there is more than one way to realize a behavior. But how many different structures can realize identical input–output functions? Can the space of functionally equivalent structures be characterized? In linear system theory, this question has an exact answer (Kalman, 1962). Consider a linear time-invariant system with hidden state ∈ℝnx ^n, input ∈ℝmu ^m, and output ∈ℝℓy : ˙ x =+, =Ax+Bu, =. =Cx. The system’s input–output behavior is fully determined by its impulse response (t)≔eAt ζ(t) AtB. Two systems with different weights (,,)(A,B,C) and (~,~,~)( A, B, C) are functionally equivalent if they produce identical outputs for all inputs – that is, if their impulse responses match ((t)=~(t) ζ(t)= ζ(t)). For systems of equal dimension, functional equivalence corresponds exactly to coordinate change: (,,)(A,B,C) and (−1,,−1)(VAV^-1,VB,CV^-1) share the same impulse response for any invertible V. But systems of different sizes can also be functionally equivalent if some internal states are either unreachable (unaffected by input) or unobservable (irrelevant to the output). The Kalman decomposition makes this precise, partitioning any system’s state space into four subspaces according to reachability and observability (Kalman, 1962, 1963; Anderson et al., 1966; Kalman et al., 1969). Only states that are both reachable and observable contribute to input–output behavior; the rest represent degrees of freedom that can vary without affecting function. This decomposition guarantees the existence of a minimal realization – the smallest-dimensional system that reproduces an input–output map, unique up to coordinate change – and enables extracting it. These results from system theory conceptually motivate the methods developed in this manuscript. Algorithmic core extraction. The goal here is an analogous decomposition for transformers. The Kalman decomposition provides an exact algebraic characterization for linear systems; for transformers, no such closed-form decomposition exists, but the principle can be applied empirically: identify directions that are both input-driven (active) and output-relevant (relevant). If the system were linear, this would reduce to balanced truncation (Moore, 1981), a technique in model reduction that finds coordinates in which reachability and observability are aligned. Here, ACE (Algorithmic Core Extraction) operationalizes this approach for artificial neural networks. Let ∈ℝN×DH ^N× D denote mean-centered hidden activations at a transformer layer of interest, with rows i⊤∈ℝ1×Dh_i ^1× D for each of N inputs, to define active directions. To quantify relevant directions, let f:ℝD→ℝKf ^D ^K map activations to task-relevant outputs and let ∈ℝNK×DJ ^NK× D stack the N Jacobians ∂f/∂i∈ℝK×D∂ f/ _i ^K× D as row blocks. To find directions that are jointly active and relevant, ACE computes the SVD of their interaction: ⊤=⊤.HJ \;=\;U . The singular values quantify the joint importance of each direction, providing a principled criterion for rank selection. Let rU_r denote the first r columns of U. The algorithmic core is the subspace obtained by projecting these interaction modes back into activation space, ≔span(r⊤).C\, \, span\! (H U_r ). The core’s orthonormal basis ∈ℝD×rQ ^D× r is given by the QR decomposition ⊤r=,H U_r\;=\;QR, and thus, the core projector is defined as ≔⊤P . Note on implementation. Computing ⊤∈ℝN×NKHJ ^N× NK is unnecessary (and inefficient when NK≫DNK D). Instead form the activation covariance ≔⊤∈ℝD×DA H ^D× D and sensitivity matrix ≔⊤∈ℝD×DS J ^D× D. Take square-root factors +ε=⊤A+ =LL and =⊤S= , then compute the SVD of the resulting D×D× D matrix ⊤=⊤,L \;=\;U , which yields the core subspace: span(r). span(LU_r). Causal validation. The core is validated through ablation, with ~ h denoting the activation after intervention: Core-only (to test sufficiency): ~=, h=Ph, Core-removed (to test necessity): ~=−. h=h-Ph. A subspace is deemed sufficient if core-only preserves task performance, and necessary if core-removed reduces performance to approximately chance. The energy-based rank can be refined by finding the minimal r such that keeping only the core maintains baseline accuracy and removing it drops accuracy to near chance. When activity and relevance align. It is worth noting that sometimes ACE reduces to standard PCA – when activity and relevance coincide. This is even expected for simple tasks, when there is no inherent pressure for models to “hide” computations in low-variance subspaces. In more complex models, however, high-variance directions are unlikely to cleanly align with target tasks. Still, the distinction matters even when the subspaces agree: PCA identifies where variance concentrates; ACE identifies where the input–output map flows, by construction and intervention, certifying causal relevance. In other words, PCA is descriptive and statistical, whereas ACE is also causal, licensing downstream treatment and interpretation of the returned subspace and its fitted operator as a dynamical system realizing a causal algorithm. Appendix B Markov Chain Experiment Table A1: Transformer Markov-chain test accuracy: Full Model: no ablations; Core-only: ablating non-core dimensions; Core-removed: ablating core (Methods). Data are plotted in Fig. 1B. Full Model Core-only Core-removed M1 0.748 0.748 0.261 M2 0.748 0.748 0.237 M3 0.748 0.748 0.247 Table A2: Pairwise core geometry and CCA similarity. Projector overlap is the squared Frobenius overlap between core subspaces; angles are principal angles (degrees); CCA lists canonical correlations. Overlap and mean CCA are visualized in Fig. 1C,D. Pair Proj. Overlap Principal Angles CCA M1–M2 0.027 [78, 80, 85] [0.999, 0.999, 0.927] M1–M3 0.031 [76, 80, 85] [0.999, 0.999, 0.949] M2–M3 0.027 [76, 82, 89] [0.999, 0.999, 0.958] Table A3: Eigenvalues λi \ _i \ from operators fit in transformer cores compared with those from the Markov transition probability matrix T (excluding the Perron–Frobenius eigenvalue). Spectral overlap is visualized in Fig. 1E. Markov chain Core1 Core2 Core3 λ1 _1 0.75+0.25i0.75+0.25 i 0.75+0.25i0.75+0.25 i 0.75+0.25i0.75+0.25 i 0.75+0.25i0.75+0.25 i λ2 _2 0.75−0.25i0.75-0.25 i 0.75−0.25i0.75-0.25 i 0.75−0.25i0.75-0.25 i 0.75−0.25i0.75-0.25 i λ3 _3 0.500.50 0.510.51 0.490.49 0.480.48 Three single-layer transformers (dmodel=64d_model=64, dff=256d_f=256, |V|=4|V|=4) with causal attention masking were trained with independent random seeds on next-token prediction for sequences generated by a four-state Markov chain. The Markov chain transition probability matrix, ≔(αβ000αβ000αβ00α), pmatrixα&β&0&0\\ 0&α&β&0\\ 0&0&α&β\\ β&0&0&α pmatrix, was instantiated with α=0.75α=0.75 and β=0.25β=0.25, yielding eigenvalues (spectrum) λ∈1,0.75+0.25i,0.75−0.25i,0.5λ∈\1,0.75+0.25 i,0.75-0.25 i,0.5\, and has stationary distribution =[0.25,0.25,0.25,0.25] π= [0.25,0.25,0.25,0.25 ]. Training used AdamW with learning rate 10−310^-3 and no weight decay for 40 epochs on 3,000 sequences of length 32 generated by T, with batch size 64. Trained model performance is compared against two baselines: Chance: max(), ( π), Bayes-optimal: ∑iπimaxjTij, _i _i _jT_ij, Chance accuracy reflects always predicting the most common token; Bayes-optimal accuracy reflects the best possible one-step prediction given the stochastic nature of the chain. Algorithmic cores were extracted using a 99.9% rank energy threshold without ablation-refinement, H was computed for all test activations, and J was defined by the target function f()≔logits()f(h) (h). Fitting dynamics. Hidden state (mean-centered) sequences were projected into core coordinates t=⊤tz_t=Q h_t and a linear operator was fit by least squares to predict next-step dynamics, t+1≈t.z_t+1 _t. The spectrum of A was used to characterize the learned dynamics. When comparing fitted operators in the core to ground truth, the Perron–Frobenius eigenvalue λ=1λ=1 of T (corresponding to the stationary distribution) is excluded, as it reflects normalization. To calibrate, core operator fits were compared against an oracle ceiling for next-token prediction: Roracle2≔1−1|V| 1⊤(⊘),R_ oracle^2 1- 1|V|\,1 \! (m ), where ≔diag()(⊙(−))⊤,≔⊙(−),m ( π) (T (1-T) ) 1, π ( 1- π), and ⊙ and ⊘ denote elementwise product and division. Appendix C Modular Addition Experiment Three two-layer transformers (dmodel=128d_model=128, dff=512d_f=512, |V|=53|V|=53) were trained on a+b≡c(mod53)a+b≡ c 53. The dataset consists of all 532=280953^2=2809 input pairs, split evenly into train and test sets with a fixed random seed. Input sequences are [a,b][a,b] with target [b,c][b,c]. Training used AdamW with learning rate 10−310^-3, batch size 512, and weight decay ω=1ω=1. Models were trained for 2×1042× 10^4 epochs, with core extraction performed every 100 epochs. The grokking epoch was defined as the first analysis time point at which all three models achieved perfect test accuracy, which occurred at epoch 800. To study the effect of continued weight decay after grokking, at epoch 900, transformers were “branched” – duplicated and split into two regimes – where weight decay was either maintained at ω=1ω=1 or disabled (ω→0ω→ 0) for the remainder of training. For core extraction, H was computed over all test-set activations, and J was estimated using 64 Jacobian samples, defined by the target function f()≔logits()f(h) (h). Core rank was selected first via the 99% energy threshold, and then refined with ablations to ensure causal importance. For operator fitting, centroids ¯c r_c were computed as the centered mean core activation over all test examples with answer token c. A linear shift operator A satisfying, ¯(c+1)mod53≈¯c r_(c+1) 53 r_c was fit by ridge-regularized least squares after dimensionality reduction with SVD. Generalization was evaluated by holding out cycle transitions rather than examples: the 53 answer classes were split into disjoint calibrate/evaluate sets by selecting a contiguous block of classes for the evaluate set, and the fit was performed only on transitions c→c+1c\!→\!c+1 whose endpoints both lie in the calibration class set; evaluation used only transitions whose endpoints both lie in the evaluate class set and fit is denoted as Rh2R_h^2. For descriptive fits, R2R^2 is reported without holding out transitions or ridge-regularization. To summarize spectral structure, eigenvalues of A with magnitude close to 1 were identified as rotational modes, and each such mode was assigned a frequency bin by rounding its angle to the nearest integer multiple of 2π/532π/53. Because complex-conjugate eigenvalue pairs correspond to the same oscillation up to direction, bins k and 53−k53-k were mapped to the same bin. This implies a maximum of ⌊53/2⌋+1=27 53/2 +1=27 distinct bins: one k=0k=0 bin and 26 nonzero oscillatory bins. Mode count is defined as the number of occupied nonzero bins, and derives from operators fit without holding out transitions, since the goal is descriptive characterization rather than generalization evaluation. Appendix D Subject–Verb Agreement Experiment Table A4: Similarity of core coordinates (from Fig. 5B) across six language models. Spearman’s ρ measures rank correlation. Pearson’s correlation r r measures linear relatedness; its magnitude equals CCA for one dimension. Pair Spearman’s ρ Pearson’s r r (CCA1D) (CCA_1D ) Qwen2.5 × LLaMA-3.1 0.921 0.924 Qwen2.5 × Gemma-2 0.934 0.943 Qwen2.5 × GPT-2 Large 0.844 0.857 Qwen2.5 × GPT-2 Medium 0.829 0.842 Qwen2.5 × GPT-2 Small 0.701 0.757 LLaMA-3.1 × Gemma-2 0.888 0.880 LLaMA-3.1 × GPT-2 Large 0.760 0.752 LLaMA-3.1 × GPT-2 Medium 0.727 0.718 LLaMA-3.1 × GPT-2 Small 0.585 0.616 Gemma-2 × GPT-2 Large 0.893 0.909 Gemma-2 × GPT-2 Medium 0.894 0.911 Gemma-2 × GPT-2 Small 0.790 0.846 GPT-2 Large × GPT-2 Medium 0.923 0.951 GPT-2 Large × GPT-2 Small 0.878 0.924 GPT-2 Medium × GPT-2 Small 0.919 0.968 Table A5: A one-dimensional subject–verb agreement core was extracted from each model, despite massive variation in training, model parameterizations, and architectures (model dimension, number of layers). Extracted core size (dcored_ core) is supported by the large spectral gap (ratio of largest two singular value squares σ12/σ22 _1^2/ _2^2). The large spectral gaps indicate that these subspaces are effectively one-dimensional. Model Parameters Layers dmodeld_ model dcored_ core Spectral gap Core location (layer) GPT-2 Small 117 M 12 768 1 4040 11 GPT-2 Medium 345 M 24 1024 1 4444 22 GPT-2 Large 774 M 36 1280 1 2.8×10102.8× 10^10 36 LLaMA-3.1 8 B 32 4096 1 266266 28 Gemma-2 9 B 42 3584 1 133133 42 Qwen2.5 32 B 64 5120 1 531531 62 Table A6: Agreement performance (AUC; 1 = perfect, 0.5 = chance, 0 = inverted) under core ablations across all six models. Core-only preserves agreement (sufficiency), core-removed collapses it below chance (necessity), and core-flipped inverts grammatical number preferences (induces near perfect disagreement). Model Baseline Core-only Core-removed Core-flipped GPT-2 Small 0.911 0.994 0.241 0.038 GPT-2 Medium 0.934 0.997 0.217 0.023 GPT-2 Large 0.921 0.975 0.244 0.021 LLaMA-3.1 0.808 0.918 0.209 0.092 Gemma-2 0.886 0.978 0.102 0.035 Qwen2.5 0.890 0.949 0.213 0.076 GPT-2 Small (117M parameters, 12 layers), Medium (345M, 24 layers), Large (774M, 36 layers), LLaMA-3.1 (8B, 32 layers), Gemma-2 (9B, 42 layers), and Qwen2.5 (32B, 64 layers) were analyzed on subject–verb number agreement. Prompts. A dataset of 1,200 prompts (600 singular, 600 plural) was constructed by combining head nouns (for example, “key”/“keys”, “child”/“children”) with attractor nouns of opposite number (e.g., “cabinets”/“cabinet”) via connectors (“to the”, “near the”, “next to the”, etc.). Five syntactic templates were used: base (“The key to the cabinets”), front-padded (“In this ancient kingdom, the key to the cabinets”), back-padded (“The key to the cabinets in the old kingdom”), existential (“There key near the boxes”), and relative clause (“The key that guards the cabinets”). Half of prompts were prefixed with “In the past,” to vary tense context. The dataset was split evenly into train and test sets. Note: some prompts deliberately employ ungrammatical word order (e.g., “There key near the boxes”) to assess whether the agreement core remains robust to structural violations, forcing the model to resolve agreement based on the head noun rather than positional heuristics. Target function. The number margin was defined on the final-token hidden state h: f()≔(logitare+logitwere)−(logitis+logitwas). f(h) (logit_are+logit_were)-(logit_is+logit_was). Layer sweep. Candidate cores were extracted at each layer and evaluated via ablation. For each model, the layer with maximal flip effect was selected as the core location. Adaptive generation steering. For open-ended generation, a per-token adaptive intervention was applied during autoregressive decoding. Let ∈ℝDq ^D denote the (unit-norm) core axis and μ the mean activation at the intervention layer. The intervention reflects the hidden state h at the last token position through the hyperplane orthogonal to the core axis: ~=−2s[(−)⊤], h=h-2s [(h- μ) q ]q, where s is a per-token steering strength determined adaptively. At each decoding step, three forward passes are performed. First, a gating check: a clean forward pass (with s=0s=0) computes the softmax probability mass on the agreement-relevant verb tokens (is, are, was, were). If this mass falls below a threshold, the token is unlikely to involve an agreement decision and no intervention is applied (s∗=0s =0). Otherwise, the steering strength is calibrated to produce a minimal margin flip. Define the generation margin as m≔log∑v∈are,wereeℓv−log∑v∈is,waseℓvm \! _v∈\are,\,were\e _v- \! _v∈\is,\,was\e _v, where ℓv _v denotes the logit for token v. This logsumexp margin more accurately reflects the probability-space competition between singular and plural verb groups than the linear logit sum used for core extraction, where operating-point independence of the Jacobian is preferred. The calibration proceeds as: (1) the current margin m0m_0 is measured under the clean pass; (2) a small probing perturbation at strength s0s_0 estimates the local gain g=(m1−m0)/s0g=(m_1-m_0)/s_0; (3) the intervention strength is set to s∗=(mtarget−m0)/gs =(m_target-m_0)/g, where mtarget=−sign(m0)εm_target=-sign(m_0)\, targets the minimal margin crossing with buffer ε . An optional cap |s∗|≤scap|s |≤ s_cap prevents extreme extrapolation. This adaptive approach produces grammatical inversions while minimizing collateral disruption to non-agreement tokens. Appendix E Grokking Dynamics E.1 Mathematical Model Let (t)∈ℝμ α(t) ^μ denote the mode coefficients and let ∈ℝμ ψ ^μ be fixed with ψk≔1−cos(2πk/p) _k 1- (2π k/p). Define the (test-relevant) margin m(t)≔⟨(t),⟩m(t) α(t), ψ . Post-memorization, training loss is approximately zero. Updates are driven by the weight decay penalty −ω(t)-ω α(t) and a minimal corrective motion γ(t)γ(t) ψ needed to remain on the zero-loss manifold, plus zero-mean stochasticity ξ(t)ξ(t) (optimizer noise). Direction of ψ ψ. Among all infinitesimal updates Δ α that increase the margin by one unit, the minimum-norm update solves argminΔ‖Δ‖2s.t.⟨Δ,⟩=1. _ α\| α\|_2 .t. α, ψ =1. By the Cauchy–Schwarz inequality, the solution is Δ=/‖22 α= ψ/\| ψ\|_2^2. Thus, the corrective gradient direction is strictly parallel to ψ. Margin dynamics. Differentiating m(t)=⟨(t),⟩m(t)= α(t), ψ and isolating the noise-free deterministic trajectory yields the scalar ODE: m˙(t)=−ωm(t)+γ(t)‖22. m(t)\;=\;-ω\,m(t)\;+\;γ(t)\,\| ψ\|_2^2. Because weight decay is the only systematic drift pulling the network off the margin, the mean corrective force is taken to scale proportionally to maintain zero loss: γ(t)≈cωγ(t)≈ c\,ω for some constant c. Assuming sufficient dimensional capacity (p<dmodelp<d_model), the initial memorized state is unstructured, meaning it carries negligible margin (m(0)≈0m(0)≈ 0). Substituting the corrective force yields a simple linear relaxation equation: m˙(t)=−ωm(t)+cω‖22,m(0)≈0. m(t)=-ω m(t)+c\,ω\| ψ\|_2^2, m(0)≈ 0. Solving this ODE yields the exact continuous-time margin trajectory: m(t)=m∗(1−e−ωt),wherem∗=c‖22=κp.m(t)=m (1-e^-ω t ), m =c\,\| ψ\|_2^2=κ p. Predicting grokking time. Grokking occurs at the first-hitting time τ when the margin reaches the generalization threshold δ. Solving m(τ)=δm(τ)=δ yields the continuous gradient-flow time: τ(p)=−1ωlog(1−δκp).τ(p)=- 1ω \! (1- δκ p ). To map this idealized ODE to discrete training steps, the physical constants are decoupled. The scaling rate becomes: τgrok(p)=−Ωlog(1−pcritp),(pcrit<p<dmodel). _grok(p)=- \! (1- p_critp ), (p_crit<p<d_model). Here, pcrit≔δ/κp_crit δ/κ is the architectural constant, defining the absolute capacity floor limit independent of the optimizer. Conversely, Ω∝(ηω)−1 (ηω)^-1 is the optimizer constant, an empirical parameter that captures the characteristic relaxation time while absorbing discrete step-size dynamics, learning rate, momentum, and adaptive preconditioning from the AdamW optimizer. E.2 Grokking Sweeps and Scaling Fits To measure scaling laws for the grokking delay in modular addition a+b(modp)a+b p, one-layer transformers (dmodel=128d_model=128, dff=512d_f=512) were trained on input pairs using AdamW (lr=1e-3). The data, p2p^2 input pairs, were randomly partitioned into a 50/5050/50 train/test split. Memorization (τmem _mem) and generalization (τgen _gen) times were defined as the first optimizer steps at which train and test accuracy reach 0.990.99, respectively. The grokking delay is evaluated as the difference τgrok≔τgen−τmem _grok _gen- _mem. Accuracy was evaluated every step to avoid quantization artifacts. Two sweeps were performed, averaging over 12 random seeds per condition: (1) Weight decay: fixing p=53p=53 and sweeping ω∈0.3,0.5,1,1.5,2,3ω∈\0.3,0.5,1,1.5,2,3\. To simulate standard training stochasticity on a fixed-size dataset, this sweep utilized minibatch gradient descent with a batch size of B=512B=512. (2) Modulus: fixing ω=1ω=1 and sweeping primes p∈31,43,53,61,79,89,101p∈\31,43,53,61,79,89,101\. Because the dataset size grows quadratically with p, this sweep utilized full-batch gradient descent to ensure the empirical hitting time was isolated from dataset-dependent minibatch noise. Scaling exponents for the asymptotic limits were obtained by fitting power laws y=Cxβy=Cx^β via ordinary least squares in log–log space. The macroscopic constants Ω and pcritp_crit were obtained by fitting the exact deterministic ODE solution τ(p)=−Ωlog(1−pcrit/p)τ(p)=- (1-p_crit/p) to the empirical delay using non-linear least squares (scipy.optimize.curve_fit). Goodness-of-fit for all curves is reported by R2R^2.