Linear Transformers Implicitly Discover Unified Numerical Algorithms
Authors: Patrick Lutz
Aditya Gangrade
Hadi Daneshmand
Venkatesh Saligrama
Venue: NeurIPS 2025 poster ()
Keywords: in-context learning
transformers
generalization
matrix completion
Primary area: deep_learning

Abstract: A transformer is merely a stack of learned data–to–data maps—yet those maps can hide rich algorithms. We train a linear, attention-only transformer on millions of masked-block completion tasks: each prompt is a masked low-rank matrix whose missing block may be (i) a scalar prediction target or (ii) an unseen kernel slice for Nyström extrapolation. The model sees only input–output pairs and a mean-squared loss; it is given no normal equations, no handcrafted iterations, and no hint that the tasks are related. Surprisingly, after training, algebraic unrolling reveals the same parameter-free update rule across all three resource regimes (full visibility, bandwidth-limited heads, rank-limited attention). We prove that this rule achieves second-order convergence on full-batch problems, cuts distributed iteration complexity, and remains accurate with compute-limited attention. Thus, a transformer trained solely to patch missing blocks implicitly discovers a unified, resource-adaptive iterative solver spanning prediction, estimation, and Nyström extrapolation—highlighting a powerful capability of in-context learning.

Linear Transformers Implicitly Discover Unified Numerical Algorithms
Patrick Lutz" Boston University plutz@bu. edu
Aditya Gangrade" Boston University gangrade@bu.
Hadi Daneshmand University of Virginia thadi@virginte edu
Venkatesh Saligrama Boston University srvŌbu. edu
Abstract
A transformer is merely stack of learned datatodata maps-yet those maps can hide rich algorithms. We train a linear, attention-only transformer on millions of masked-bite completion tasks: each prompt is masked low-rank matrix whose missing block may be (i) a scalar prediction target or (ii) unseen kernel slice for Nyström extrapolation. The model sees only input-outpil output pairs and a mqeanes loss; it is given normal equations, no handcrafted iterations, and no hint that the tasks are related. Surprisingly, after training, algebraic unrolling reveals the same parameter-free update rule across all three resource regimes (full visibility, bandvedth heads, rank-limited attention). We prove that this rule achieves second-order convergence iteration complexity, and remains accurate with compute-limited attention. Thus, a transformer trained solely to patch missing blocks implicitly discovers unified, sssource-adaptiv iterative solver spanning prediction, estimation, and Nyström extrapolationhihhighhnng powerful capability of in-context learning.
full-batch problems, cuts distributed
1
Introduction
Models trained on next-token prediction achieve strong performance across various NLP tasks, including question answering, summarization, and translation (Radford et al. 2019). This multitask capability suggests intriguing possibility: within structured mathematical contexts transformers might implicitly learn generalizable numerical algorithms solely through next-token prediction.
Next-token prediction naturally be viewed form of matrix completioninferring missing entries by exploiting dependencies observed data. Prior work extensively explores matrix completion under computational constraints, such distributed data access Jordan al. (2019), limited communication bandwidth Ma and Chen (2020), and low-rank recovery from partial observations Canddè and Recht (2009); y Candès and Tao 2010); Davenport and Romberg (2016). This raises natural question: how do transformers implicitly handle such computational constraints to perform multitask learning?
Can a neural network invent numerical algorithm simply by learning fill in missing data?
Transformers excel at in-context learning (ICL), adapting to new tasks from a short prompt of ex- (2020); Xie al. (2022). and ICL serves as a simple subdomain investigate the representational and learning properties of transformers. Recent studies suggest that for trans- formers trained on ICL tasks, computations encoded by the weights resemble gradient-based methods Akyürek et al. (2023); Von Oswald et al. (2023); Ahn et al. (2023), typically in single-task scenarios. However, this analogy limited: grudient descent operates explicitly in parameter space, whereas transformers perform data-lo-date trunsformations without direct access or updates to
amples [Brwnn
"These authors contributed equally to this work; code is available on gggthuu..com
Research performed under NSF Tripods Foundations of Data Science Institute Grant Boston University.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

parameters. Additionally, much of the exploratory focus in such work has been on the complexity of the regression task, while keeping the computational regime fixed. Thus, it remains unclear whether transformers implicitly develop distinct, unified numerical algorithms suited to diverse tasks and resource constraints.
Masked-block Completion and Architectural Masks: To probe this question, we train a linear transformer on masked-block completion tasks that hide one block of low-rank matrix mirroring the classical Nyström completion problem. The transformer's 's task is to infer these missing blocks based solely on observable data and ureanssquared error objectivewithout explicit guidance about the underlying parameters or relationships between tasks. We impose three distinct architectural visibility constraints on transformers, each reflecting practical computational limitations common in large-scale optimization: centralized (full visibility). distributed (restricted communication), and computation- caaputation-immitee (restricted complexity via low-dimensional attention). Fach regime employs the same underlying transformer architecture, differing only minimally in their attention masks.
Remarkably, despite training independently under these distinct computational constraints, we find that transformers implicitly uncover the same concise, two-line iterative update rule. This unified algorithm, termed EAGLE, emerges consistently across con- straints, exhibiting strong theoretical and empirical properties: it achieves second-order convergence on ball-batch problems, matches classical methods in centralized settings, significantly reduces communication complexity in distributed settings and remains accurate under compute constraints.
Emergence of Unified Algorithm:
Summary of Contributions. The main contributions of this work are: Unified masked-blatt block completion benchmark: A single training framework integrating multiple prediction inference tasks into unified interface.
Resource-aware transformer architectures: Transformer-based models naturally adapt to central- ized, distributed, and computation- -imited environments through simple architectural constraints. Implicil numerical algorithm discovery: Transformers implicitly discover a unified, efficient numer- ical solver, EAGLE, that achieves second-order convergence on full-batch problems and consistently matches or outperforms classical methods across completion, extrapolation, and estimation tasks under varying resource constraints. Our theoretical results uniformly apply across all these tasks, highlighting the broad upplicability of the discovered algorithm.
These findings position transformers trained on block completion powerful tools for uncovering numerical algorithms, offering a promising route toward developing adaptive, data-driven solvers.
2
Related Work
We focus literature most pertinent viewing transformers as fixed data-to-data transforms whose forward pass exposes emergent algorithms.
Implicit algorithu learning. Transformers have been shown to recover _ and 2nd -rder methods for least squares, dual GD for optimal transport, and TD updates for RL (Von Oswald et al. 2023
Ahn et al. Weight-level circuit studies evvrreeengines copy-and-addittin induction heads (Elhage et al. 2021
2023 Daneshmand 2026 Wang et al.) 2025 Fu al. 2024] Giannou et al. 2023).
Olsson et al. 2022). but are still confined to ooken-manipulation algorithms. RNNs and MLPs display related behavior (Siegelmann and Sontag [992] Tong and Pehlevan (024). yet attention yields
depth-aligned, easily inspected updates (Garg et al.) (0022). Most prior work is therefore limited to first-order rules and a single, centralized computational regime.
Task vs. regime. Earlier in-context studies (Akyürek et al. 2023 Bai et al. 2023 Min et al. 2021 vary the regression tasksampling while the hardware regime stays fixed. We take the opposite view: we fix one low-rank matrix-complede task and show that the same transformer discovers rule that adapts automatically as compute, memory, or communication budgets are tightened.
Data-space vs. parameter-space Prior analyses interpret the forward pass of a transformer trained on an ICL task as updating an underlying hidden parameter for some model of the training data relationships. However, in reality, transformer is data-to-data map, with no explicit supervision of underlying model of parameters, which is the viewpoint we adopt. Weight-level circuit studies (Elhage et al. 2021 Olsson et al.) (202)) of trained transformers also take this view, but are focused on tasks like copying or addition; we instead are focused on numerical linear algebra tasks in varying computational settings.
2

---

Ji Draibood saley cralHooo prontareee tor tordene
(8) Pareracire sutg
MIA
Marts And
junca
Figure I: Architectural regimes studied in this work. (A) Computation-liiited. single-head attention- only transformer whose query. key. value projections are restricted to low-rank embedding of dimension == & nI ie., an explicit bottleneck that reduces the per-layer cost from ((n²d) to ((rdd), Unconstruined. The embedding dimension is set to h = n + n'. (B) Distributed: a multi head transformer in which each head operates data stored on separate machine; the heads run local attention and their outputs aggregated. (B') Per-machine routing: detailed view the distributed setting showing how each machine ş forms its local projection Pk: and contributes to the next-layer representation Z
In-context regression limits. Theoretical lenses that view in-context learning as linear regression (Akyürek et al. 2023 Bai et al. 2023 recover gradient-descent dynamics with vk dependence on the condition number. Vlaymmyrov et al. (2024) show that, with suitable parametrization, preconditioning yields second-order dynamics and log convergence. We show that transformer training naturally converges such parametrizations in practice, with this behavior robust to computate and memory constraints. Theoretically, we prove that the recovered EAGLE method enjoys favorable guarantees across all studied settings, and further validate these empirically.
Links to classical numerical algorithms. Structurally, the recovered EAGLE iterations bear intimate relationship the Newton-Schult method for inverting positive matrices (Higham 1997), but extend this to solve non-positive linear systems. In the centralised and sketched settings, this underlying relationship lends EAGLE a second-order convergence rate, with iteration complexity depending only logarithmically on the condition number. This is particularly interesting since usual implementations of Nyström approximations for kernel methods all focus on Krylov-based solvers (Williams and Seeger, 2001 Halko et al. 2011 Gittens and Mahoney 2016). and EAGLE may lend new approaches the same.
Distributed and sketched solvers. Block-CG, communication-avoiding Krylov (Demmel et al. 2013
Hoemmen 2010)) and sketch-and-solve techniques (Clarkson and Woodruff 2017] Woodruff 2014)
dominate practice but remain first-order, and have net communication requirements scaling with the (joint) covariance of the data. We characterise the distributed EAGLE performance in terms of 'data diversity index' a, and find significantly smaller communication (and iteration) complexity when a -1 is much smaller that this joint condition number.
3 Method
We explain how data are generated and encoded, how architectural constraints enforce three resource regimes, and how we extract an explicit numerical solver from trained weights.
Block Completion Setup. We generate a low-rank matrix task by first drawing a matrix of the form
AA
C
X
[gd++d")((+n"),, rank(X) rank(A),
B
D
and then feeding the transformer a prompt matrix Zo which the lower-right grey block masked, ie., set to zero. Given the visible blocks A, B, C, the model must reconstruct D with low & error. This single prompt template generalizes Nyström extrapolation and scalar regression (Appx.B). We train on both exact samples and noisy variants obtained by adding Gaussian noise with variance o² E {0, 10 2). so the low-rank assumption is approximate but realistic. Note that the Nyström approximation is the minimum-rank completion that matches observed entries, even if X is full rank.
Data generation. We construct rank-s matrices X E (dd+d")x(utn) by sampling Ri E (dd=2'xx and R2 E [g(++n)xx )xa each with rows drawn i..dd from (0, E), and setting X R2 R, /5. To
3

---

control the difficulty of the problem, we choose Σ to be diagonal with entries En a², where a < 1, inducing strong anisotropy. We use the following parameters throughout: n 18, n' d' 10 and a = 0.7. When applicable, we add Gaussian noise of variance 0.01 to X.
Linear-attention transformer. We employ the linear-attention variant of transformers (Ahn et al. 2023
Von Oswald et al. 2023 Wang et al. 2025). Each layer e applies multi-head attention with
C
[ Ce via the iterative structure B
residual skip connection, to transform Zo=
to Z, :
B 0
Dt
Ze+1
Ze +
Σ KE[1:H]
Attn (Ze),
Attn (Z) = (ZW2, (ZW ( R o M ZW (1)
where Wo, Wx, Wv, Wp Bank are query, key, value, and projection matrices, and the fixed
mask Md [(a+o") 1 0[a+d')xa'] blocks the flow of information from the incomplete (Ce: D₂)
column in Zy towards the visible column. Models are trained to minimize mean-squared error on D; training details are [B.2] A schematic of the architecture appears in Fig. [
Computational Regimes via Architectural Constraints. We study three distinct computational regimesunconstrained (or centralized), computation limited now-dimensional attention), and dis- tributed (restricted memory per machine)each encoded into the architecture through explicit attention and dimensionality constraints. Figure summarizes these architectural regimes.
Unconstruined. No visibility or dimension constraints are imposed; each token attends to all others. We set the embedding dimension k n'.
Computation- Limitedd To emulate memory- or datency-boundee hardware we enforce a low-rank attention constraint: the query and key matrices have embedding size k < n, and the value and projection size r + n' (we use cuts the per-layer cost of recovered algorithms from (ndd) to ((ndr). (See $B.3 for details).
ss n/4;; see Fig. [A). This bottlemeck compresses the input and
Distributed Computation. The prompt Zu distributes X across M machines through a block structure
AA B
C M M1]
Zo [... X
Xa+1
...] R (g+dd')) ((xxt'))))))"))) where X
D
Each head is assigned to a machine μu, and all but the pth block in its query, key matrices are zeroed, enforcing local self-attention. Each μ then 'transmits' information to others through the value, projection matrices, and incoming transmissions are summed to generate Z4+1- The full algebraic form, matrix partitioning and communication setups provided in $B.4
Algorithm extraction. We generate explicit algorithms from a trained transformer by progressively simplifying its weights until concrete update rule emerges. The steps taken are (also see $B.5):
Weight quantization. Cluster coefficients and sparsify by dropping the values below S We then evaluate the performance with simplified weights to confirm no loss in performance.
Matrix Property Tests. Check whether resulting weight matrices random, sparse, or low-rank. Scaling Laws. Identify how the transformer scales weight matrices layer-by-layer.
1: function UPDATE(A, B,C, D)
1: Input: [A5 D query/output (Co, Do)
2: AS,
B BS
2: for i= 0...., L- do
3: p + [|Al||2
3: Coum +- O, Daum 4-
4: < A A- ĀÄT
AÄT ÄST
4: for 1,..., M in parallel do
5: +- B --pphi BA ÄST
5: (Ai B C", D <- UPDATE(A(, Bi, Cr, De)
6. C' C C
6: end for
T: D' e D DATPBA'
7: Ci+1 1 E, C", Das EE, D
8: return (4', (d,B",C,,D"
8: end for
9: end function
9: Return Dr.
(A) Update: Identical for all regimes
(B) Fusion: Lightweight Communications
Algorithm 1: Emergent Algorithm for Global Low-rank Estimation (EAGLE), Left: the numerical kemel UPDATE runs unchanged on every machine. Right: the outer loop calls UPDATE, and in the distributed regime (M>1). it averages the resulting query updates (C") and outputs (D'). M = 1 in the unconstrained and computation-limite settings. In the former, S = I.., while in the latter. S E IR (d+ is composed of random orthogonal rows. The values y, 7 are global constants (see below for their values).

---

In the unconstrained and compute- compute-limited settings, an ar- chitecture with one head per layer suffices achieve com-
layer
Ux104
Dx101
C-Lx104
petitive performance, while in the distributed setting, we
0.00
0.00
0.00
use one head per machine. After these simplifications ev- ery regime collapses to the same two-line, parameter-free update rule (Eq, in g providing a direct algorithmic interpretation of the transformer's in-context computation.
2.36
0.62
0.01
2
5.27
1.40
0.02
3
3.70
1.28
0.13
4
2.16
1.32
0.14
4 Emergent Algorithm
Table 1: Differences across iterations be-
tween transformer and extracted algorithm (squared Frobenius norm divided by size of 2) in the unconstrained (U), distributed (D), and compute-limited (C-L.) settings. Mean times 10* across seeds reported.
One method, three resource regimes. Across all three architectural constraints, the trained 4-layer transformer achieves 1050 times reduction in test MSE (Fig. B a). Moreover, algebraic interpretation of the learned weights
shows it implements the same two-line scaling transformation across constraints. Algorithm visualises the result, which we call EAGLE (Emergent Algorithm for Global Low-rank Estimation). The blue UPDATE box is identical in all scenarios and across all machines, while the outer grey look captures the information fusion in the distributed setting (i.e., if M > 1). The matrix S in UPDATE is an orthogonal sketching matrix (see below). Table 3 shows that this extracted algorithm reproduces the exact layer-wise activations of the transformer to within 6 10 error, demonstrating its fidelity, We now discuss how this algorithm is derived in the three settings, focusing the noiseless case. However, our findings remain valid under modest data noise (o² 0.01, see SC.2). Moreover, SD shows how Algorithm can be adopted for parameter estimation.
Unconstrained setting. Because token updates are not restricted, the weight patterns are the cleanest interpret here, and will re-appear in the distributed and compute-limited regimes.
Emergent weight structure. During the extraction procedure we progressively pruned attention heads and quantised weights while tracking validation loss. Remarkably, one head per layer is already sufficient: pruning from eight to a single head has little effect on test MSE. After pruning, the weight products Wo.eW K,t and Wv..W.T. P.z collapse to (almost) diagonal matrices.
WoK,e= WQ, diag(a) 0w)
Wvp,: Wv, diag(a a.
where only the three sealars 0g 1.2.3 vary from layer to layer. Figure (left) visualises typical pair; layer-wise statistics across 10 seeds are reported in $C.3 This near- wear--look-ciiggnn structure lets us algebraically reduce the transformers forward step to the update shown in Fig.
unconstrained
distributed
compute-limited
4
4.
10
10-
12
12
16.
16
18
18
12
15
1 2
7
3
12 15
18
2-
2.
4-
8
8
o
10
10
12
12
-1
16
16
3
6
12 15 18
3
6 9 12 1821
3
12 15 18
Figure 2: Block structure of W Wok.e (top) and Wyr. (bottom) learned in the three regimes. Example shown for (n,, d, nd') = (18, 18, 2, 2): distributed run uses M 3 workers with per- worker
6. and all of head,
the [W (WPp)=i,, while omitting (null) non-u blocks W Off-diagonal blocks appear only W p.s. and are identical suggesting structure of messaging. Statistics over 10 seeds are in CC
collated together with each of the three block-wise
corresponding
5

---

From weights to UPDATE. The additive steps in UPDATE are due to the residual connection. Further, the diagonal forms of W W and W., imply that (see $C.4), the attention block acts as
Attne TalaZAAATA. Ar L B Ae o.o
(2)
The form of the iterative update in Alg. with 1,5 I, is then immediate.
Rationale for Weights. Further, the values a. 1.2.3 see two persistent structural effects: throughout, ot0ž x Aell,, and ad a AelA, where ||ll|2 is the spectral norm of Ae (i.e., largest singular value, tuned to the largest such value typically seen in training batch), and 27 the constants and ς z 1.9. Thus, the scale of these weights is determined by || which is used in fixed way by the method. Figure I1 shows that a |ell is constant across layers.
Understanding the method. Note that the transformation + (I pq ) contracts large singular values of A more than small singular values. The net effect of this repeated action Aş is that for large l, all initially nonzero singular values of Ag converge to one another, i.c., the matrix A becomes well-conditioned (see Fig. 3). This aspect of the method is reminiscent of the Newton-Schulz method for matrix inversion (Schulz) 1933/ Ben-Israel] 1965) Higham (0088] Fu et al.) 2024). The overall structure of the iterations
'continuous conditioning update for Ar. The iterations for Cr, D, are reminiscent of gradient descent, adapted to the varying Az. Relation to prior work. Von Oswald et al. (2023) dubbed the same update GD++ in the scalar case (d' n'd==''11 1) and interpreted it as pre-conditioned gradient descent. In our opinion, there are two deficiencies in this viewing of the method. Firstly, gradient based methods implicitly assume that one optimising parameter (i.e., searching for W such that B W A, and imputing D = WCT C). whereas the transformer has not been supervised in a manner that reveals the existence of an underlying parameter. other words, a gradient descent type interpretation not a behavioural (in the sense of Willems (1991)) description of the recovered method. Secondly, our analysis shows that unlike gradient descent, the recovered method is closer to a NewtonSchulz conditioning loop wrapped around direct prediction. For these reasons, we refer to it as the EAGLE method in the rest of the paper. ??? contrasts the two viewpoints in detail.
then be
Distributed setting. Following the unconstrained setting, we train with a single head per machine.
Communication Structure As detailed in $ for each head, the block structure of the (n + π²') x M(n + n²) value-projection matrix W.P. governs what machine sends to other machines. Fig.
(middle) reveals two striking regularities in this matrix.
blocks. Every diagonal block equals diag(a) af Jexactly the same form as in the unconstrained model, with the same pair (a), az) across all machines.
Across-machins blocks. Off-diagonal blocks are diag(O,. a: ): iec only the incomplete columns (Ce, Dz) are transmitted, and with the very same factor oi as the within machine blocks. Resulting algorithm. Since W @K, and the local block of W have the same structure as the un- constrained setting, each machine executes the same local iteration as in the UPDATE method. Further, via the off-diagonal blocks W W.P.g each machine transmits only the O(n' + d')) entries of its update to C₁2, Dr, captured in the C", D' outputs in Alg. I Since attention heads are simply added, ev- ery machine averages its received blocks, leading to C(+1 M ' C. De+1 M M-'E, D
The shared columns are thus identical across all j at all f, recovering Alg. with 5 = I, and > 1.
unconstrained
distributed
compute-limited
cordition (rescaled)
10a
0.050
parameter of
10a
gmae
0.020
a/
0.010
valuee
0.005
10-1
0.002
0.001
10-4
layer
layar
layer
(a) Median test performance
(b) Median condition number
(c) Abstracted parameters
Figure 3: The trained transformer solves matrix completion with unified algorithm all three computational settings. The evolution of key quantities throughout transformer layers illustrate the remarkable similarity between the latent algorithms. Mean across 10 training seeds is reported.

---

Woc
W
orthonormal
norm sp ectrr
1500
froberius
1000
750
0.6
1.0
-0.5 0.0
0.0
0.5
layer
1.0
1.5
2.0
Figure The cippputation-immiiee transformer implements pseudo-random sketching. The figure matches the candidate sketch matrices across all layers against common sketch charac-
Figure 5: The transformer learns to malize the batch-maximum spectral plot Ar observed during training. teristics (randomness, clustered eigenvalues). Left: Distribution reports agaî max [8] [| ||², where the eigenvalues. Right: Distribution confficients.
maximum is taken over each batch.
Communication cost. Under the star topology in Alg. each machine communicates ((d+ d') n/)
floats per round; for scalar prediction (d' = 1) this is the advertised 2d/-float message. Notably, the transformer recovers this sparse pattern without any explicit communication constraint, Exploring how stronger topological bandwidth limits shape the learned algorithm is an open problem.
Computation-limited setting, Here the query, key, value and projection matrices are rank- constrained: Wo,e, WK, & (dd++')xxr and Wv.,, Wp. E R(d-d')x(r+n') with r < n. Similar low-rank constraints commonly arise when attention head dimensions smaller than the trans- former's overall embedding size.
Emergent weight structure. Again, only one head per layer is needed. The (n + n²) (n n') matrices W and W share three universal properties (Fig. 2
Block-diagonal form. Both are block- diack-diagonal; the cross blocks (n x n²) and (n' x n) vanish as in the unconstrained regime. The (n' x n') block W e is a scaled identity, exactly as before. Rank-r top left block. The leading n n block of each matrix has numerical rank r (< n) and the two blocks are similar up to sign when rescaled to spectral norm I (relative difference: 0.28). Random sketch spectrum. The eigenvalues and entry distribution of the top-left block match those S where a random n r orthogonal matrix (Fig. 4
Sketching interpretation. The transformer therefore materializes an orthogonal row sketch S Raxr within the 'top left' blocks of its Wox and Wvp matrices. This sketch acts upon the columns within the Be) blocks of the state Ze. and the output of the attention block is structured as
aa AAS) (AS)ST (A(AS)(AS)CC]
(2)
a (BS) (AS)ST (AS)
where a1, a² are constants depending on ||A||-2 In other words, the update first computes the sketches AS, BS E Rdxr of the 'complete' columns, and then proceeds with the update in the unconstrained case, lifting them back to n x by the terminal S Setting M = and S S; in Algorithm recovers recovvrr the exact layer dynamies.
What the sketch buys. Under the sketch, A, AS E Rdxr is only r columns wide, reducing the per-layer flop count from (n)2d) to O(nrd). Although this sketched update causes the iteration complexity of the method to increase, the overall hope is that the lower per-iteration cost may give better total runtime. This aligns with the original motivation in Vaswani et al.) (2017). where low-rank constraints on Wo.e, WK. Wv.. and Wr. were explicitly introduced for computational efficiency.
Evaluation of EAGLE
Having identified the EAGLE update underlying the transformer weights, we move on to studying how well the extracted algorithm performs in the three regimes, both theoretically and empirically. All ablations are deferred to E while theoretical proofs appear in $F Unless stated otherwise, all numerical evaluations in this section use data sizes
240 and d' 2, condition number x(A) 1 rank(A) 240, average performance across 50 runs is reported.

---

EAGLE
Gradient Descent
Conjugate Gradient
SVD-based
10-1
16-1
LE-1
10
20-10
22-4
1r10
10-1
20
22-44
10
33.15
20
05
25
2.
05
2.6
255
300
255
0.5
05
1.5
(a) rank 30
(b) rank 60
(c) rank 120
(d) rank 240
1.9
187
157
18
20
16
10"
10
16-4
15"1
33-03
3353
33:54
33
310
130
100
199e
1996
(e) K = 100
(1) K
1 1000
(g) = 100000
(h) K 100 000
Figure 6: EAGLE shows second order convergence in the unconstrained regime. extends to low-rank completion and scales logarithmically with the condition number K.
5.1 Unconstrained (centralized) setting
Second-order convergence guarantee. We discuss the theoretical convergence properties of EAGLE in the centralized case. Let (MM) and g(M) denote the largest and smallest positive singular values matrix M and let ni ()) = a (M)/a(() denote the effective condition number.
Theorem 1. For any X, let i. be the Nyström estimate for D and k = m(Ao)- 1f = 1 1/3, 7 then under EAGLE with Ln, M for any > 0, there exists
L = =OOlog log log(s VV [[W.|ee|[CC(]))
such that Ve 2 L. D. |F S where W. BA(AAT)* the Nyström parameter B.
Mechanism. We show Thm. in FF.1 by arguing that each iteration shrinks Kę := k(Ac) by least constant factor, and that once K2 2, then Kę - 1 decays supergeomettrialll. The error is controlled by developing telescoping series with terms decaying with (N₂ - 1). The log log(s
dependence in T., i.e., the quadratic rate, is related to the quadratic decay in K2 - 1, which in turn arises since the EAGLE iterations for Ag hear strong relationship to the classical Newton-Schie -Shuulz method (Higham) 2008
Ch. 5,7). We note that a simple Dr De- < stopping rule suffices in practice.
Positioning classics. Direct inversion via QR-decomposition once, independent of k. Krylov methods such as the Conjugate Gradient method (CG, Hestenes and Stiefel 1952 and gradient descent (GD) need x. logg and k² loogE=1> iterations and ((dd)) time per-iteration to compute matrixvector products. By contrast, EAGLE achieves quadratic convergence with iteration counts scaling only with log(k), albeit with min(ne 2n, d)) cost per iteration.
SVD costs O ndd, n)n)
Empirical protocol. To evaluate the empirical performance of the recovered algorithm, we benchmark EAGLE against a SVD-based solver (torch. wwn.llinllgii. 1stsq), 2 the Conjugate Gradient method and gradient descent (GD) on synthetic A A Rxn (see Appx. E for details). In particular, we study
1. Error vs. wall-clock time. Visually confirms the second-order convergence predicted by Theorem k-sweep. Time to reach = 10 20 for x {102, ...0........ 10²); visualises the log vs. v gap. Rank-deficient check. Time to reach e 10-20 highlights robustness to rank-deficiency.
Figure 6
presents (1)(3;; further experiments including extended size sweeps are in $F.3
Take-away. The transformer--ettraceed update converges quadratically with only a log penalty- while retaining O(min perr))prtittration cost. This positions EAGLE in a previously unexplored region of the speedaccuracy trade-off: an iterative solver that matches SVD-based performance up to log factor.
100x fewer iterations than CG
5.2 Distributed setting
Diversity drives the rate. Besides the local condition numbers =*=)(AA), convergence depends on how distinct the column spaces on different workers are. We capture this via:

---

EAGLE
Distributed Gradient Descent
237
107
UPP
20
130 20
(a)
(b) 0.83
(c) 0.35
(d) & 0.01
Figure 7: In the distributed setting, EAGLE shows significant improvements in iteration count. The dotted line indicates the iteration where the reaches the theory suggest, the convergence depends linearly on a.
Definition 1. (Diversity index) Let Pa be the orthogonal projection onto the range of A", The diver- sity index of (AM) J== is the smallest positive eigenvalue of the average projection 1. Ep=1 P
The diversity index a is large when the projectors cover their joint range uniformly. Conversely, is small when the ranges are highly overlapping yet enlarge the joint span, leaving "innovation" directions weakly covered. The diversity index a captures the convergence scale of distributed EAGLE via the following result shown in $F2
Theorem 2. Let Emax max, R²", In the noiseless case, i.e. when EW, such that B DJ W. [A C7. then under EAGLE with 7j = 1/3,7 S= In M N. for any & > 0, there exists
= 1=00log(Mmas + a- log( loggVV"
such that ve 2.
Mechanism. After log Kimx iterations, every local A has near-unit singular spectrum. Subsequent progress is controlled by the condition number of the average energy matrix Es M- -, Aj A which tends to a Second-order convergence re-emerges in practice when a = 1. Let us note that a-1 is always smaller than k(Eg) and may be much smaller than the same.
Baselines. Our main interest is in comparing with distributed baselines that operate with O(d) units of communication per round-any more, and we may simply share the data across all machines into one central server. For this reason, QR/SVD decompositions and the Conjugate Gradient method are unavailable, which leaves gradient descent (GD) as the main competitor. GD and EAGLE have identical communication costs, transmitting O((d + d') d')) floats per iteration. However, the iteration complexity of GD scales with k< <(E)), which is always greater than a and may be much larger. As scenarios.
a result, EAGLE has a strong advantage practical strongly communication-line
Empirical protocol. We vary the two nate-detemmningg param eters M and a, and report error against iteration in Figs.
EAGLE
Distributed Gradient Descent
10a
1. M-sweep. 111,3, 5, 8 workers among which a fixed 1000, Kmax 103 R 1).
of data is split (d
104
2. a-sweep. Four synthetic datasets, Kmax 103, M
error
3; cross-machine left-subspace alignment is tuned yielding varying centered around means ā E (1, 0.83, 0.35, 0.01).
10a
10-1
Take-away. Distributed EAGLE generally converges at first order and its iteration count grows linearly with a (Fig. 7) Further, for controlled o, the resulting iteration complexity is independent of the number of workers M (Fig.B) [). This confirms our theoretical result. In contrast, gradient based methods require x more rounds. Experiments on best- of-both hybrid (EAGLE until A stabilizes, then GD) are detailed E.
10
100
iteration
Figure 8: Iteration complexity of dis- trihuted EAGLE is independent of the number of workers M.
5.3 Compatation-limite sketching
Sketch-aware iterations. We now turn to iterations driven by the sketched columns Á ST. A, Bi S. B via the iid. orthogonal sketches St Baxr (r < n). The isotropicity

---

of the sketch ensures that the singular of A is conditioned in roughly the same as the centralized setting ( up slowdown of r/ factor.
Baselines Protocol. The standard competitor is Stochastic Gradient Descent (SGD): GD run on an r-column sketch. It shares the n/r spectral slowdown but keeps its native x² dependence. For efficiency, our implementation of the resource-constraile EAGLE uses random row sampling matrices as S. We vary the rank, the single rafe-controllinng parameter
1. r-rank sweep. We vary the rank r E (n/8, n/4, n/2,1i n in Figure 2. Per-iteration cost. The per-iteration wall-clock time decreases with sketch size: 21, 16, and milliseconds for 8 = 240, 120, 60 and 30, respectively.
Take-away. Iteration count grows linearly with n/r. while time per iteration falls up 7x. Against SGD, the method reaches 1=10-2 roughly 2.5x faster for n/4, demonstrating that sketching pluy the second-order update can lead to significant runtime improvements.
5.4 Evaluation on non-synthetic data
We evaluate the unconstrained EAGLE variant on seven low-rank matrices from the SuiteSparse Matrix Collection Davis and Hu (201T) spanning 1100 million coefficients. For each matrix of rank r, and columns permuted such that the top-right r block A is full-rank. For compatability with existing CUDA solvers, we reduce the block matrices B, D to a single row, yielding a setting akin to the linear system W. B.
We benchmark EAGLE against three GPU-accelerated baselines implemented in cupyx.scipy.. wwwp:ysare.....re...... linalg: an SVD solver, conjugate gradients on the normal equa- tions (CG; Hestenes and Stiefel 1952). and LSMR (Fong and Saunders (011). We run each iterative solver until residual of 10 B is reached. Table 2 lists wall-clock times in seconds. First, we observe that EAGLE consistently outperforms the SVD-based solver despite its nominal log overhead: attribute this to its heavy of efficient, massively parallel matrixmatrix multiplication. Second, EAGLE is competitive with sparse iterative methods- beating LSMR on every matrix and CG on all but one-onooeee highlighting EAGLEs practical strength for large, ill-conditioned linear systems. $E.4 corroborates both results for varying residual tolerances.
Task
EAGLE
SVD
CG
LSMR
Maragal_ 4/Maragal
0.625
0.849
1.55
2.34
Baragal_S/Marrgal__
0.632
1.40
1.75
1.40
Mes7aros/pf2177
3.43
4.23
1.75
3.51
HB/dwt_100
0.640
0.902
4.65
9.73
HB/bcsstm13
0.697
0.955
3.17
1.16
Pricbel/162bit
0.847
2.35
16.7
57.1
Schulthess/N_pii
0.675
1.15
1.27
1.96
average rank
1.14
2.29
2.86
3.71
Table 2: Runtime in seconds of different numerical solvers on real-world matrix data. EAGLE is strong competitor for solving linear systems.
Conclusions
We view transformers as fixed data-to-data transforms whose forward pass exposes emergent algo- rithms. We fix one low-rank compiixooo task and explore various regimes. Our proposed rule, EAGLE is transformer-indus distributed and sketching regimes, achieving second-order convergence with only log x (where k is the condition number), a -l (distributed data diversity) n/r (sketch) slow-downs. Empirically, it outperforms Conjugate Gradient and Gradient Descent by 12 orders-of-ne ofders-of-magniuudee in both iteration complexity and net communication, and matches the behaviour of QR-based solvers up to log terms.
method that emerges as unifying algorithm in centralized,
Limitations. Numeric stability beyond <>> 10 is untested. Second, we show one pre-training recipe that yields EAGLE; extensions to non-linear regimes that do so open.
10

---

Acknowledgements
This research was supported by the Army Research Office Grant W911NF2110246 AFRL Grant FA8650-22-C1033 and the National Science Foundation grants CPS-2317079, CCF-2007350, DMS- 2022446, DMS-2022448, and CCF- 1955919...
References
Ahn, K., Cheng, X., Daneshmand, H., and Sra, S. (2023). Transformers learn implement pre- conditioned gradient descent for in-context learning. Advances Neural Information Processing Systems, 36:45614
Akyürek, E., Schuurmans, D., Andreas, J., Ma. T., and Zhou, D. (2023). What learning algorithm is in-context learning? investigations with linear models. International Conference on Learning Representations.
Bai, Y., Chen, F., Wang, H., Xiong, C., and Mei. S. (2023). Transformers statisticians: Provable in-context learning with in-context algorithm selection. Advances in neural information processing systems, 36:57125-57211
Ben-Israel, A. (1965). An iterative method for computing the generalized inverse of an arbitrary matrix. Mathematics of Computation, pages 452-455
Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, J. D., Dhariwal, P., Neelakantan, A., Shyam,
P., Sastry, G., Askell, A., Agarwal, S., Herbert-Voss, A., Krueger, G., Henighan, T., Child, R., Ramesh, A., Ziegler, D., Wu, J., Winter, C., Hesse, C., Chen, M., Sigler, E., Litwin, M., Gray, S., Chess, B., Clark, J., Berner, C., McCandlish, S., Radford, A., Sutskeeer I., and Amodei, D. (2020). Language models are few-shot learners. In Advances Neural Information Processing Systems, volume 33.
Candès, E. J. and Recht, B. (2009). Exact matrix completion via convex optimization. Foundations of Computational Mathematics, 7677777-
Candès, E. J. and Tao, T. (2010). The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5):2053-208) 2080.
Clarkson, K. L. and Woodruff, D. P. (2017). Low-rank approximation and regression in input sparsity time. In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing. STOC 17, pages 81-90.
Daneshmand, H. (2026). In-context learning for discrete optimal transport: Can transformers sort? International Conference Artificial Intelligence and Statistics.
Davenport, M. A. and Romberg, J. (2016). An overview of low-rank matrix recovery from incomplete observations. IEEE Journal of Selected Topics in Signal Processing, 10(4):608-622.
Davis, A. and Hu, Y. (2011). The university of florida sparse matrix collection. ACM Trans. Math. Softe, 38(1).
Demmel, J., Grigori, L.. Hoemmen, M., and Langou, J. (2013). Communicationoopiinaal parallel and sequential krylov subspace methods. SIAM Journal on Scientific Computing,
Elhage, N., Nanda, N., Olsson, C., Henighan, T., Joseph, N., Mann, B., Askell, A., Bai, Y., Chen, A., Conerly, T., DasSarma, N., Drain, D.,, Ganguli, D.,, Hatfield-Dodds. Z., Hernandez, D., Jones, A., Kernion, J., Lovitt, L., Ndousse, K., Amodei, D., Brown, T., Clark, J., Kaplan, J., McCandlish, S., and Olah, C. (2021). A mathematical framework for transformer circuits. Transformer Circuits http://foooooorrriiruuuuusuup/oe h/w.///wwrkiindex..innt
Thread.
Fong, D. C.-L. and Saunders, M. (2011). LSMR: An iterative algorithm for sparse least-squares problems. SIAM Journal on Scientific Computing, 33(5):2950-2971.
Fu, D., Chen, T.-q., Jia, R., and Sharan, V. (2024). Transformers learn to achieve second-order convergence rates for in-context linear regression. Advances Neural Information Processing Systems, 37:98675-98716,

---

Garg, S., Tsipras, D., Liang, P., and Valiant, G. (2022). What can transformers learn in in-context? case study of simple function classes. Advances Neural Information Processing Systems (NeurIPS).
Giannou, A., Rajput, S., Sohn, J.-y, Lee, K., Lee, J. D., and Papailiopoulos, D. (2023). Looped transformers as programmable computers. International Conference on Machine Learning, pages 1 1194-1 PMLR.
Gittens, A. and Mahoney, M. W. (2016). Revisiting the nyström method for improved large-scale machine learning. Journal of Machine Learning Research, 17:1-65.
Halko, N., Martinsson, P., and Tropp. J. A. 1). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217-288
Hestenes, M. R. and Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6):409-436
Higham, N. J. (1997). Stable iterations for the matrix square root. Numerical Algorithms, 15:227-242.
Higham, N. J. (2008). Functions of matrices: theory and computation. SIAM.
Hoemmen, M. (2010). Communication-avoiding krylov subspace methods. University of California, Berkeley.
Jordan, M. I., Lee, J. D., and Yang, Y. (2019). Communication-effiiieenn distributed statistical inference. Journal of the American Statistical Association, 1 114((26):668-6811
Ma, C. and Chen, Y. (2020). Communication-effcienn distributed statistical estimation with optimal convergence rate. In Proceedings of Machine Learning Research (PMLR), volume 108, pages
Min, S., Lewis, M., Zettlemoyer, L.,, and Hajishirzi, H. (2021). Metaicl: Learning to learn in-context. In ACL Findings.
Musco, C. and Musco, C. (2017). Recursive sampling for the nyström method. In Advances in Neural Information Processing Systems (NeurIPS).
Naim, O. and Asher, N. (2025). Two in context learning tasks with complex functions. arXiv preprint arXiv:2502.03503.
Olsson, C., Elhage, N., Nanda, N., Joseph, N., DasSarma, N., Henighan, T., Mann, B., Askell, A., Bai, Y., Chen, A., Conerly, T., Drain, D., Ganguli. D., Hatfield Doddss Z., Hernandez, D., Johnston, S., Jones, A., Kernion, J., Lovitt, L., Ndousse, K., Amodei, D., Brown, T., Clark, J., Kaplan, J., McCandlish, S., and Olah, C. (2022). In-context learning and induction heads.
Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., Sutskever, I., et al. (2019). Language models are unsupervised multitask learners. OpenAl blog, 1(8):9.
Saad, Y. (2003). Iterative Methods for Sparse Linear Systems. SIAM, 2 edition.
Schulz, G. (1933). Iterative berechung der reziproken matrix. ZAMM-Journid of Applied Mathematics and Mechanics/Zeitschrifi für Angewandte Mathematik und Mechanik, 13(1):57-59.
Siegelmann, H. T. and Sontag, E. D. (1992). On the computational power of neural nets. In Proceedings of the fifth annual workshop on Computational learning theory, pages 440-449.
Tong, W. L. and Pehlevan, C. (2024). MLPs learn in-context on regression and classification tasks. arXiv preprint wxXiv22405.16618
Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. (2017). Attention is all you need. In Advances in Neural Information Processing Systems (NIPS).
Vladymyrov, M., Von Oswald, J., Sandler, M., and Ge, (2024). Linear transformers are versatile in-context learners. Advances in Neural Information Processing Systems, 37:48784-48800
12

---

Von Oswald, J., Niklasson, E., Randazzo, E., Sacramento, J., Mordvintsev, A., Zhmoginov, A., and Vladymyrov, M. (2023). Transformers learn in-context by gradient descent. In International Conference on Machine Learning, pages 35151-35174, PMLR.
Wang, J., Blaser, E., Daneshmand, H., and Zhang. S. (2025). Transformers can learn temporal differ- ence methods for in-context reinforcement learning. The Thirteenth International Conference Learning Representations.
Willems, C. (1991). Paradigms and puzzles in the theory of dynamical systems. IEEE Transactions on automatic control, 36(3):259-294.
Williams, C. K. I. and Seeger, M. (2001). Using the nyström method to speed up kernel machines. In Advances in Neural Information Processing Systems, NeurIPS 13, pages 682-688.
Woodruff, D. (2014). Sketching as tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, 10(12):1-157.
Xie, S. M., Raghunathan, A., Liang, P., and Ma, T. (2022). An explanation of in-context learning as implicit bayesian inference. In International Conference on Learning Representations.
13

---

A Related Work
Neural networks as implicit algorithm learners. Large models can emerge classical iterative procedures when trained on synthetic tasks. For transformers this spans gradient descent for least squares (Von Oswald et al. 2023 |Ahn et al., 2023), dual gradient descent for optimal transport (Daneshmand (2026), and lemporaldifferent
planning in RL (Wang et al., (025). RNNs and MLPs exhibit similar behavior (Siegelmann and Sontag 1992 Tong and Pehlevan (002)) though attention yields depth-aligned, cospection-rriendde updates (Garg et al. 2022). Gap. All of these works uncover first-order rules tied to a single regime. Our advance. EAGLE is the first second-order update that generalizes unchanged across centralized, distributed, and sketched settings.
Data transformers vs. parameter updates. Unlike classical solvers that iterate in purameter spuce, a transformers computation is purely data-to-data: weights stay fixed at inference time, and all state lives inside the forward activations. Prior analyses of in-context learning focus mapping this data flow to gradient descent on hidden parameters (Akyürek et al., 2023 Von Oswald et al. )023). EAGLE shows stronger resultan explicit, second-order data-space projector that converges quadratically without ever touching model parameters. This clarifies that algorithmic power can emerge without internal parameter update loop.
In-context learning as regressionrates and limits. Theoretical analyses cast in-context learning as linear regression prompt examples (Akyürek et al. 2023 Bai et al. 2023 Min al. 102TT) Transformers so far reproduce gradient descent or SGD, yielding linear rates and the usual vk condition-numbe dependence. Our advance. We extract emergent NewtonSchulz-style update with quadratic convergence and log & dependence that remains numerically stable under three disparate resource constraints.
Classical second-order solvers, Nyström, and low-rank approximation. NewtonSchulz iterations deliver quadratic convergence for matrix inversion (Higham 1997) but require explicit matrix products and have not been analyzed in noisy, data-dependent regimes. Nyström methods scale kernel learning and matrix completion (Williams and Seeger (001:: Halko et al. 2011 Gittens and Mahoney 2016);
state-of-the-ar variants use QR or Krylov subspace solvers (Saad) 2003 Musco and Musco 2017).
Gap. These algorithms are hand-designed for fixed settings and degrade to first-order (vK) rates when run under communication or memory limits. Our advance. The transformer learns the same NewtonSchull-type projector that attains Nyström accuracy in only log steps and adapts automatically to distributed or sketched execution.
Distributed and sketched solvers. Block-CG, communication-opiini Krylov methods (Demmel at 2013 Hoemmen 2010) and randomized sketch-and-sole techniques (Clarkson and Woodruff 2017 Woodrufi 2014) are the de-facto choices at scale, yet remain first-order and Ii-limeeed. We introduce diversiry index a that predicts our distributed rate and achieves provably fewer rounds when worker subspaces overlap (ox E Vk).
Meta-learning across tasks. Transformers pretrained heterogeneous corpora behave meta learners (Min et al.) 2021 Bai et al. 2023 Naim and Asher (2025); prior work, however, stops at behavioral observation. EAGLE offers mechanistic account: pre-training can imprint a single algorithm that flexes with compute, memory, and communication budgets.
Scope. This survey is necessarily selective; we focus works most relevant to the data-space, second-order perspective. A broader taxonomy of numerical solvers and in-context phenomena is beyond our scope but complementary to the questions addressed here.
Positioning. To our knowledge this is the first demonstration that pre-training alone can induce one second-order, matrix-completion rule that (i) bridges centralized, distributed, and sketch regimes, (ii) admits tight non-asymptotic guarantees, and (iii) out-iterates classical Krylov baselines and cuts communication under typical ML-scale conditioning.
14

---

10
setting
mean square difference
weight
unconstrained
103
Wok.
10a
distributed
Wyp.
compute-iigg limited
-0
10
MSE
10
10
10-3
10-3
10
5000
10000
15000
20000
step
number clusters
Figure 9: Training loss, median across 10 inde- pendent runs.
Figure 10: Prediction difference when sparsify- ing weights in centralized setting. For reference, the black dotted line marks the trained trans- former's mean squared error. Mean across 10 training runs reported.
B Methodological Details
B.1 Nyström Approximation
A
C
The Nyström approximation of a block matrix X
with missing block D the matrix
D
A
C
8,
D.
where
D. = BA(AAT)CC.
Here, Zt is the Moore-Penrose pseudoinverse of Z. Note that this approximation ensures that the rank of the approximate X.. is the same as the rank of A-in this way, the Nystrom approximation reduces both the rank of X, as well as the Frobenius seen from the fact that the approximation just amounts to joint execution of ordinary least squares. Indeed, if we treat the ith row of B, B' as 'response variable' and the jth column of C, cf query point, then
of the error. The latter property can be
(..;; = B' B'A(AT
preciscly the estimate when regreeses the data (A, B3) onto e². For this reason, we may also interpret the method as computing the estimation parameters
W, = BA(AA E ddd
and then computing the estimate
b D.WWC
much in ordinary linear regression. This value W. will appear later in our theoretical analyses.
We note that the approximation above is defined whether the matrix X noise-corrupted or not. In the presence of noise, it inherits many of the statistical properties of lincar regression when assume that the noise is restricted to the 'response variables' B D. As a special case, if the rank of X is equal to that of [A C), then in fact holds that [ D] W, [A C1 i.e., the Nyström approximation is exact. We will sometimes refer to this as the 'noiseless' case.
B.2 Training specifications
The transformer is trained using the Adam optimizer with constant learning rate of 0.001 and a batch size of 1024 for 20,000 iterations. To stabilize training, gradient clipping with a 2-norm threshold of 0.1 is applied at each step. At every iteration, a new batch is independently sampled
15

---

as detailed in Section 3 ensuring that the model never sees the same data twice during training or inference. The block masking pattern remains fixed throughout. We observe occasional spikes in the training loss (Figure ② and interpret those as occasional failures of the learned solver to converge, typically occurring when ||||2 attains untypically large values. We recall that, under our data sampling procedure, ||/ is random and unbounded, allowing for such outlier cases.
B.3 Runtime attention map
Recall that the dominant computational cost transformer layer arises from the attention mecha- nism:
(ZWo(ZWx) 2 @ @ WWWWW.
Aun(2) 2)
where ((d++d")(((+)) (x(+n') the input matrix and Wq. WK, Wv, Wp E B2nk are learned projec- tion matrices. Assuming that the size of the submatrix D remains constant, we have d' = ((1) and n' = O1))
The computation of ZWq. ZWx, and ZWv requires (nndk) time. The inner product and masking operation (2Wq) (ZWx) k)) @ M incurs cost of (dd)), and the subsequent multiplication with ZWywp requires an additional O + ndk). Thus, the total time complexity of the attention computation (dd2 + ndk).
For the unconstrained case = n, this yields a cost of O(n?² + dn²). Therefore, selecting k < significantly reduce computational overhead.
Of course, in our recovered algorithms, we can "precompute' WoW and Wy WTT, and these are further scaled versions of block-identity matrices, and the computation of ZWQ etc. can be avoided. This reduces the cost to k)),), which for the unconstrained setting of & works out to OlPn).
In the sketched versions of the iterations, these matrices become nontrivial again, and their multipli- cation must be incorporated into the accounting of the costs of the recovered algorithm. Of course, in this regime, k = < n, so the domainting cost is ((drr).
B.4 Details setup for distributed setting
We design the transformer emulate a distributed algorithm by constraining the attention mechanism. Specifically, the query, key, and value matrices are restricted to data local to single machine, and due to symmetry across machines, a single set of transformation is shared across all attention heads. At every layer, the input data is partitioned as:
A B
C D
Zp = [... X
X
.... E R (d+ M(n+nn) where X
e [1: M]
with blocks A E Raxd and Ba IRR''x xd distributed across machines, while the shared blocks C' E Raxdd and D E IR" xd are replicated identically in each X". This setup reflects common scenarios where the missing block is significantly smaller than the full matrix A, often remaining constant in size (e.g., d' =n' in regression), making the distributed partitioning both efficient and natural.
The projection matrix after attention is unconstrained, enabling unrestricted communication between heads and implicitly modeling a fully connected communication topology. This architectural design encourages local, maahine-specifii computation in the attention mechanism. while allowing for global coordination in the projection step. Empirically, we find that the model leverages this flexibility efficiently: although capable of learning full data sharing, it consistently limits itself to communication scaling with DI) bits.
We find that transposing the per-machine data X is necessary in the distributed setup. That is, we b xa+1. ...] and train it on the mean-squared
give the transformer inputs Z = [...
error to DT.
16

---

layer
Ux104
Dx10a
C-Lx100
romm
spacital Trobenlus
0.00
0.00
0.00
2.36
0.62
0.01
2
5.27
1.40
0.02
3.70
1.28
0.13
4
2.16
1.32
0.14
Table 3: Differences across iterations between transformer and extracted algorithm (squared Frobenius norm divided by size of Z) in the un- constrained (U), distributed (D), and compute- limited (C-L.) settings. Mean times 10a across 10 seeds is reported.
Figure 11: The transformer learns to normal- ize the batch-maximum spectral norm of A observed during training. The plot reports af max || where the maximum is taken over each batch.
B.5 Details on algorithm extraction
We derive closed-form iterative updates from the trained transformer weights WoK. and Wvr.t. This involves thresholding and clustering the weight matrices in both the unconstrained and distributed settings:
The sparsification threshold is set dynamically as 1.5 ( #elemets,,,nssss. accounting for varying scales across weight matrices. The constant 1.5 selected empirically and is not further tuned.
The number of quantization levels is determined by analyzing the mean prediction error induced by weight sparsification (see Figure I0). To maintain performance comparable to the original transformer, we cluster W into 2 values and Wvp. into 3.
In the compute-limise setting, a sketch matrix extracted and its structural properties analyzed (see Section 4. Parameter choices for the update steps are abstracted through empirically observed scaling laws. Specifically, we identify a') 0' Rs 1/ //A/ as supported by results the main text. Furthermore, we approximate a1 & 1.92 03, where the resulting squared prediction error is 3.5 10-4 (averaged over 10 training runs), remaining below baseline error approximately 1 x 10-3
C Emergent Algorithm
C.1 Numerical evidence for faithful extraction
Table the difference between the transformer embeddings and abstracted iterations for every transformer layer. The largest error of x 10-4 is still substantially smaller than the transformer' ''s overall predictive performance after 4 layers (around 1 x 10-3. This shows that the abstracted iteration truthfully reflects the numerical procedure implemented by the transformer.
Figure 11 shows the product a. a max (1.45) I| for difference choices of norms. In this expres- sion, the maximum is taken over training bach of size B. The plot reveals that the transformer scales the constant a with the maximum spectral norm of A, expected at layer during training.
C.2 Transformer behaviour under data noise
In this section, we demonstrate that the emergent behavior identified Section remains robust under moderate data noise with variance の² = 0.01.
Figure 13 compares the learned weights representative transformer trained on noiseless data and on data corrupted with noise. Additionally, Figure reports key diagnostic quantities across layers. In both analyses, the observed behavior closely aligns with that of the noiseless case, indicating stability of the learned algorithm under modest perturbations.
17

---

Layer
Layer
distributed
vanilla
distributed
compute-
022
0.50
8-
10
10
0.5
12
15
-0.0
02
-04
16.
16.
16.
3
15
16
3
22151921
15 18
3 6
12
1211 1821
Layer
Layer
vanilla
distributed
compute-limited
varilla
distributed
compute-limited
6
B
8-
10
10-
10
12
12
18
6
0.1
12
9 12 18
6 12
2.
10.
10
10.
12.
1921
Figure 12: Transformer weights by layer, averaged over 10 training The top panels display the weight products WoK,,, and the bottom panels show Wyp.t. Prior to averaging, sign ambiguities were resolved by aligning the dominant coefficients: large entries in Wok, were set to be positive, and those in WyP. to be negative.
noise (a² 0)
modest noise (o² 0.01)
distributed
computo-lmitor
compute-limite
10
8
10
10
10
16
12
18
18.
15 18
12
12 15
0
12
0.0
30
10-
10.
-0.5
16
18
18
Figure 13: Comparison of learned model weights without noise vs. with modest data noise on representative training run. The top row shows W QQ. and the bottom row shows WvP.
18

---

C.3 Transformer behaviour across training runs
We show that the emergent behavior identified in Section 3 is robust across training runs with varying sources of randomness, including initialization and training data. Figure 12 displays transformer weights averaged over 10 independent runs, reproducing the characteristic patterns reported in Section Additionally, the sketching sub-matrix in the compute-limited setting is, on average, close to the identity, further supporting out interpretation as a random orthonormal sketch.
C.4 From weights updates
We derive the blockwise update rule implied by the abstracted weight parametrization
[an I
0 and WVP.
a
WoK,-
0
and show that recovers the update structure presented in Section 4 for the unconstrained setting. Recall that the layer representation takes the block form
A Ce Be Dr
Ze
and single-head transformer layer performs the update
Ze+1 = Ze + (Z.WoKeZe @ Me) Ze Wvp.C.
(4)
Substituting the block structure into the attention term yields
AeA] (.A
ZeWok.ZZ =
where the entries marked- are not needed due masking. Applying the attention mask,
TeA 0
Z.Wo o Mr
0n B
0
Further, we compute
TazA agCi az agDi
ZWvP
Multiplying these terms and adding to Ze. we obtain the update:
[1122A4A] A B1A A a1a3 Co
a103AAA Ci
Z
Ze
This clean separation of updates across the blocks A,, Br, Ct, and De posteriori motivates the choice of block decomposition of Z; used in notation throughout the model layers.
unconstrained
distributed
compute-llmited
cordition (rexcaled)
10a
10a
gomma
7
volue
10-1
x
2x10-
10-2
(a) Median test performance
(b) Median condition number
(c) Abstracted parameters
The trained transformer solves matrix completion (modest noise. o² 0.01) with unified algorithm over all three computational settings. The evolution of key quantities throughout the transformer layers illustrate the remarkable similarity between the latent algorithms. Mean across 10 training seeds reported.
Figure 14:
19

---

EAGLE
Distributed Gradient Descent
accelerated FAGLE
10a
10
10-5 5
10-10
10
error
error
10
1
10 20
10-2
10
10-3
0.0
0.5
1.0
1.5
time
time
Figure 15: Convergence of EAGLE adapted the estimation problem: unconstrained setting, with parameters d 240, n' d' 2, <(A) 100. Mean over 50 runs is reported.
Figure 16: The accelerated EAGLE for dis- tributed computing largely outperforms EAGLE in runtime. In this data generation, the data di- versity a = 0.1 is relatively low. Mean over 50 runs is reported, <(A) = 100
d 240, n' d' 2.
C.5 Interpreting the emergent algorithm
Prior work interprets EAGLE gradient descent with preconditioning (Von Oswald al. 2023 Ahn
et al. 2023 Vladymyrov et al. 2024). In particular, these works frame the method as comprising two distinet phases: preconditioning step involving only the updates to and B (this is achieved by setting = 0), followed by a single gradient descent step using the data (Ar, B,, Co). This perspective treats the conditioning phase as auxiliary and external to the actual optimization step.
However, this interpretation does not accurately reflect the behavior of trained transformer.
The model does not explicitly separate conditioning and update phases; instead, it interleaves them continuously throughout the computation. Rather than implementing classical preconditioning followed by gradient descent, the learned procedure functions as continuous conditioning mechanism that shapes the dynamics of the optimization process at every layer.
Moreover, this prior viewpoint underemphasizes the internal structure and complexity of the condi- tioning dynamics themselves. fact, the iteration bears close resemblance the Nowton-Schulle method for matrix inversion, suggesting fundamentally different mechanism at play.
To make this connection explicit, consider the unconstrained version of the iteration. Assume ||A||2 1, and set 1, y = 3, as in the statement of Theorem | Then, re-normalize A₃ = Ae/ll Aell at every iteration. The update to the iterates A in this setup can then be expressed in the compact form
Å 1 (31 -31- -31-
which is exactly the Newton-Sehhe iteration for computing the matrix sign function, as presented in (Higham 2008, Section 5.3). This reformulation reveals that the driving force behind the iteration is more accurately characterized as implicit matrix inversion rather than traditional gradient descent updates.
D Parameter Estimation
Beyond the standard prediction task in which the hidden block D is to be recovered, one may also be interested estimating the matrix W BA(AAT)) Nyström formulationwhere D reduces to scalar 1 x blockthis corresponds to solving the classical least-squares problem. While it is straightforward to recover W by setting C (Algorithm 1 so that the prediction becomes B BA(AA = BA(AAT) W*, this approach is inefficient terms of memory and runtime. It thus natural ask whether the estimation task can be achieved more efficiently, with reduced overhead.
the least-squares variant of our
EAGLE
20

---

(a) W
(b) WB.C
Wwr.
Figure 17: Learned transformer weights in the estimation setting. The transformer architecture is modified and bias term introduced: Ze+1 Ze + ((zWoQka QK, W Me) WvP.8 where WB E gaxd a rank-2 bias term whose first n 1 rows are identical by design.
To investigate this question, we consider an adapted transformer architecture trained on matrix completion problem posed as
AA ?
B 0
Z =
where the transformer is tasked with filling the missing block "?" with W BA(AAT)1. We find that minor architectural modification is necessary to achieve competitive performance on this estimation task. Specifically, applied when updating the first n tokens, and the second applied when updating the last tokens.
introduce two biases to the query transformation. The first bias is
We train transformer on the least-squares setting (i.e., with D E RR l and extract the learned weights, shown in Figure [7] This leads modified estimation version of EAGLE, defined by the following iterative updates:
A4+1 A- - Ac,
Be+1 B -2- B A
We+1 = We YPe We BB)AT.
These updates can be adapted to all computational settings considered the main body. Importantly, this iteration is theoretically equivalent to the original EAGLE algorithm in the sense that W₂C DE. where We evolves according to the update rule above and D& evolves according to EAGLE. This equivalence can be formally shown by expanding the recursive updates and applying a telescoping argument, as in the proof in $ F
Figure 1S provides empirical evidence that the proposed estimation iteration (with the same parameter setup Pe |/(3| Ae (12) 1 1, 7 3) converges as expected and recovers the same second order convergence EAGLE.
E Further Details on Evaluation of EAGLE
E.1 Details on the implementation of baselines
We provide additional implementation details for the baseline methods used in comparison with EAGLE, All three implementations compute an approximation of X* BA' and return D R X*C.
Exact closed-form baseline (mp.linal,, lstsq) This native numpy impementation of the least- squares solver is base on highly optimized LAPACK routine. It determines solution to the over/under- determined linear system minx X|X BI The closed-form minimiser is X* BA(AAT)) returns D = X.C.
21

---

Gradient Descent. We minimize the reconstruction loss
f(X) - XA - B
1=1
Parameters. The learning rate is set to ne m/ (A)), ensuring monotonic decrease of f. Spectral norms are estimated via two power iterations, which incur negligible memory and runtime overhead. To avoid instabilities caused by near-singular matrices, we add fixed ridge penalty 10 when low-rank, modifying the gradient as + A When is full rank, we use = O.
Iteration and Variants. The method extends naturally distributed stochastic data acess variants. Each worker computes either exact column-sketched gradients (Áž the data worker i with subsampled columns) which are aggregated on a central node, followed by a synchronous gradient (A)))) This yields the following iterative updates:
step (?μ- = 1/ max (
Gk = Σ ==1
B + AXx
(gradient + ridge)
T]k Gk
X = Xx
(synchronous update)
m
Training stops when the iterate X reaches predefined tolerance after fixed maximum number of iterations.
Conjugate Gradient Method. We implement the Conjugate Gradient (CG) method to solve the normal equations arising in Nyström estimation. The objective is to find X minimizing the empirical loss 2 Ei= || |XA- B. We define the Gram matrix G AAT and the right-hand side B BA
and solve the system XG = B using CG.
Iterations. The standard CG updates are as follows:
R4 = B - XRG,
k 0: Po = Ro, LTo (Ro, Ro)
0g = rTk/ P&G),
X = =X, + tapp
Rk+1 R4 P&G,
M+++1 (R R
Bk 1T+41/TTK
Pk+1 Ret1 + Pr
The algorithm halts when |X&C' - Xx-1C1Ipp falls below pre-specified tolerance or after fixed number of iterations.
Limitations. The CG method is not directly suited for distributed or stochastic settings. In distributed environments, computing G and B requires sharing all local data across machines, incurring com- munication costs that exceed the desired (dd) complexity. In stochastic settings, CG requires a consistent system matrix every iteration guarantee convergence, which is not available when data is sampled independently at each step. Empirically, CG under stochastic updates exhibits unstable and non-convergent behavior.
E.2 Accelerated EAGLE for distributed setting
In the distributed setting, convergence is impeded by low data diversity, characterized by a < 1. Recall that the continuous conditioning performed by EAGLE operates exclusively on per-machine data. As a result, the condition number of the global matrix A is lower bounded by a a-1 even if the data each individual machine is perfectly conditioned.
Once this bottleneck is reached, further updates to 4he most computationally intensive component of EAGGEccasse to be effective. This motivates the question of whether runtime can be reduced by terminating updates to A and 8" once per-machine conditioning has converged.
To this end, we implement an accelerated version of EAGLE for the distributed setting. At each iteration, evaluate the criterion min,, || At A I| < 10-10 and halt updates to A and 8
22

---

for all jl once this threshold is satisfied. This condition ensures that per-minn per-machine data is sufficiently well-conditioned.
Figure T6 reports the runtime performance of this accelerated variant. As expected, the modified version achieves faster execution compared standard EAGLE. We note that the number of iterations required remains unchanged; the acceleration only affects wall-clock efficiency. The theoretical guaranteed of the accelerated EAGLE remain unchanged as we ensure that x(AA)) (to machine precision) before freezing the updates to A and B
E.3 Ablations for EAGLE
Data distribution. We perform ablations on various real-world data settings. As in the main body, we place ourselvels the setting n d 240, n' d' 2 and r 240.
SVD construction. Sample Asj. Bi, orthogonal and set = AEB /VT with E diagonal with entries log-uniform in [le 2, 1. This is the reference distribution in the simulation part.
Gaussian. Sample Aij, Bij (N( 1) and set Z AB /vF. Gaussian data is widely made modeling assumption, and this is the distribution the transformer was trained on.
Student-ty factors. Sample Ag,. Bij t///v/(w 2) for = 4 and set Z = AB /VF. Heavy- tailed entries create sporadic outliers, checking that the solver is stable beyond aub-Gaussian data.
Correlated Gaussian factors. Draw each column of and B from (0, Σ) with Σn أ (p = 0.8). Strong row/column correlations are typical in spatio-temporal panels and challenge methods assuming independent observations.
Sparse Rademacher factors. With sparsity p = 0.1, take Au = vF where Eιy ~ Bernoullipp). Sij {+1} (analogous for B), then Z AB ABT/VF. Entry-wise sparsity yields highly coherent singular vectorsa worst-case regime for many theoretical guarantees.
Block / factors,, Assign each of the d rows one of k: clusters (k 5), {0, 1jaxk is one-hot, and B e E. axk holds cluster centroids; structure seen in recommender systems, violating global incoherence yet preserving low rank.
ABT /vk. Produces piece-wise-cre wise-constant
Larger systems. We increase the size of the data matrix X to d 1920 and n' d' = 2: results are shown in Figure 20 These experiments further reveal that the Conjugate Gradient method does not converge reliably when applied low-rank matrices. As result, its performance is omitted from the main results in the low-rank setting.
Data noise. To assess robustness noise, add moderate ii.d. Gaussian noise to the matrix X and evaluate EAGLE in all three computational settings; results are presented in Figure 2T We observe that EAGLE exhibits the same form of implicit regularization during early iterations that well-documented for gradient descent. Across all settings, EAGLE with early stopping achieves performance comparable to gradient descent, recovering regularized solutions that more stable than the Nyström solution Additionally, we note that the addition of Gaussian noise to X inadvertently improves the conditioning of the matrix A, effectively reducing x(A).
Sensitivity to step sizes. To evaluate the sensitivity of EAGLE to the choice of step sizes, we con- duct series of controlled experiments in the unconstrained setting, using a small-scale configuration with n = d = 15, n' = d' = 2, and rank 15.
EAGLE
Stochastic Gradient Descent
Conjugate Grradient
QR-based
162
112
167
10-
10*
35"
32*
10-
10-1
0.2
6.5
50
0.55
025
Ure
(a) 30
(b) r 60
(c) r = 120
(d) r = 240
Figure 18: the stochastic data access regimea variant of the compute- limited seltingEAGLE outperforms SGD in wall-clock time for sketch sizes r 2 n//..
23

---

EAGLE
Gradient Descent
Conjugate Gradient
SVD-based
10"
10
10-4
10-17
10-22
10"4
10-24
10:41
10"14
10
(a) SVD construction
(b) gaussian
student
10a
10
10-1
10-3
10
10
10-44
10-42
10-12
10
10
10"
10
10
10
10
10
10
10"
157
10
10"
10"H
10
10"
10"
Umc
time
tim
(d) correlated gaussian
(e) rademacher
(f) block
Figure 19: Performance with different data distribution, d 240, ni' d' 2, rank 240.
Gradient
Conjugate Gradient
10.5
10
10-20
10-4
10-25
10--5
10-20
10-3.
10-2
10
10
20
time
time
(a) rank 240
(b) rank 480
100
10a
10-5
30-30
10-44
20-25
10-45
10
10
10 m
rank 940
(d) rank 1920
Figure 20: Performance on larger problems, n d 1920, ni' = d' = 2.
EAGLE
Gradient
Conjugate Gradient
10a
10.
30*
10/
10'
301
101
204
10
16
10a
10
16a
Iteration
Figure 21: Performance with i.i.d. Gaussian noise (standard deviation o = 0.01) in the unconstrained (left), distributed (center) and stochastic data access (right) settings.

---

10
10
10-
10-11
10-14
10-17
10-19
10-21
10
19-29
10-2
3
2
3
scaling
ratio
Figure 22: Best performance of EAGLE 200 iterations as the step sizes are varied. Left: y and are both scaled by a common factor. Right: y is scaled by Vratio while 1j scaled by 1 //vratio such that //n is scaled by ratio.
Task
Size
Nonzeros
rank
conditioning
Maragal_ Maragal_4/Maragaa
4/Maragal
1964 x 1034
1.3 %
801
1.38
107
4654 3320
0.6 %
2147
1.22
106
Meszaros/pf2177
9728 x 10178
0.031%
9662
2.19 x 103
HB/dwt_1007
1007 x 1007
0.85%
1000
2.97 x 10a
HBBbbssttm3
2003 2003
0.53%
1241
2.70
104
Priehel/162bit
3606 x 3597
0.29%
3460
1.68 x 10a
Schulthess/N _
3625 x 3923
0.057%
2048
2.47 x 102
Table Relevant characteristics for each SuiteSparse Matrix Collection task.
First, scale both step sizes and by common factor, shown in Figure 22 (left). This setup models scenarios in which the estimate of the scaling parameter = -1/l| is inaccurate. The results indicate that scaling 7 and y) jointly factor in the range [0.2, 2] has minimal effect performance.
Second, vary the ratio 7/n test the algorithms robustness to discrepancies between these two parameters, as illustrated in Figure (2 (right). find that EAGLE remains stable over a wide range of ratios; specifically, varying 7/1 within [0.2, 2| has little impact on overall performance.
E.4 Evaluation SuiteSparse Matrix Collection data
Table summarizes the characteristics of all evaluation tasks. Moreover, figure 23 gives the runtime of common linear solvers on real-world data as function of the target residual. We observe three trends. First, EAGLE outperforms the SVD-based solver when very high accuracies are required. Second, as the theory suggests, the runtime of gradient-based methods out scales that of EAGLE when very high accuracies required. Finally, we observe tradeoff where gradient based solvers can outperform EAGLE for large, sparse and well-conditioned systems such as Mcszaros/pf2117.
F Theoretical Analyses
F1 Analysis of the Unconstrained or Centralized Iteration
We begin by analysing the behaviour of Algorithm in the 'centralized' setting of S = M For simplicity, let us set ne n||-Ae|l 1|= and N'
M / 1.
will incorporate the specific
25

---

Maragal 4/Maragal
Maragal 5/Maragal
Meszaros/p12177
20
10
time
5
5
2
0.5
10
10-7
10
10
10
10
HB/dwt
1007
HB/bcsstm 13
Priebel/162bit
20
time
5
1
2
0.5
0.5
10
10
10
10
10-4
10
Schulthess/N pid
epsilon
epsilon
time
10-4
10
epsilon
meth od EAGLE
Conjugate Gradlent
LSMR
SVD-based
Figure 23: Runtime of linear solvers on real-world data function of the target residual c BB
structure of 71 later. Then write the iterations
Ac+1 = Ar(I - neAd Ae),
(5)
Bi+1 = Be(1 - neAt As),
C++1 (1--AAA) )Ce,,
D&+1 De + 7 Ce.
Of course, we set Ag. Bo, Co to be the observed blocks of the input matrix X, and Do Useful Definitions. We will analyze these iterations through the follow two objects Definition 2. For l 2 0, define the "signal energy' E,, and the "signal correlation' V2 as
& = A.A E ddd and V = B:A E dd''d xd
We further note the critical relationship that
Bg W.
Vo W.Eo,
where W, is the Nyström parameter for the data matrix (see B.I An immediate consequence of the iterations above is the following dynamics.
Lemma 3. Under the iterations of (5) we have the following dynamics
E111 = E7(1 - n2E2)
(6)
V/+1 = Ve(T - ET²
C4+1 (I C
D+11 = D: NoVe
where I is the d x d identity matrix.
26

---

Proof.
The De, Ce iterations follow immediately by definition. For the former iterations, first observe
that
Ve+1 = B1+144+1 By - n.A A1) neA
= =Be(1 - -TTAA 7A A(AT)
=BelL - ) (1-
(1 ( NEE)))
where the first few inequalities arise by multiplying A? from the right, and then factoring it out by the left, and the final is by definition. Similarly,
Eet1 1 A AelL mA Ac)(1 Ac)
A (1 mAA,
= Ee(T - M)EE))
where we reuse this transfer of A from the right to the left, and then use the definition of Eg.
Note from the above expression that for any l' > l, Eg is a polynomial in &g. As a result, & and Eg commute for all pairs ()).
Telescoping the Error. For the sake of conciseness the further expressions, we further introduce the following notation.
Definition 3. Define Mt == Il - mTEi) and N == Ilieel - mEi), where we use the convention that
II = U2U- ... U-e-1 t<
and Mo = No I.
By Lemma we have & Mr and V Vo Further, Ct = NeCo
This leads the following estimate of the value of De.
Lemma 4. Let W, be the Nyström regression parameter for the data (A, B). For any l 2 0, we have
D2 = W. No)Co.
Proof. Since Dn 0, we have
Dj = Σ
nuc = V4
Σ
ETMiiN Co
Lee
Ice
EncaMiN) Co lse
= W.
Σ
Co W,
W.
Σ
HENNN
Co.
I<t
where we used the definition Mp. But now, using the definition of N₃. notice that for any 4
NoT N, (I EE)
N411 N₁ -nE NENN No N
Consequently, we can write
C(Ni-Ni+1) C. = WI 1<B
D. = W.
"Note that the order of multiplication is flipped here, which follows notationally since we defined the transpose N above. Of course, all the & commute, so this is not very important in this case, but this subtle distinction will matter more subsequent analysis.
27

---

Recall that D. = W.Cu, meaning that the error D D, is W,NCCo Thus, this lemma captures basic fact: the size of the matrix N₂ decays, the output D3 gets closer to the target output D.. Thus, our main focus is characterize the decay of this matrix. We shall show this by arguing that the matrices & quickly become well-conditioned through the course of the iterations, and consequently Ne decays quickly.
Before proceeding, note that in general, Eμ may have zero eigenvalues, which are left unchanged by the main dynamics (see below), meaning that I I-N; Ng may even asymptotically have large eigenvalue corresponding to vectors in the kernel of Eg. However, we observe that W. Vo and thus any energy in Co that lies in the kernel of &n is irrelevant to the prediction D.. As such, it is equivalent for us to analyse the two-norm of the matrix N₂ == (T No) where Po projects onto the kernel of Eg. Instead of this notational complication, we will henceforth just assume that Eo is full rank (i.e., all of its eigenvalues are nonzero), and mention where changes need to be made to handle the general case.
Conditioning of Ece Let us then begin with some notation: for a positive (semi-) symmetric matrix M, we let Xi (M) denote its ith largest eigenvalue, X(MM) denote its largest eigenvalue, and X(M) denote its smallest (nonzero) eigenvalue. Notice that X(-) X²)),, and X(M = (MM)., where is the rank of M (we will work with full rank M. but be cognizant of rank sensitivity). Let v² (M) denote the corresponding eigenvector. We further introduce the notation
X = (E)), Xe = -(Ee), and X = A(EE).
Our first observation is that the iterations leave the eigenstructure of & undisturbed. Lemma 5. Suppose that is an eigenvector of &g. Then it is also an eigenvector of for any E >> 0. Further, if v E kcr(EE) then ker(Ee).
Proof. We have E = c.. But note that & is some polynomial in Ep, i..., for some finite number of coefficients On, Egv En anEg = (Enzo and")o ee), and the claim follows. Crucially, observe that since & McEo; we have agg = 0, and so if Egv O, then Egr
well.
Since the cigenstructure of E is preserved, we can begin to study the behaviour of its eigenvalues in simplified way. The critical relationship for our analysis is the following iterative structure on the largest and smallest (nonzero) eigenvalues of &
Lemma 6. Suppose that AY, m) S Xx /3. Then for all E
Ae+1 = - 7=>) Xe
and det1 = (1-1)-d))2
Proof. Suppose u², 22 are two eigenvectors of Eg. with positive eigenvalues pa, No. Since w u² are also eigenvectors of &,, denote the eigenvalues for the latter as pe, ve. Note that these remain nonnegative. Indeed, by multiplying the dynamics for Ex+1 by, say, u we have
Pet1u² = Ec+1u² Ee(1 - E = (1 mefe)?
and similarly for ve. We will first inductively show that μN > vo + ue > vs for all f. other words, not only are the eigenvectors of En stable under the iterations, but also the ordering of the eigenvalues.
We thus note that the claim follows directly from this result. Indeed, let e be the eigenvector corresponding to X Then we see that it remains eigenvector corresponding to X₃ for all f. But then using our observation above with u² = v', we get the result for λ441- A similar argument works for \ Note that nothing per se demands that these eigenvalues have unit multiplicity, and the argument is completely insensitive to this.
To see the ordering claim on μt, Mz, then, we observe that
pe+1 Vt+1 = - (1-eee4)210 pM (1 - 17212)2 2
= [we 11e) - 2n - v²² + 1²² (1z 133 = (Me vz)) (1 - 2n u + 1z) + (1 + v2 + Hevi))
28

---

Now, if the term multiplying (ne - v2)) is nonnegative, then the claim follows. To see when this occurs, let us set Te = niie and pull A within the brackets. We are left with an expression the form
(1 2n( + y) + 1 + + xy)), where (z,y) (0, For 7) s 2/3, the minimum over this range occurs when y 1, and takes the value
4n + 31
is a triviality to show that this function is nonnegative for 1 S 1/3, and so we are done.
With this in hand, we will argue that 2 - X at quadratic convergence rate after an initial burn-in period. Of course, this is equivalent to saying that Ke := Xe/ Ar + Notice that this Ke is precisely the condition number of & (which in turn is the square of the condition number of A₃ mentioned in our discussion in $5)
Lemma 7. Define Ki. = Xe/ and 02 =Ke- 1. If' S X /3, then Oc+1 S Og 0xexp(-55x2e/
Further, if me 1/3, then 02+1 30//44
Proof. For the sake of succinctness, let write ne 1,, N, d.0 0, and 0c+1 0+. Then, using the previous lemma, we can write
2
1- mnd
K8+1 = 0+ +1= (0+1)
1 ndd
But notice that
1 1-11 7/
7(1 - 1/(0 / + 1))
θng
=1 17
1 n/
1- 7d
1 - / +1)
0 1 T
Thus,
Oyd
0A1p2²
04 + 1=+1=000 101=00011
11)
1-2 (0 ++ -17)
+ (0 1 n/
n
On
= 0 0
2
n/
(0 -1-1111 n7,)
But the bracketed term is nonincreasing in nर (in the negative term, the numerator increases and the denominator decreases with this value), and so
Tn
0
50
4
5
2
0+1
+
30
3
Thus, we have
5/ 5/3meAe
5me
O1+1 < 02 Or
O.+1
S 02 exp
< Og exp(-
n
1 m/de
On the other hand (going back to the notation that suppresses E), if we set TTĀ L, then since
1- L
n (0+1-L)) ()/())
we find by doing the long multiplication that
(0+ + 1) (0 + 1) -1)2
(0+1-- 1-
0+ A(1 - 4L 3L2) + 02 - 6L + 3L2) + 03 L)2 (0 L)2
Setting L = πe 1/3, the linear term in 0 vanishes, and we end up with
0 = 0 / 403/9 02
+
302
(30 + 2)² S
(0 2/3)/
4
29

---

The statement of this Lemma has shown two decay modes of 0₂ (Ke 1). Firstly, as long as "e S >xe /3, this quantity decays least at a linear rate. Further, ne is set to exactly X /3, then once 8g dips below 1, which occurs after about 3 log(8o) iterations, bound in to recover a stronger convergence rate. Together, this yields the following convergence behaviour for Kg.
exploit the quadratic
Lemma 8. Let Kig = de for all E, ne = 1/3, then 1 S for all
2 L 2 + (3355) log(<o)] + [log2 (log(4/33)))
Proof. First using the geometric decay in Lemma 7 we know that 0g S (Ko 1) exp(-5(e- 1)/3). Set Eo + log 1). Then for E 2 Eu. we have 6 < 1/c. Further, for iterations beyond this fo. the supergeometric decay 8/41 - 36} /4 in Lemma implies that
30 (30:/4)² = Te+1 I 2Tc,
where Te := 1og(33//4). Thus,
Te S 2 To -22--6 = 0g S exp
Setting this < gives the claim.
We again note that this Ke is the square of the condition number of A₁ however, this distinction is easily accommodated due to the logarithmic dependence on it the only change is that the 3/5 above increases 6/5.
Back to No- The conditioning of & yields a direct control on the behaviour of Ne. as encapsulated below. We note that if Eg is not full rank, the statement holds for Nr = (I - Po) where Po projects onto the kernel of Eg.
Lemma 9. For all e, Nelll S exp((- Ei-eteendi Further, if m1 X /3 and v x then
||Nella < for all é 2 L + 1, where
L = 2 + [3/5 log( + [loga 4/3=))].
Proof. Since N441 (1 ( we immediately have
||Nel12 S II i<E
711112.
But note that ||4 EEill2 = (1 11)) Thus, we immediately have ||Nell s II(1 -nd) S
exp(Eiee nd)) Further, recall from Lemma that if me = x /3 for all f, then d 2 Te/(1 + for all f L. But then, if me x
have
- ned 5 -1
</
yielding the claimed bound on ||-eel2.
Finishing the Error Control. With all the pieces in place, conclude the main argument.
Proof of Theorem By Lemma 4 we have Dj D. = W,N&Co and consequently,
D - D.IlF I || |WWNe|l=/|Cooll= S No WWW//////l
But, since W.N, Rd xd the rank about is most d'. Finally, [W.Nell ||Nella||W.11. Thus, the claim follows as soon as ||Nell2 < Ill for which we may invoke Lemmag
Again, recall that for the purposes of error, if & were not full rank, then we could instead replace Ne by N₂ in all statements above, and thus the claim also extends to this situation.
30

---

Comment on Rates 7. Going back to the notation of Algorithm |
we reparametrise "/ = Pe X || and Ne pe/3. is worth discussing briefly how we may go about estimating this
Pe.
simple observation in this setting that if assume that we know A = ||Eollo = ||All begin with, then it is simple matter to compute the subsequent values of Xe, since if we set Te T as per the above, this is directly computed iteratively via
Te-1 Te(11 1/3)a = 4Xe/9.
In fact, under this assumption, we may avoid further numerical stability issues by first recaling to have 2-norm 1, and subsequently simply rescaling the matrices up Ae, Be by 3/2 after each update to maintain the invariant Xτ Âh-the behaviour of N₂ remains the same, and convergence rates are only driven by the conditioning of Aş. So, equivalently, the question we must concern ourselves with finding ||A||2- Of course, this is quite cheap, since we can compute this simply by power iteration, which involves only matrix-vector products rather than matrix-matriy products. Nevertheless, note that power iteration is only linearly convergent (i.e., the number of iterations needed to find a e-approximation of the top eigenvalue is ((log(I/s)) so in full generality, this procedure would destroy the second order convergence.
In our simulations of $ we indeed implemented these iterations by carrying out power iteration. Note that practically, then, the method essentially retains its second order behaviour despite this approximation. One aspect of this lies in the fact that even if we underestimate the top eigenvalue by small amount, and so undertune 7 by a slight amount, the second order bound of Lemma decays gracefully, and so retains practical resilience. Of course, characterising exactly how loose this can be requires more precise analysis, which we leave for future work.
F.2 Distributed Analysis
Beginning again with Algorithm 1 with S = Tn, M 1, we need to analyse the iterations
V4, A A/ 1/ Art A/D
B, B - 1/ Ar Ab).
C++1 = Ct
1 Σ
re Ap TT.
Ce.
m
P
D1+1 De+ m
Σ BB A Ce.
We note that here we have set "/ to be constant across all machines, and nominally it is not pegged to ||4//| which could vary across machines. However, we will analyze this method under the normalization assumption that
F/, 11,6112 = 1.
In this setting, we will show that the choices v/ = 11,4/11a2/3 are also the same for each machine, and then set Ye 3rf" (which in turn is also the same for every machine). Thus, in the setting we analyze, the iterations above are faitful the structure of Algorithm
Note, of course, that this distributed data computes the Nyström approximation
A
C
B
D
where C, D are as in the iteration, while
A [A]
A2
Am) B
[Ba
B²
Bm]
We will work in the noise-free regime, wherein there exists matrix W, such that
[B DJ
W. [A C),
31

---

and C E (olumn-span(AA) Of course, this W. also equals the Nyström parameter for this matrix (B..). The second condition ensures that the rank of this matrix the same that of A, and is needed to ensure that we can actually infer W. in the directions constituted by the columns of C (without which one cannot recover D)).
Per-machine and Global Signal Energies. We begin with defining the energy and correlation matrices in analogy the previous section. As before, Ajj A", Di 1 Co C, Do 0.
Definition 4. We define the per-machine objects
& = A AT and VV = B AMT
and the per- per-machine values
X = (² = (
where r(u) is the rank of e Further, we define the global signal energy
E, == E
Note that by our assumption, ||E0"l22 1 for all μ.. Further, the behaviour of E V/ is identical to that in the centralised case, i.e.,
Lemma 10. For all p. e
Et = E - - and Vitt -f T//
As a consequence, for every l, μu, it holds that
Vit =W.E.
Finally. for any L
C+11 = (1-
Proof. The iterations for E 14 Vi follow identically to the proof of Lemma For the second claim, notice that we have
B W.
y W.E.
by multiplying by T on both sides. Then the claimed invariance follows inductively.
Finally, noting that E Ali Ap: we immediately have
C++1 I- Σ E -
In analogy with the centralized case, this leads us to define the following matrices. Definition 5. We define
M
II - 8/ and N = II nEE,),
ICF
where again [[ interpreted multiplication from the right, and M² No Id-
Succinetly, then, can write E = Mr and C₁ N2Co-
Telescoping the Error. This notation allows us set up the following analogue of Lemma Lemma 11. Under the distributed iterations without noise, it holds for all é that
Dj W.(T No)Co.
Proof. Observe that
DI
Σ Dr
2 Σ
Ct.
m
1<t
32

---

Now, V/" = W.Ef" and so
Σ
VI W,
Σ
e = W.EL.
m
m
11
/
Thus, we have
De W. Σ
nELNi Co.
lef
But
N41 = (I- N - nEiN -
and so the term in the brackets telescopes to No - Ne.
Conditioning. Of course, again, D = W.C W.(T)Co. So, to gain error control, we only need to argue (as before) that N [| vanishes quickly with E. As before, this relies strongly on the condition number of Eg. We note, again, that we will simply assume that & full-rank, since rank ranf-defic is rendered moot in this case by the fact that C' lies in the column span of A (and hence, is orthogonal to the kernel of Eo). However, the individual Eg may not be full rank, and so the distinction between X(E") and the smallest eigenvalue of 9d (which is usually 0) should not be forgotten, although it will not matter very much for our expressions.
To begin with, note that since the per-machine A", B" iterations identical to the central case, the corresponding Ed are conditioned at the same quadratic rate we saw previously. We formally state this below.
Lemma 12. For all f, set T/ == (x*) /3, and define Ni² = =x//2
If L(e)= = 2 [3/5 log mmax xO)1 + [log2 (10g(//g(443)) then max ni S +1+
Proof. Apply Lemma g
Of course, &,, as an average, will not have the same conditioning behaviour. In fact, this is strongly captured below.
sensitive the diversity index from Definition ①
Lemma 13. Suppose that N0 is constant across machines, and for all μ. l, n = (x1)-1/0 113)-/33.0 /3. Define
R = X(E)//EE))
1
Ift 2 L(e), then Tie S
a
where L(E) is the expression in Lemma 12 Further, for the same range of for all μ.
X a 1+
X(E
Proof. Firstly, by Weyl's inequality, notice that
1 Σ A(E))).
X(E)) S m
Further, recall that X(g(') for μ, and we set 1/62 = (Xx)) Thus, each of these X are infact identical (and equal to (4/9)²). Thus,
Yu, X(,) < X
Now, let Pa be the projection onto the column space of AH. Then we note that for any vector
a T 8 ()*)) (P")).
Indeed, any component in e orthogonal the column space must lie the kernel of E A# A and we know that this kernel is invariant across the iterations. Further, note that since PPv lies in this
33

---

column space, it is orthogonal to any eigenvector of E with zero eigenvalue, and so we can then conclude that
(P* E (P 2 A P* ull.
But then, if 2 [(e), then for any unit vector v we have (for any j in the last inequalities) that
v E
1
Σ
PH
X
1+
X
1+c
where we used the equality of the Ne: the fact that K/ S +1+5 E, and finally the definition of a. Since is unit vector, conclude that any Rayleigh quotient of & is so lower bound, ergo
X
Vp, A(E) 2 1+E
a
Putting this together with the upper bound on X(Ex), immediately conclude that
R4 1+
At high level, a= is the limiting condition number of the E,s, where deviation from perfect conditioning occurs only due to how the various E energise distinct subspaces for large μ.. We note that this dependence is tight-if all the Ed share the top cigenvector, then the upper bound on A(E,) is exact, while the lowest nonzero eigenvector of & is precisely the direction that achieves the minimum in the definition of the diversity index.
Concluding the argument. With this in hand, we can argue the decay of N₂, and so attain error control.
Recall that D - D. W.N&CO. We again assume that X = for all jr, and /3. We set 7/ (xx)-1 S X(E,)-1 Then, as before, we have (restricted to the appropriate subspace if Eo is not full rank)
Proof of Theorem 2
that v/" = 1
|| T S (1- (1-1/Re) exp(-1/Re)
By Lemma 13 for l Z L(L), we have Rez < 2/a, and so
||&el| < exp(-
1/Re) < T(1))/2) Σ
1=L(1)
which is smaller than 18 if
1
8> L(1 log
Of course, L(1) 3+ 13/5 ~ and setting so small that error follows in the same way as the proof of Theorem concludes the argument.
Comment on Rates. We note that, in the distributed setting, within- machine computation typically much cheaper than across machine communication. As result, the protocol we have analyzed above, wherein each machine begins with the same value of XX" is cheap to follow by using power-iteration at each machine. Given this choice, the value X is simply equal (4/9) for every machine, and setting the learning rate of p /3 for P (9/4)* is trivial, and does not require any communication.
34

---

We also note that setting P in each machine, and then averaging the results of the updates together sufficient, since o constant across all machines. Thus, we recover exactly the structure presented in Algorithm [ with = 1, 1/3. Note, again, that for the sake of stability, instead of working with decaying A, and exploding ne Tt: we can again rescale each A Bg by 3/2 after
updates.
In general, if this normalization is not carried out, the machines may actually their values for y as distinct s/s. The net effect would be that need analyze slightly different version of & that
sensitive to inter-machine variations in yt which in turn would rely on how strongly X varies across machines. We leave the study of such scenarios to future work.
35

---

1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: [NA]
Guidelines:
36

---

The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution
specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility, In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification:
The simulations and transformer training is lightweight and easy to implement.
37

---

Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (attps://nips.co cc/
cablic/guides///
for more details.
While we encourage the release of code and data, understand that this might not be possible, so No is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips cc/ wttpc/pbblccgguuudeeseeeeeeeeeee wtw.a//ddee/oeeeeunssssssssses
stlobmmiisoooooccy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if' applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer:
[Yes]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
38

---

For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: [NA]
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- cration due to laws or regulations their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deploymentt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinfornation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
39

---

The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswithcode com/ datasete has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: [NA]
40

---

Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation,
or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with
human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: [NA]
Guidelines:
41

---

The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy for what should or should not be described.
42


Reviews:
[{"review_id": "VaZJDdYQo6", "reviewer": "Reviewer_BaDn", "summary": "This paper explores the discovery of numerical linear algebra algorithms through training a linear transformer to complete masked matrix data. The authors generate a dataset by sampling matrices and masking parts of them, then train a linear transformer to recover the missing entries. Subsequently, a distillation procedure is applied to extract key components from the learned model. Remarkably, this process yields an emergent algorithm that, in a basic setting, exhibits second-order convergence. Moreover, the distilled algorithm generalizes to other contexts such as distributed computation and matrix sketching, suggesting it may serve as a general framework for algorithm design.", "questions": "Could you clarify the number of transformer layers $L$ used during training? I could not find this detail in the main text or Appendix B. Based on Figure 2, it appears that only 3–4 layers are sufficient—does this mean that the model is trained with just 3–4 layers (which correspond to iterations in the context of numerical algorithms), yet the distilled algorithm generalizes to hundreds of iterations at inference time? If so, this is quite impressive and worth emphasizing.", "limitations": "N/A", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "This work is insightful, original, and clearly presented. It demonstrates strong potential for high impact in the domain of scientific discovery through AI. I did not identify any major concerns with the methodology or results.", "was_revised": true, "final_justification": "The authors have addressed my question, so I will keep my positive rate on this paper.", "rebuttal": "We sincerely thank the reviewer for their exceptionally positive and insightful feedback. We are thrilled that they found our work to be insightful, original, and of high potential impact, and we appreciate their careful reading of our paper.\n\nTo answer your excellent question regarding the number of transformer layers: **you are precisely correct.**\n\nThe transformer was indeed trained with a small number of layers (typically **L=4**), as shown in the layer-wise performance plots like Figure 2. This corresponds to the model performing only 4 \"in-context\" iterations of the learned algorithm.\n\nYour observation that the distilled algorithm, EAGLE, can then be unrolled and iterated for hundreds of iterations at inference time—far beyond what the transformer observed during training—is a key result that highlights the true generalization power of the discovered algorithm.\n\nThank you for pointing out that this is an impressive aspect of our work. Following your suggestion, we will be sure to emphasize this point more clearly in the final version of the paper. We are very grateful for your careful review and encouraging assessment."}, {"review_id": "DIHuA3QMEN", "reviewer": "Reviewer_rz5S", "summary": "The authors explored learning numerical algorithms for low-rank matrix completion with Transformers. The authors train transformers, without explicit algorithmic guidance, on tasks encompassing prediction, estimation, and Nyström extrapolation, imposing three computational regimes (centralized, computation-limited, and distributed) through architectural constraints. Analysis of the trained weights reveals a consistent, parameter-free, two-line iterative update rule applicable across all regimes.  EAGLE theoretically and empirically matches or surpasses classical solvers, automatically adapts to resource limitations, and demonstrates strong convergence properties.  The results highlight the potential for in-context algorithmic learning within transformers.", "questions": "- How does EAGLE's performance scale with increasing matrix size, different data distributions (e.g., sparsity, noise levels), and matrix structures?  What are its limitations in high-dimensional or highly structured scenarios? \n- Given the inherent parallelism in matrix operations, how does EAGLE's performance compare to optimized parallel implementations of classical solvers on multi-core systems?  What is the potential for further parallelization within EAGLE itself? \n- Is there any chance to evaluate EAGLE's capability on well established matrix completion benchmarks (e.g., SuiteSparse Matrix Collection)", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\n- EAGLE is a single, interpretable iterative method learned by the transformer for various matrix completion tasks (prediction, estimation, extrapolation). This unification is remarkable given the differing computational constraints.\n- EAGLE's inherent adaptability to limited computation, memory, or communication resources is a significant strength, suggesting practical applicability in resource-constrained environments.\n- The paper provides a strong foundation through both theoretical analysis and consistent empirical validation across the three computational regimes. \n\n## Weakness\n- **Scalibility** The experiments primarily use small, fixed-size matrices.  The lack of scalability analysis raises concerns about performance and the applicability of EAGLE to large-scale, real-world problems which often necessitate approximate solutions. \n- **Synthetic Data** The reliance on synthetic data, with controlled noise and parameters, limits the generalizability of the findings.  The paper lacks analysis on real-world, potentially noisy and less structured datasets, hindering the assessment of EAGLE's real world feasibility.", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for their thoughtful feedback and their interest in our work. We appreciate the concerns raised about scalability and generalizability, and we are grateful for the opportunity to clarify these points and the central contributions of our paper. Thanks also for suggesting experiments on the *SuiteSparse Matrix Collection*, which we have done below. We hope this response will clearly illuminate the key aspects of our submission.\n\nThe primary discovery of our paper is that a transformer, trained simply on a matrix completion task, can **implicitly learn** a single, sophisticated numerical algorithm (EAGLE) and **automatically adapt** its execution strategy for **centralized, distributed, and computation-limited** hardware regimes. The emergence of this unified, resource-adaptive solver is, we believe, a surprising and powerful demonstration of in-context learning. While we show that EAGLE itself is a strong algorithm, our **core focus** is on the transformer's emergent ability to derive and unify it.\n\n-----\n\n### Weaknesses\n\n**1. Scalability:** *The experiments primarily use small, fixed-size matrices. The lack of scalability analysis raises concerns about performance and the applicability of EAGLE to large-scale, real-world problems...*\n\nWe respectfully point out that our work provides both strong theoretical and empirical evidence for EAGLE's scalability.\n\n  - **Theoretical Guarantees:** Our main theoretical contributions, Theorems 1, 2, and 3, provide explicit convergence guarantees that hold for matrices of **any size** and under **any** distribution (i.e., result is data independent). The convergence rates depend on fundamental matrix properties like condition number ($\\\\kappa$) and data diversity ($\\\\alpha$), not on the matrix dimensions or the distribution themselves. Specifically, EAGLE's logarithmic dependence on condition number ($O(\\\\log \\\\kappa)$) makes it theoretically well-suited for large, ill-conditioned problems where first-order methods struggle.\n\n  - **Existing Large-Scale Experiments:** While the main text reports experiments on matrices up to $240 \\\\times 240$ to clearly illustrate the dynamics, we provide further scalability analysis in the appendix. **Figure 19 (p. 28)** explicitly tests EAGLE on problems up to size **$1920 \\\\times 1920$**, and **Figure 7(h) (p. 9)** demonstrates its robustness on matrices with a condition number of $\\\\kappa = 100,000$. These results confirm the favorable scaling predicted by our theory.\n\n**2. Synthetic Data:** *The reliance on synthetic data, with controlled noise and parameters, limits the generalizability of the findings.*\n\nWe argue that the use of synthetic data is a methodological strength for our study's goals, and our theory ensures the findings generalize.\n\n  - **Controlled Stress-Testing:** Our primary goal is to understand the *mechanism* by which the transformer adapts its learned algorithm. Synthetic data allows us to precisely control the parameters that our theory identifies as critical performance drivers—such as condition number and data structure—enabling a rigorous validation and stress-testing of the emergent algorithm in ways that would be infeasible with uncontrolled real-world data.\n\n  - **Distribution-Agnostic Theory:** Crucially, our theoretical guarantees are **distribution-agnostic**. The convergence of EAGLE depends on the spectral properties of the data matrix, not on the underlying distribution of its entries. To demonstrate this robustness empirically, **Figure 18 (p. 28)** shows that EAGLE performs consistently across **six different data distributions**, including heavy-tailed (Student-t), correlated, and sparse factors.\n\n  - **New Real-World Benchmarks:** Nevertheless, we agree that evaluation on real-world data is valuable. Following your suggestion, we have benchmarked EAGLE on several matrices from the **SuiteSparse Matrix Collection**. The results, presented at the end of this response, confirm EAGLE's strong performance.\n\n-----\n\n### Questions\n\n1.  *How does EAGLE's performance scale with increasing matrix size, different data distributions..., and what are its limitations?*\n\n    As discussed above, our paper analyzes this from two perspectives:\n\n      - **Computational Performance (Speed):** Our theory and experiments show that EAGLE's iteration count scales logarithmically with the condition number ($O(\\\\log \\\\kappa)$) and polynomially with data diversity ($\\\\alpha^{-1}$) or sketch size ($n/r$). This is a significant advantage over first-order methods like Gradient Descent, whose complexity scales with $\\\\kappa$ or $\\\\kappa^2$. This makes EAGLE particularly effective for the high-dimensional, ill-conditioned problems often found in practice.\n\n      - **Statistical Performance (Accuracy):** EAGLE is proven to converge to the optimal Nyström approximation. The suitability of this solution depends on the task; for problems with different underlying structures (e.g., entry-wise sparsity), a different learning target would be needed. Our work's novel contribution is to fix the statistical task to isolate how a transformer adapts its *computational strategy* to different *hardware regimes* —a fundamentally different axis of investigation from prior work that studies adaptation to different statistical tasks.\n\n2.  *How does EAGLE compare to optimized parallel implementations... and what is its potential for parallelization?*\n\n    EAGLE is very well parallelizable, \n\n      - **High-Level Parallelism:** The **distributed version of EAGLE is itself a parallel algorithm** designed for multi-worker systems. As detailed in Sections 3 and 5.2, it performs computation-intensive work locally on each worker and requires only lightweight communication of the partial results for the shared columns ($O((d+d')n')$ floats per round), making it naturally suited for modern parallel hardware.\n\n      - **Low-Level Parallelism:** The core operations within EAGLE's update rule are matrix-matrix products. These operations are highly parallelizable and have been extensively optimized in standard libraries (e.g., BLAS/LAPACK).\n\n    A direct performance comparison with a library like ScaLAPACK would require an equally optimized C++/CUDA implementation of EAGLE, which we believe is an exciting piece of future engineering work but is beyond the scope of this paper's focus on the transformer-based discovery of the algorithm.\n\n3.  *Is there any chance to evaluate EAGLE on... the SuiteSparse Matrix Collection?*\n\n    Yes, this is an excellent suggestion. We have run these experiments.\n\n----\n\n### Experiments on SparseSuite Matrix Collection (SSMC)\nWe benchmark EAGLE against several baselines on matrices from the SSMC. These were chosen from **diverse domains**, and chosen to be **big and sparse**, to handle both the concerns raised in the review. Their condition numbers also vary a lot. In line with the paper's setting, we take these low rank matrices, and generate a matrix completion task in the $A,B,C,D$ format from the paper. We then hide $D$ and run each solver. \n\nFor iterative methods, we say an iterative method converges when the update norm satisfies $\\|D_{t+1}-D_t\\|\\le10^{-6}$ within a maximum wall-clock time of 20 min. We implement EAGLE, CG, GD, QR/SVD, and also stochastic EAGLE ($r = n/5$) and SGD for small batch experiments.  Implementation details for each baseline are in Section E.1 (p. 26). \n\n**Relative Convergence Times** Because convergence times and errors vary greatly across tasks, we report all results relative to NumPy’s native `np.linalg.lstsq` implementation (which internally uses SVD/QR). So the SVD/QR time is set to $1$ in each row.  Each solver is run three times with different random seeds, and we report the median runtime. We write DNC for \"did not converge.\" In these cases, errors were all large. Otherwise, errors all fell below 1 and below $10^{-4}$ for EAGLE, stochastic EAGLE and SVD/QR methods.\n\n| **Task**            |   **Size**   | **Nonzeros** | **rank** | **conditioning** | **SVD/QR** | **CG** | **GD** | **EAGLE** | **Stoch. EAGLE** | **SGD** |\n|---------------------|:------------:|:------------:|:--------:|:--------------------:|:----------:|:------:|:------:|:---------:|:----------------:|:-------:|\n| Maragal_4/Maragal_4 |  1964 × 1034 |     1.3 %    |    801   | $1.38\\times10^{7}$ |      1     |   8.4  |   DNC  |    0.16   |        5.2       |   DNC   |\n| Maragal_5/Maragal_5 |  4654 × 3320 |     0.6 %    |   2147   | $1.22\\times10^{6}$ |      1     |   DNC  |   DNC  |    0.21   |        DNC       |   DNC   |\n| Meszaros/pf2177     | 9728 × 10178 |    0.031%    |   9662   |  $2.19\\times10^3$  |      1     |   1.4  |   DNC  |    1.3    |        DNC       |   DNC   |\n| HB/dwt_1007         |  1007 × 1007 |     0.85%    |   1000   |  $2.97\\times10^3$  |      1     |  0.21  |  0.97  |    0.64   |        6.3       |    10   |\n| HB/bcsstm13         |  2003 × 2003 |     0.53%    |   1241   |  $2.70\\times10^4$  |      1     |   DNC  |  0.024 |   0.058   |       0.048      |  0.024  |\n| Priebel/162bit      |  3606 × 3597 |     0.29%    |   3460   |  $1.68\\times10^4$  |      1     |   1.8  |   DNC  |    0.31   |        DNC       |   0.84  |\n| Schulthess/N_pid               |  3625 × 3923 |    0.057%    |   2048   |  $2.47\\times10^2$  |      1     |   2.4  |   DNC  |    0.26   |        DNC       |   DNC   |\n\n**Summary of New Results:**\n\n - EAGLE is **robustly accurate and efficient** on these SuiteSparse matrix tasks.\n    - Further, it converges even on problems where CG/GD fail to converge.\n - For **large ill-conditioned systems**, EAGLE's iterative approach is compelling, and gets both speed and accuracy, validating its practical potential.\n- Stochastic EAGLE improves upon SGD in large ill-conditioned problems. It is also more robust.\n\nWe will add the above experiments to the paper."}, {"review_id": "LP1x4G4uTb", "reviewer": "Reviewer_fRWb", "summary": "The paper train a linear transformer under 3 different regime (unconstrained, computed limited and distributed) on masked-block completion tasks showing that transformer is able to discover solution for the task under all regimes without any input on the task.", "questions": "- Which other tasks do you think can be interested to explore?\n- Why do you choose linear transformer? How the method can be extended to non linear regime?\n- Which real-world tasks do you think EAGLE can be implemented for?\n- How stable numerically is the proposed algorithm? Do you found places where it diverge?\n- Do you plan to open-source your code?", "limitations": "The paper focus only on linear transformer and only in one task - which limited the contribution on real world scenarios.\nI'll happy to hear which additional ideas the authors have to extend it to non linear regime and additional tasks.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "I think that the general idea of the paper is interesting and it is not a common area of research - I think it is a great bonus.\nThe paper regime is very limited and extension for additional models and regimes will be appreciated. It's unclear how transferable the EAGLE method is to real-world or non-linear tasks.", "was_revised": true, "final_justification": "I believe the entire regime is limited, so I decided to keep my score.", "rebuttal": "We thank the reviewer for their positive feedback on our paper's core idea and for recognizing its novelty. We appreciate the thoughtful questions regarding the choice of a linear transformer and the path to extending our work. We are happy to provide additional details and context.\n\n### On the Methodological Framework (Linear Transformer & Task Scope)\n\nThe reviewer's main concern is that our focus on a single task with a linear transformer limits the contribution. We would like to clarify that these were deliberate methodological choices aimed at making the discovery of the transformer's implicit algorithm tractable and interpretable.\n\n- **The simple architecture was a result of principled simplification, not a starting assumption.** Our goal was to isolate the computational essence of the learned algorithm. To do this, we performed a systematic simplification process, as described in the **Algorithm extraction** section (p. 4) and **Appendix B.5** (p. 22). For example, we found that performance was maintained even after pruning the model from eight attention heads down to a single head per layer. The resulting \"linear, attention-only\" model is therefore not an arbitrary limitation, but rather the outcome of a principled approach to uncover the fundamental algorithm learned by the transformer.\n\n- **The transformer learns a simple, structured hyperparameter scheme.** A crucial part of algorithm discovery is identifying the hyperparameters. A surprising result of our work is that the transformer simplifies this process immensely. [cite_start]After the simplification procedure, we found that the learned weights revealed a clear, interpretable rule for the step-size parameters of EAGLE. As shown in **Section 4** (\"Rationale for Weights\", p. 5) and **Figure 5** (p. 7), the learned weights correspond to fixed constants ($\\eta,\\gamma$) combined with an automatic normalization based on the data's spectral norm ($\\\\|A\\_\\ell\\\\|\\_{2}^{-2}$). This discovery transforms a potentially complex hyperparameter tuning problem into a simple, parameter-free update rule, providing strong evidence that the transformer is uncovering a robust and elegant numerical method.\n\n- **A \"linear transformer\" is not a linear model.** This is a critical point of clarification. A \"linear transformer\" is **not a linear function of its input data**. The attention update is fundamentally non-linear—specifically, it is **cubic in the data activations $Z$**. This is shown explicitly in the derivation in **Appendix C.3 (Eq. 4)**. The term 'linear' simply refers to the absence of an additional, point-wise non-linear activation (like a ReLU), which our simplification process found to be unnecessary for this task.\n\n### On Generalizing the \"EAGLE Method\" and Future Work\n\nWe agree with the reviewer that extending this work is an exciting direction. We believe reviewer comments may stem from a potential ambiguity in the term \"EAGLE method,\" which can refer to either the specific algorithm or our discovery process.\n\n- **The `EAGLE` algorithm is task-specific.** The specific two-line update rule we call EAGLE is the algorithm the transformer discovered for the Nyström-type matrix completion problem. It is a powerful second-order solver for this task, but we would not expect it to solve a different problem like sparse regression.\n\n- **The discovery *process* is highly general.** The core transferable contribution of our work is the **methodology for discovering the implicit algorithm**. We demonstrate this concretely in **Appendix D** (p. 25), where we apply this same process to a related but distinct **parameter estimation** task, revealing a slightly modified, highly efficient version of EAGLE. This shows the process is not a one-off trick. We are actively applying this framework to other domains like sparse regression and high-dimensional classification.\n\n### Answers to Specific Questions\n\n- **Why a linear transformer? How to extend to the non-linear regime?**\nAs explained above, the choice was a result of a principled simplification to enable algorithm discovery. The model is already non-linear in its data dependence. Our ongoing work on other tasks (e.g., classification) shows that this same discovery process reveals algorithms where non-linear activation functions (like ReLU or GeLU) are kept, as they are essential for those tasks.\n\n- **What real-world tasks can EAGLE be implemented for?**\nThe specific `EAGLE` algorithm is a second-order solver for least-squares problems based on Nyström approximation. It is therefore well-suited for real-world scenarios that rely on this technique, such as large-scale kernel methods, recommender systems, and statistical inference.\n\n- **How numerically stable is the algorithm? Do you find places where it diverges?**\nPractically, our method is stable. As we mention in the paper (lines 87, 188, 248 and appx C), the inner loop of EAGLE is closely related to the classical Newton-Schulz method for positive definite matrices. As such, we expect numerical stability similar to Newton-Schulz. Theoretically, this can get problematic when condition numbers start approaching 1/floating-point-precision, but empirically we expect stable behavior. This is demonstrated well in the Sparse-Suite Matrix Completion experiments (see response to reviewer rz5S), where we see very good performance.\n\n- **Do you plan to open-source your code?**\nYes, we plan to release our code upon publication to encourage further research in this area."}]

Author rebuttal:

Meta-review: (a) Summary: the paper shows that linear transformers can discover algorithms given input-output data, and the algorithms it discovers can adapt to compute and resource constraints. The study focuses on low-rank matrix completion problems, and the discovered algorithm, called EAGLE, is a simple and interpretable second-order algorithm. The paper establishes a quadratic convergence guarantee for EAGLE, and shows that in practice, it converges faster than linear convergence baselines in certain cases. 

(b) Strengths
- The topic of algorithm discovery by transformers is very relevant to the community, and the paper’s findings are original. 
- Different from the single-regime, multi-task approach that is common in prior works, the paper focuses on one task and incorporates several regimes. Specifically, it studies three regimes of algorithm learning: unconstrained learning, compute-constrained learning, and distributed learning, and finds that the transformer learns a unified algorithm that can adapt to all regimes. 
- The paper shows the significance and applicability of the discovered algorithm both in theory (establishes convergence guarantees) and in practice (benchmarking on several tasks). 

(c) Weaknesses: the main weaknesses of the work are the following:
- Limitation to linear transformers: even though making simplifying assumptions leads to a cleaner investigation, it’s unclear to what degree the findings transfer from linear transformers to the standard softmax attention architecture. It would be interesting to include results on non-linear transformers to increase the impact of the work.
- Limitation to one task: if showing algorithmic discovery in transformers is a central contribution of the paper, increasing the number of tasks will improve the generality of the result. 
- Scale of experiments: the paper could conduct larger scale experiments on the algorithm discovery side, or provide evidence whether the same algorithm is learned for higher-dimensional problems. This can provide further insight into the nature of algorithm discovery in transformers.

(d) Rationale: I recommend a borderline accept if there is room. The paper studies an interesting problem and the findings are novel, though there are a few limitations. The analysis on the discovered algorithm EAGLE is solid, and during the rebuttal phase, the authors provided additional experiments on higher-dimensional problems. However, the investigation uses linear transformers and is restricted to one task. Furthermore, the algorithm discovery aspect of the study can be more developed, for example by including more tasks and varying the size of the problem.

(e) Rebuttal discussions: after the rebuttal period, one reviewer maintained their very positive score, while the two remaining reviewers seem to have kept their initial borderline assessment. The main weaknesses raised were limitations to linear transformers, limitations to one task, and the performance of EAGLE. During the rebuttal, the performance of EAGLE was partially addressed after the authors provided additional experiments, but the reviewer would still like a more optimized implementation of EAGLE. The other reviewer fRWb finds the limited scope to warrant a borderline accept.
Decision: Accept (poster)