Paper deep dive

Optimizing Large Language Models with Chained LMOs

Sungyoon Kim, Kaan Ozkara, Youngsuk Park

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 163

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%

Last extracted: 10/9/2026, 8:08:48 AM

Summary

This paper introduces 'chained linear minimization oracles' (chained LMOs) as a unified framework for optimizing large language models, specifically addressing the fragmentation of Muon-based optimizers. The authors demonstrate that many existing Muon variants, which compose multiple normalizations, fall outside the standard LMO framework and can diverge on smooth convex objectives. However, they prove that chaining provides a provable learning advantage in linear associative memory settings with anisotropic embeddings. Empirically, they propose 'TensorChain', a novel optimizer that stacks weight matrices across layers and normalizes the resulting 3D tensor. TensorChain outperforms existing chained baselines in pretraining Qwen3 models, achieving 9.6% token savings over Muon.

Entities (8)

Muon · optimizer · 99%Chained LMOs · methodframework · 98%TensorChain · optimizer · 97%Qwen3 · model · 95%Matrix Sign (msign) · normalization-technique · 92%Row Normalization (RN) · normalization-technique · 92%Linear Associative Memory · theoretical-concept · 90%Aurora · optimizer · 85%

Relation Signals (7)

TensorChain → isevaluatedon → Qwen3

confidence 95% · In Qwen3 0.6B and 1.7B pretraining, TensorChain outperforms all chained baselines...

TensorChain → isinstanceof → Chained LMOs

confidence 95% · Empirically, we propose TensorChain, a novel optimizer within the framework...

TensorChain → outperforms → Muon

confidence 95% · TensorChain outperforms all chained baselines in average token efficiency, with average token savings of 9.6% over Muon at matched validation loss.

Chained LMOs → unifies → Muon variants

confidence 95% · We introduce chained linear minimization oracles (chained LMOs), which cast these methods as compositions of LMOs.

Chained LMOs → improvesover → Muon

confidence 92% · chaining can improve over Muon under anisotropic embeddings.

Chained LMOs → candivergeon → smooth convex objectives

confidence 90% · Despite their empirical success, many chains fall outside the standard LMO framework and can diverge on smooth convex objectives.

Row Normalization → iscomponentof → Chained LMOs

confidence 90% · Other norms recover other familiar normalizations: ... row normalization (RN) ... Recent Muon variants augment the matrix-sign step with exactly these normalizations...

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Muon has motivated a growing family of optimizers that compose multiple matrix normalizations, but these methods remain fragmented and lack a unified perspective. We introduce chained linear minimization oracles (chained LMOs), which cast these methods as compositions of LMOs. Despite their empirical success, many chains fall outside the standard LMO framework and can diverge on smooth convex objectives. To explain why composition can nevertheless help, we turn to linear associative memory and show that chaining can improve over Muon under anisotropic embeddings. Empirically, we propose TensorChain, a novel optimizer within the framework that stacks compatible weight matrices across different layers and normalizes the 3d tensor across its axes. In Qwen3 0.6B and 1.7B pretraining, TensorChain outperforms all chained baselines in average token efficiency, with average token savings of 9.6% over Muon at matched validation loss.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

161,009 characters extracted from source content.

Expand or collapse full text

Optimizing Large Language Models with Chained LMOs Sungyoon Kim Affiliation: Department of Electrical Engineering, Stanford University Email: sykim777@stanford.edu Kaan Ozkara Affiliation: Amazon Annapurna Labs Email: kaanozka@amazon.com Youngsuk Park Affiliation: Amazon Annapurna Labs Email: pyoungsu@amazon.com Abstract Muon has motivated a growing family of optimizers that compose multiple matrix normalizations, but these methods remain fragmented and lack a unified perspective. We introduce chained linear minimization oracles (chained LMOs), which cast these methods as compositions of LMOs. Despite their empirical success, many chains fall outside the standard LMO framework and can diverge on smooth convex objectives. To explain why composition can nevertheless help, we turn to linear associative memory and show that chaining can improve over Muon under anisotropic embeddings. Empirically, we propose TensorChain, a novel optimizer within the framework that stacks compatible weight matrices across different layers and normalizes the 3d tensor across its axes. In Qwen3 0.6B and 1.7B pretraining, TensorChain outperforms all chained baselines in average token efficiency, with average token savings of 9.6%9.6\% over Muon at matched validation loss. **footnotetext: Work done during an internship at Amazon Annapurna Labs. 1 Introduction After the success of the Muon optimizer (modded_nanogpt_2024) in pretraining large language models, numerous theoretical and empirical follow-up works have been proposed. Consider minimizing a differentiable objective f:ℝm×n→ℝf:R^m× n over a weight matrix W∈ℝm×nW ^m× n. Muon maintains a momentum matrix Mt∈ℝm×nM_t ^m× n, the exponential moving average (EMA) of the gradient, and updates the weights along its matrix sign: Mt=βMt−1+(1−β)∇f(Wt),Wt+1=Wt−ηtmsign(Mt).M_t=β M_t-1+(1-β)∇ f(W_t), W_t+1=W_t- _t (M_t). (1) Here, β∈[0,1)β∈[0,1) is the momentum coefficient and ηt>0 _t>0 is the learning rate. For a matrix M with singular value decomposition M=UΣV⊤M=U V , the matrix sign is defined as msign⁡(M)=UV⊤ (M)=UV ; in practice, it is approximated by a few Newton–Schulz iterations rather than by an explicit SVD. A principled way to understand this update is through the linear minimization oracle (LMO) (pethick2025training; kovalev2025understanding). For a norm ball =X:‖X‖≤ρD=\X:\|X\|≤ρ\, define lmo⁡(G)∈arg⁡minX∈⟨G,X⟩.lmo_D(G)∈ _X G,X . Over the spectral-norm ball, lmo⁡(G)=−ρmsign⁡(G)lmo_D(G)=-ρ\, (G), so Muon is exactly the iteration Mt=βMt−1+(1−β)∇f(Wt),Wt+1=Wt+ηtlmo(Mt),M_t=β M_t-1+(1-β)∇ f(W_t), W_t+1=W_t+ _tlmo_D(M_t), (2) with ρ=1ρ=1. This iteration updates in the steepest descent direction under the geometry of D, and inherits the standard convergence guarantees under the corresponding norm-smoothness assumptions (pethick2025training). Other norms recover other familiar normalizations: up to sign and scale, the max-row- and max-column-ℓ2 _2 balls yield row normalization (RNRN) and column normalization (CNCN), and the entrywise-ℓ∞ _∞ ball yields the elementwise sign (SgnSgn). Recent Muon variants augment the matrix-sign step with exactly these normalizations, applied before the matrix sign (chang2026muoneq; liu2026muon), after it (zhang2026muon+; li2025normuon), on both sides (si2025adamuon), or in repeated alternation (dewulf2026aurora). Although these methods were proposed separately, each with its own motivation and method-specific analysis, they share a common structure: each is a composition of LMOs. We introduce chained LMOs as a unifying framework for these variants and evaluate them in a unified experimental protocol. Given constraint sets 1,…,kD_1,…,D_k, we define the chained LMO T(G):=(−1)k+1(lmok∘⋯∘lmo1)(G),T(G):=(-1)^k+1 (lmo_D_k ·s _D_1 )(G), (3) and the corresponding iteration Mt=βMt−1+(1−β)∇f(Wt),Wt+1=Wt+ηtT(Mt).M_t=β M_t-1+(1-β)∇ f(W_t), W_t+1=W_t+ _tT(M_t). (4) The factor (−1)k+1(-1)^k+1 accounts for the negation introduced by each LMO. Muon is the single-LMO case k=1k=1 with the spectral-norm ball, while many recent variants become instances of eq. 4, either directly or, for methods with Adam-like state, after turning off the EMA terms following bernstein2024old (). Having identified chained LMOs as a common representation of Muon variants, a natural question is whether the favorable theory of individual LMOs survives composition. We show that for proposed Muon variants in the literature (zhang2026muon+; liu2026muon; chang2026muoneq; si2025adamuon; li2025normuon; dewulf2026aurora), it does not. This exposes the limitations of the existing LMO perspective and motivates a complementary view through linear associative memory, where we show that chaining can itself provide a learning advantage. Empirically, we compare these variants under a common protocol. Since each was originally evaluated in a different setup, their reported results are not directly comparable; a controlled comparison is therefore valuable in that it isolates the effect of the chain itself from differences in architecture, scale, and tuning. We further introduce a novel normalization axis that outperforms existing chained baselines under the same protocol, suggesting that the chained-LMO design space remains far from exhausted. Specifically, • A unified theory of chained LMOs. We introduce chained LMOs as a common representation of recent Muon variants and characterize the limits of the standard LMO perspective. Many practical chains are not LMOs of any fixed constraint set (Section 2.1). For positively aligned chains such as RN∘msignRN , a ν-approximate LMO framework recovers convergence guarantees, but with unavoidable dimension dependence (Proposition 4). Other chains admit negatively aligned directions and can diverge even on smooth convex objectives (Proposition 1). • A provable learning advantage from chaining. Building on recent analyses of Muon in linear associative memory (wang2026muon; kim2026sharp), we show that different normalizations can address distinct sources of anisotropy. Under anisotropic answer embeddings and power-law pair frequencies, msign∘RN can achieve asymptotically perfect retrieval while both Muon and gradient descent fail. • A new normalization axis and unified empirical comparison. We introduce TensorChain, which treats the layer dimension of stacked parameter matrices as a novel axis of normalization. Under a common training setup (Section 5.1), TensorChain outperforms existing chained Muon baselines in the 1M-token batch regime, and remains competitive in the 4M-token batch regime. The rest of the paper is organized as follows. In Section 2, we show that many chains in are not LMOs of any fixed constraint set, and analyze their convergence and divergence on smooth objectives. Section 3 offers a complementary perspective through linear associative memory, where chaining provably improves learning. Section 4 shows our exploration of different chains and introduces TensorChain, which is compared with existing chained-LMO methods under a matched training protocol Section 5. Section 6 concludes the paper. Related optimization frameworks and their relation to chained LMOs are discussed in Section A.1. 2 Optimizing with chained LMOs We work with the stateless chains in (Section A.2) and write N=Nk∘⋯∘N1N=N_k ·s N_1 for the corresponding chain of normalization maps; since each NiN_i is odd and invariant to positive rescaling, the chained LMO (equation 3) is T=−NT=-N up to scale. As −N⁡(G)∈argminX∈⁡⟨G,X⟩-N(G)∈ _X G,X exactly when N⁡(G)∈argmaxX∈−⁡⟨G,X⟩N(G)∈ _X∈-K G,X , chained LMO is an LMO over some set if and only if N is a linear maximization oracle over some set C, that is, N⁡(G)∈argmaxX∈⁡⟨G,X⟩N(G)∈ _X G,X for every G. For each chain, we therefore ask whether N is a linear maximization oracle over any set. 2.1 Chained LMOs are not necessarily LMOs Section A.2 shows that many Muon variants optimize with chained LMOs. We now make precise the claim from the introduction that these chains generally leave the LMO class, so that the standard LMO theory does not apply to them. The claim has to be established chain by chain, since some compositions do stay inside the class: for instance, positive row scaling preserves entrywise signs, so Sgn∘RN=SgnSgn =Sgn and the chain −(−Sgn)∘(−RN)=−Sgn-(-Sgn) (-RN)=-Sgn, which is indeed the LMO of the entrywise-ℓ∞ _∞ ball. We consider three groups: the chains used in practice (), finite alternations of RNRN and CNCN, and their limit ℰE. Theorem 1. Let m,n≥5m,n≥ 5, and let normalizations map zero rows and columns to zero. There is no fixed set ⊆ℝm×nK ^m× n that satisfies: for every G in its domain T⁡(G)∈argmaxX∈⁡⟨G,X⟩T(G)∈ _X G,X for each map T listed below: (i) msign∘RN , msign∘RN∘CN , RN∘msignRN , CN∘RN∘msignCN , and RN∘CN∘msignRN , for m>nm>n. (i) Alternation with the matrix sign: (msign∘RN)∘k( ) k and (RN∘msign)∘k(RN ) k, for every k≥1k≥ 1 and m>nm>n. (i) Every finite composition of RNRN and CNCN that contains both, in particular (RN∘CN)∘k(RN ) k; its limit ℰE; and msign∘ℰ , the last two on matrices with no zero entries. (iv) Every finite composition of msign and SgnSgn that contains both. Hence the map −T-T fails to be a LMO of any fixed set. The statement concerns the exact matrix sign and the stateless chains, and it covers every chain in up to transpose. Normalizing after msign can be vacuous: on tall inputs of full column rank, the columns of msign⁡(G) (G) are already orthonormal, so CN∘msign=msignCN = , and symmetrically, on wide inputs of full row rank, RN∘msign=msignRN = . The latter is why the restriction m>nm>n exists. The proof relies on properties that a linear maximization oracle over a fixed set C must satisfy. Lemma 1. Let N⁡(G)∈argmaxX∈⁡⟨G,X⟩N(G)∈ _X G,X for every G, for a fixed set C, and let Y=N⁡(G)Y=N(G) for some G. If ‖N⁡(Y)‖F=‖Y‖F\|N(Y)\|_F=\|Y\|_F, then N⁡(Y)=YN(Y)=Y. Lemma 2. Let N⁡(G)∈argmaxX∈⁡⟨G,X⟩N(G)∈ _X G,X for every G, for a fixed set C. If N⁡(−G)=−N⁡(G)N(-G)=-N(G), then ⟨G,N⁡(G)⟩≥0 G,N(G) ≥ 0. For (i)–(i) and the finite compositions in (i), the normalizations keep undoing each other: on a suitable input, applying the chain a second time changes the output but not its norm, contradicting Lemma 1. For msign∘ℰ and (iv), the first map is invariant to a family of transformations of its input: diagonal rescalings for ℰE, and changes of singular values or entry magnitudes for msign and SgnSgn. Transforming G within this family leaves the output unchanged but can make it negatively aligned with G; since these chains are odd, this contradicts Lemma 2. The limit ℰE is idempotent, odd, and positively aligned, so both lemmas hold for it; we instead violate the LMO condition directly, by finding G and H with ⟨G,ℰ⁡(G)⟩<⟨G,ℰ⁡(H)⟩ G,E(G) < G,E(H) . 2.2 Convergence behavior of chained LMOs An LMO update always has a negative alignment with the gradient, which makes it a descent direction. When there exists a positive alignment, which is possible for chained LMOs, we can construct a counterexample where an optimization method diverges even for smooth functions (Proposition 1). Proposition 1. Let T be positively scale invariant with ⟨G,T⁡(G)⟩>0 G,T(G) >0 for some G. For the convex, 11-smooth quadratic f⁡(X)=⟨G,X⟩2/(2‖G‖F2)f(X)= G,X ^2/(2\|G\|_F^2) with X0=G/‖G‖F2X_0=G/\|G\|_F^2, the iterates Xt+1=Xt+ηtT(∇f(Xt))X_t+1=X_t+ _tT(∇ f(X_t)) satisfy f⁡(Xt)→∞f(X_t)→∞ for every step size sequence with ∑tηt=∞ _t _t=∞. The counterexamples of Theorem 1 give such inputs for every composition of msign and SgnSgn that contains both, including the stateless chains of Muon2 and AdaMuon, and for msign∘ℰ \ E in all sufficiently large shapes. The length-4 chain of Aurora, (msign∘RN)∘2( ) 2, also admits such input for n=2n=2 and m≥258m≥ 258 (Lemma B.11). None of these methods can therefore have a convergence guarantee even over smooth functions. Approximate LMOs. For chains with negative alignment, we measure how far they are from an exact LMO. We formalize this notion using approximate LMOs. Definition 1 (Approximate LMO). Let K be a compact set symmetric with respect to the origin, and define m(G)=minX∈⁡⟨G,X⟩m_K(G)= _X G,X . For ν∈(0,1]ν∈(0,1], a map U is a ν-approximate LMO of K if U⁡(G)∈U(G) and ⟨G,U⁡(G)⟩≤νm(G) G,U(G) ≤ν\,m_K(G) for every G. The feasibility condition U⁡(G)∈U(G) prevents the approximation factor from being increased trivially by rescaling U. It also implies ν≤1ν≤ 1, where ν=1ν=1 implies exact LMO. Using the same machinery from the LMO framework, we can prove convergence of approximate-LMOs. Proposition 2. Let f have an L-Lipschitz gradient, let Δ=f⁡(X0)−inff =f(X_0)- f, and let U be a ν-approximate LMO of K. Define R:=maxX∈⁡‖X‖FR:= _X \|X\|_F. (i) Gradient descent. With Xt+1=Xt+ηU(∇f(Xt))X_t+1=X_t+η\,U(∇ f(X_t)) and η=2Δ/(LR2T)η= 2 /(LR^2T), 1T∑t=0T−1h(∇f(Xt))≤Rν2LΔT. 1T _t=0^T-1h_K(∇ f(X_t))≤ Rν 2L T. (i) SGD with momentum. With stochastic gradients of variance at most σ2σ^2, Xt+1=Xt+ηU(Mt)X_t+1=X_t+η\,U(M_t) with the momentum of Equation 4, 1−β=T−1/21-β=T^-1/2 and η=cT−3/4η=c\,T^-3/4, the same average in expectation is O(ν−1[Δ/c+R(LRc+σ)]T−1/4)O (ν^-1[ /c+R(LRc+σ)]T^-1/4 ). Compared with an exact LMO, both rates lose exactly a factor 1/ν1/ν. Two chains from , msign∘RN and RN∘msignRN , fall into the class of approximate-LMOs (Proposition 3). Proposition 3. On all m×nm× n matrices, −RN∘msign-RN is a 1/m1/ m-approximate LMO of row=X:maxi∥Xi:∥2≤1K_row=\X: _i\|X_i:\|_2≤ 1\, and −msign∘RN- is a 1/m1/ m-approximate LMO of the spectral-norm ball. By transposition, the same holds for CNCN with 1/n1/ n. This gives convergence guarantees for the stateless versions of NorMuon, Muon++ with a single normalization and MuonEq, with additional m m appearing in the upper bound. The additional dimension factor cannot be removed within this analysis: any set for which a map is feasible must contain the convex hull of its image, and a larger set only makes the approximation factor smaller. Proposition 4. Let n=2n=2. For every closed convex set K containing the image, the best approximation factor satisfies ν≤3(m−1)−1/4ν≤ 3(m-1)^-1/4 for −RN∘msign-RN (for odd m) and ν≤3/m−1ν≤ 3/ m-1 for −msign∘RN- . Both bounds tend to zero, so any guarantee derived from Proposition 2 degrades with dimension. A tighter analysis may need to go beyond uniform approximation, e.g. by bounding the approximation factor along the trajectory, which is out of scope for this work. 3 Provable Benefit of Chaining in Linear Associative Memory Learning Section 2 shows that the standard worst-case guarantees of LMO-based optimization do not generally survive chaining: some chains can diverge even on smooth convex objectives, while others admit only dimension-dependent approximate-LMO guarantees. Any positive theory for chaining must therefore exploit additional problem structure beyond generic smoothness. We instead demonstrate a benefit on the linear associative memory (Section A.3), following recent theoretical work on Muon (wang2026muon; kim2026sharp; li2026muon). Muon helps there by whitening the anisotropy of the label frequencies, which suffices when the frequencies are the only source of anisotropy, as for the isotropic Gaussian embeddings of kim2026sharp. With a second source, such as anisotropic embeddings, Muon alone may not suffice, and chains such as msign∘RN can provably help. Anisotropic embeddings also make the model a closer abstraction of trained transformers, whose hidden representations are dominated by a few outlier coordinates (timkey2021all; dettmers2022gpt3). We make this argument precise in the following sections. We store N facts with queries yi∈ℝny_i ^n and answers vi=Σxi∈ℝdv_i= x_i ^d, where yi∼⁡(0,In/n)y_i (0,I_n/n) and xi∼⁡(0,Id/d)x_i (0,I_d/d) are independent, and Σ=diag⁡(σ)∈ℝd×d = (σ) ^d× d has positive entries. Fact i occurs with probability pi=i−α/∑r=1Nr−αp_i=i^-α/ _r=1^Nr^-α for a fixed α>1α>1. A memory W∈ℝd×nW ^d× n is trained with the cross-entropy loss and evaluated by its retrieval accuracy, ℒ(W)=∑i=1Npi(−vi⊤Wyi+log∑j=1Nevj⊤Wyi),Acc(W)=1N∑i=1Nvi⊤Wyi>maxj≠ivj⊤Wyi.L(W)= _i=1^Np_i (-v_i Wy_i+ _j=1^Ne^v_j Wy_i ),\ \ Acc(W)= 1N _i=1^N 1 \v_i Wy_i> _j≠ iv_j Wy_i \. Following prior theoretical analysis of linear associative memory (kim2026sharp), we study the retrieval accuracy after a single update from W0=0W_0=0. Let G=−∇ℒ(0)G=- (0) denote the negative gradient at initialization. We compare three methods, gradient descent, Muon, and row normalization followed by Muon: WGD=G,WMuon=msign⁡(G),Wchain=msign⁡(RN⁡(G)).W_GD=G, W_Muon= (G), W_chain= (RN(G)). Step sizes are irrelevant because retrieval accuracy is invariant to positive rescaling of W. This chain is tailored to the anisotropy in our setting, and other forms of anisotropy may give provable advantages to other chains. Asymptotic regime. Fix α, an infinite positive sequence σ=(σj)j≥1σ=( _j)_j≥ 1, and dimension sequences Nn,dnN_n,d_n. For each n, the instance uses N=NnN=N_n facts, query dimension n, answer dimension d=dnd=d_n, and the scales σ1:dn _1:d_n. For a method M∈GD,Muon,chainM∈\GD,Muon,chain\, write AccnM=Acc(WM;Nn,n,dn,σ1:dn)Acc_n^M=Acc (W_M;\,N_n,n,d_n, _1:d_n ) for the random accuracy of the instance. The quantity of interest is its limit in probability as n→∞n→∞. Assumption 1 (Scaling). Fix 0<ε<(β−1/2)/20< <(β-1/2)/2 and γ>1/(1−β)γ>1/(1-β). Eventually, n1+(1−β+ε)/α≤Nn≤n1+(1/2−ε)/αn^1+(1-β+ )/α≤ N_n≤ n^1+(1/2- )/α and n1/(1−β)≤dn≤nγn^1/(1-β)≤ d_n≤ n^γ. Assumption 2 (Anisotropy). For some fixed 1/2<β<11/2<β<1 and constants 0<cσ≤Cσ<∞0<c_σ≤ C_σ<∞, we have cσj−β≤σj≤Cσj−βc_σj^-β≤ _j≤ C_σj^-β for every j≥1j≥ 1. The scaling regime (Assumption 1) is roughly n≪N≪dn N d: there are more facts than query dimensions, so retrieval requires extracting a weak fact-specific signal, while the large answer dimension provides many coordinates over which that signal can accumulate. Assumption 2 ensures that ∑jσj2 _j _j^2 remains bounded while ∑j≤dσj _j≤ d _j diverges. Thus row normalization, which removes the answer-coordinate scales before the matrix sign, allows the chain to accumulate signal across coordinates in a regime where GD and Muon fail. The core result is that the asymptotic one-step accuracy of chained Muon tends to 1, whereas those of Muon and GD go to 0. Theorem 2 (One-step separation). Under Assumptions 1 and 2, as n→∞n→∞, (Accnchain,AccnMuon,AccnGD)→ℙ(1,0,0). (Acc_n^chain,Acc_n^Muon,Acc_n^GD ) P(1,0,0). The separation is driven by how the answer-coordinate scales enter the three updates. For GD and Muon, the relevant normalized score quantities depend on quadratic and quartic powers of σa _a; because ∑aσa2<∞ _a _a^2<∞ and ∑aσa4<∞ _a _a^4<∞, their normalized score vectors have only O⁡(n)O(n) total energy, while the correct-answer signal remains too weak to overcome the maximum of the N≫nN n incorrect scores. By contrast, row normalization removes the answer-coordinate scales before the matrix-sign step, since RN⁡(ΣXK)=RN⁡(XK)RN( XK)=RN(XK). The useful contributions can then add coherently in proportion to ∑aσa _a _a, while the residual fluctuations are controlled by ∑aσa2 _a _a^2. Since ∑a≤dσa→∞ _a≤ d _a→∞ while the quadratic and quartic sums remain bounded, the chained update can accumulate signal across answer coordinates in a regime where GD and Muon cannot, yielding the separation above. The full proof is deferred to Appendix C. 4 Exploring the Space of Chains Section 3 suggests that equalizing the update across coordinates is what makes chaining useful. We push this further by exploring three directions in the space of chained LMOs. First, we introduce TensorChain, which adds a new normalization axis across layers. Second, we explore Aurora-LN, a simple long-side generalization of Aurora that applies the same chain to both tall and wide matrices. Third, motivated by recent work replacing Newton–Schulz with simpler axis-wise normalization (glentis2025minimalist; deng2026rmnp), we study whether the spectral step itself can be replaced by iterated row–column normalization, which we refer to as SinkhornChain. Algorithm 1 TensorChain 1: Momenta M1:g∈(ℝm×n)gM_1:g∈(R^m× n)^g 2: T←Stack⁡(M1,…,Mg)∈ℝg×m×nT (M_1,…,M_g) ^g× m× n 3: A←ANlayer∘ANcol∘ANrowA _layer _col _row 4: T←A∘k(T)T← A k(T) 5: T←NSslice(T)T _slice(T) 6: T←A∘k(T)T← A k(T) 7: return Unstack⁡(T)Unstack(T) Algorithm 2 Aurora-LN 1: momentum M∈ℝm×nM ^m× n 2: if m≥nm≥ n then 3: U←(NS∘RN∘NS∘RN)(M)U←(NS )(M) 4: else 5: U←[(NS∘RN∘NS∘RN)(M⊤)]⊤U← [(NS )(M ) ] 6: end if 7: return U 4.1 TensorChain: Normalizing Across Layers Existing chains normalize only within individual weight matrices, through row, column, or spectral normalization, and treat each layer independently. TensorChain instead jointly normalizes parameters of the same shape across layers (Algorithm 1). Parameters with the same shape and role, for example, the query projections WQW_Q of transformer blocks, are partitioned into groups. For a group of size g, we stack their momenta into a tensor M∈ℝg×m×n,M ^g× m× n, whose first axis indexes layers. Normalizing fibers along this axis prevents individual layers from dominating the corresponding coordinates of the group’s update. Let ANlayer,ANrow,ANcolAN_layer,AN_row,AN_col denote ℓ2 _2-normalization along the corresponding tensor axes, define :=ANlayer∘ANcol∘ANrow,A:=AN_layer _col _row, and let NSsliceNS_slice apply Newton–Schulz independently to each matrix slice. TensorChain maps TensorChain⁡(M)=∘k∘NSslice∘k(M),TensorChain(M)=A k _slice k(M), whose slices give the per-layer updates. TensorChain fits directly into the chained-LMO framework of Section 2 when we assume exact msign instead of Newton–Schulz. We will denote this as idealized TensorChain. Each fiber normalization is an exact LMO over a product of ℓ2 _2 balls, while slice-wise matrix sign is an exact LMO over the corresponding product of spectral-norm balls. Proposition 5 (Idealized TensorChain as chained LMO). For k≥1k≥ 1, Idealized TensorChain is a composition of exact LMOs over product norm balls, but the resulting map is not itself an LMO of any fixed constraint set when m,n≥5m,n≥ 5. Moreover, Idealized TensorChain map admits a negative-alignment input and can diverge on a convex smooth objective under non-summable step sizes. A closely related approach is Teon (zhang2026teon), which also stacks matrices across layers but applies Newton–Schulz to a 2d matricization of the tensor. TensorChain instead applies Newton–Schulz independently to each matrix slice, with axis-wise normalization before and after. Empirical comparison with the Teon optimizer can be found in Section F.7. We can consider the overhead of TensorChain in two axes: additional flops due to additional normalization, and additional communication cost to construct the tensor. Fiber normalizations cost O⁡(gmn)O(gmn) per group, negligible compared with the O⁡(gmnmin⁡(m,n))O(gmn (m,n)) Newton–Schulz steps, so TensorChain matches Muon’s optimizer FLOPs up to lower-order terms. Under FSDP, layer-axis normalization couples matrices within a group and requires additional communication. Empirical performance of TensorChain can be found in section 5.2. 4.2 Aurora-LN: Normalizing Along the Long Side Aurora (dewulf2026aurora) includes several design choices beyond the normalization chain itself: the Aurora update is applied only to tall matrices, while wide matrices fall back to Muon, and a β parameter accumulates the diagonal row preconditioner across steps. From the perspective of chaining, we consider a simplified variant that removes this accumulation and applies the same chain along the long side of every matrix. For tall and square matrices we apply the chain directly, while for wide matrices we apply it to the transpose and transpose the resulting update back. We refer to this variant as Aurora-LN (Algorithm 2, LN stands for Longside Normalization). Long-side or wide-matrix variants of Aurora-style normalization have previously appeared in modded-nanoGPT experiments (modded_nanogpt_2024), so we do not present Aurora-LN as a novel optimizer or algorithmic primitive. Rather, we revisit this design choice as a controlled simplification of Aurora. Our experiments in Section 5.2 show that this simplified variant can be more effective than the default Aurora implementation. 4.3 Replacing Newton–Schulz with Sinkhorn Normalization A complementary direction is to ask whether chaining can eliminate the spectral step altogether. Motivated by recent work on inexpensive normalizations replacing Newton–Schulz (deng2026rmnp; glentis2025minimalist), we consider SinkhornChain, which alternates row and column rescaling as in Sinkhorn balancing, but uses ℓ2 _2 rather than ℓ1 _1 normalization. For an m×nm× n matrix with m>nm>n, each row–column normalization sweep costs O⁡(mn)O(mn) FLOPs, whereas each Newton–Schulz step costs O⁡(mn2)O(mn^2). Thus, for a fixed number of iterations, SinkhornChain has substantially lower arithmetic complexity than Newton–Schulz. However, none of our tested normalization variants or layerwise replacement configurations outperformed Muon; the strongest replaces NS only on Q,KQ,K, reaching validation loss 2.59452.5945 versus Muon’s 2.59052.5905. Nevertheless, we still report these results, as in large-scale training the computational cost of Muon may be problematic and Sinkhorn shows a possibility of replacing it. Full details and ablations are given in Appendix E. 5 Experiments 5.1 Setup Models and data. We train Qwen3 models (yang2025qwen3) of sizes 0.6B and 1.7B on FineWeb (penedo2024fineweb) with sequence length 4096, using pre-norm RMSNorm (zhang2019root), QK-norm, untied embeddings, and the default torchtitan initialization (liang2025torchtitan). Training uses a 30B-token corpus, with experiments exceeding this budget spanning multiple epochs. To assess the effect of data repetition, we additionally conduct experiments on a 112B-token corpus without repeated training tokens, obtaining similar results and consistent relative performance across best performing optimizers (see Figure 10). We report loss on a held-out FineWeb split of 26M tokens. Training. We evaluate token batch sizes of 1M and 4M. The 0.6B model uses 8-way FSDP, while the 1.7B model uses 4-way data parallelism and 8-way FSDP. Learning rate schedules of both cosine decay with 10% warmup and WSD results (10% warmup and 10% decay) were evaluated. Matrix parameters use the optimizer under comparison with Nesterov momentum, with each update U rescaled by 0.2mn/‖U‖F0.2 mn/\|U\|_F (liu2025muon). Embeddings, output, and normalization parameters use AdamW with learning rate 8×10−38× 10^-3, and all parameters use decoupled weight decay 0.10.1 (loshchilov2017decoupled). We use a fixed random seed. Baselines and tuning. We compare the methods in together with the two explorations in Section 4, TensorChain and Aurora-LN. We use the official implementation where available and our own otherwise, with the same five-step bfloat16 Newton–Schulz for all methods. For TensorChain, we grouped all the weights that have the same role as a same group for simplicity, and used k=5k=5. On 0.6B at 1×1× Chinchilla (20 tokens per non-embedding parameter; hoffmann2022training), we sweep learning rates over 1,2,3,5,8×10−3\1,2,3,5,8\× 10^-3 for WSD and 5,8,12,16,20×10−3\5,8,12,16,20\× 10^-3 for cosine decay. Across these sweeps, we observe that most baselines attain their best validation loss at the same learning rate within each schedule: 3×10−33× 10^-3 under WSD and 12×10−312× 10^-3 under cosine decay. This motivates a matched-learning-rate comparison in which the selected baselines are evaluated at a shared learning rate. AdaMuon is the only clear exception. Under WSD, it shares the same optimum and shows no clear advantage. However, under cosine decay, its optimum shifts to a substantially larger learning rate and its performance improves markedly, though it remains worse than TensorChain. We therefore exclude AdaMuon from the matched-learning-rate comparison, report its full sweep separately in Appendix F, and leave a systematic study to future work. Among the remaining baselines, we select MuonEq, Aurora-LN, and Muon2 as the best-performing representatives using one Newton–Schulz, two Newton–Schulz, and additional momentum-dependent state, respectively, while retaining Muon and Aurora as references. For 4×4× Chinchilla, we rescale it by 4−0.32≈0.644^-0.32≈ 0.64 from the learning rate tuned at 1x Chinchilla (bjorck2025scaling). 5.2 Results Figure 1: Top: Comparison of Muon, MuonEq and TensorChain; Bottom: Comparison of Aurora and Aurora-LN, both setups are training on Qwen3-0.6B with a 1M-token batch at 1×1× Chinchilla under cosine LR. For other setups see section F.6. Table 1: Final validation loss under WSD and cosine schedules. Bold marks the best method in each row. The two methods on the right (Aurora, Aurora-LN) use Newton–Schulz twice, whereas other methods apply Newton–Schulz once. Muon2 has an additional momentum state. (*) were introduced in Section 4. Schedule Model Batch Tokens Muon MuonEq TensorChain∗ Muon2 Aurora Aurora-LN∗ WSD 0.6B 1M 8.8B 2.9201 2.9151 2.9000 2.9104 2.9184 2.9116 WSD 1.7B 4M 112B 2.5310 2.5279 2.5236 2.5264 2.5301 2.5242 Cosine 0.6B 1M 8.8B 2.8639 2.8540 2.8379 2.8497 2.8624 2.8474 Cosine 0.6B 4M 35B 2.7527 2.7434 2.7344 2.7443 2.7487 2.7391 Cosine 1.7B 1M 28.1B 2.5974 2.5924 2.5857 2.5885 2.5964 2.5882 Cosine 1.7B 4M 112B 2.4823 2.4784 2.4747 2.4744 2.4810 2.4734 Main comparison. TensorChain shows a consistent advantage over the more directly comparable chained baselines. Under WSD, it achieves the best final validation loss in both reported settings; under cosine decay, it is best in three of the four settings and outperforms Muon and MuonEq throughout. The only cases in which TensorChain is surpassed are by methods that use either an additional Newton–Schulz application or additional momentum-dependent state. Thus, among variants with comparable computational and state complexity, TensorChain provides a consistent improvement across model scales, batch sizes, and learning-rate schedules. The original Aurora performs only slightly better than Muon in our experiments, whereas the simplified Aurora-LN variant consistently improves over Aurora and becomes the best-performing method in the largest cosine setting. This suggests that seemingly minor choices in how a chain is applied across matrix shapes can have a substantial effect on performance. Figure 1 visualizes these two comparisons: TensorChain against Muon and MuonEq, and Aurora-LN against the original Aurora. Table 2: Token-efficiency gain of TensorChain: percentage(%)(\%) reduction in tokens needed to reach a common validation-loss target, estimated by linear interpolation between adjacent checkpoints. Positive values favor TensorChain. Model Batch Tokens Aurora-LN Muon2 MuonEq Aurora Muon 0.6B 1M 8.8B 8.2 9.0 10.3 12.4 12.8 0.6B 4M 35B 5.4 7.6 7.4 8.9 10.0 1.7B 1M 28.1B 3.9 4.0 6.2 7.5 7.8 1.7B 4M 112B -1.9 -0.5 6.7 7.4 7.8 The token-efficiency comparison in Table 2 makes the same trend more explicit. We report the percentage reduction in training tokens required to reach the same validation loss, estimated by linear interpolation between adjacent checkpoints; within each fixed-batch setting, this is equivalently the reduction in training steps. Across the four settings, TensorChain requires on average 9.6%9.6\% fewer tokens than Muon and 7.7%7.7\% fewer than MuonEq. 5.3 TensorChain Ablations Grouping. We compare consecutive, stride, greedy load-balanced, and all- compatible-layer grouping (Table 3); implementation details are deferred to Section F.5. Consecutive grouping is consistently worse with validation loss gap of ≈0.007≈ 0.007 for WSD and ≈0.01≈ 0.01 for cosine, while stride, greedy, and all performed similarly compared to consecutive grouping. This suggests that mixing matrices across distant layers matters. Table 3: TensorChain grouping ablation on Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla. Each entry reports the best validation loss from the corresponding full learning-rate sweep. Schedule Consecutive Stride Greedy All WSD 2.9070 2.8990 2.8994 2.9000 Cosine 2.8474 2.8357 2.8370 2.8379 Cross-layer normalization and iteration count. We conduct ablation studies to investigate two aspects of TensorChain: the contribution of cross-layer normalization and the effect of repeated normalization. All experiments use Qwen3-0.6B with a 1M-token batch size and a (1×) Chinchilla training budget, under both WSD and cosine learning-rate schedules. To isolate the effect of cross-layer normalization, we compare TensorChain against Sink–NS–Sink, which applies iterative row and column normalization before and after Newton–Schulz orthogonalization without stacking parameters across layers. As shown in Table 4, TensorChain consistently achieves lower validation loss across all settings, with improvements of 0.0208–0.0220 under WSD and 0.0248–0.0250 under cosine. These results highlight the contribution of cross-layer normalization beyond within-matrix normalization alone. We next examine the effect of the number of normalization iterations k, varying it from 1 to 5 for both methods. Increasing k does not consistently improve performance: validation loss remains nearly unchanged across iteration counts, with (k=1k=1) achieving the lowest or tied-lowest loss in every setting. This suggests that the benefits of TensorChain arise primarily from cross-layer normalization rather than convergence of the normalization procedure, making smaller k an attractive choice when computational efficiency is important. Table 4: Ablation of cross-layer normalization and iteration count k on Qwen3-0.6B (1M-token batch, 1×1× Chinchilla). Each entry reports final validation loss. Schedule Method k=1k=1 k=2k=2 k=3k=3 k=4k=4 k=5k=5 WSD TensorChain 2.8992 2.9003 2.9005 2.9002 2.9000 Sink–NS–Sink 2.9210 2.9218 2.9214 2.9216 2.9213 Cosine TensorChain 2.8368 2.8381 2.8386 2.8382 2.8379 Sink–NS–Sink 2.8624 2.8630 2.8626 2.8624 2.8627 6 Conclusion We introduced the chained-LMO framework as a unified view of recent Muon variants. Although chaining generally falls outside the standard LMO framework and therefore may lose its usual smooth optimization guarantee, it can still provide substantial practical benefits. To close this gap, we analyzed linear associative memory with multiple sources of anisotropy and proved that chaining can yield a strict learning advantage. Empirically, we introduced TensorChain, which adds normalization across the layer axis, and showed that it can outperform existing chained variants in language-model pretraining. An interesting direction for future work is to better understand how TensorChain’s gains vary across model scales and training settings, including the particularly strong improvements observed at 0.6B. On the theory side, extending the associative-memory analysis to more general chains and forms of anisotropy would clarify when chaining should help. Finally, TensorChain naturally extends to additional tensor axes arising in multi-head attention and mixture-of-experts models, where understanding the potential gains is an interesting direction for future work. AI use statement In this paper, AI tools were used in the following capacities: • Assistance in proof verification and organization • Assistance in novel proofs in Section 3; specifically on probabilistic rigor • Coding assistants were used to help implement the proposed methods, run experimental pipelines, and retrieve and organize experimental results. References Appendix A Appendix: Related work A.1 Frameworks around Muon and variants bernstein2024old take a descent step that minimizes the local quadratic model argminΔW⟨∇f(W),ΔW⟩+λ2‖ΔW‖2argmin_ W ∇ f(W), W + λ2\| W\|^2 with the norm ∥⋅∥\|·\| fixed in advance. The solution splits into a step size and a direction, −(∥∇f(W)∥†/λ)⋅argmax‖X‖=1⟨∇f(W),X⟩-(\|∇ f(W)\| /λ)· _\|X\|=1 ∇ f(W),X . Taking ∥⋅∥\|·\| to be the induced operator norm ‖W‖α→β=maxx⁡‖Wx‖β/‖x‖α\|W\|_α→β= _x\|Wx\|_β/\|x\|_α gives a family of optimizers including popular adaptive algorithms. Switching off their EMAs, Adam (kingma2014adam), Shampoo (gupta2018shampoo) and Prodigy (mishchenko2023prodigy) are each steepest descent under a particular member: (α,β)=(ℓ1,ℓ∞)(α,β)=( _1, _∞) for Adam and Prodigy, which differ only in their step-size rule, and (ℓ2,ℓ2)( _2, _2) for Shampoo, which is a positive scaling of the Muon update when the EMA is switched off. pethick2025training write the same direction as a linear minimization oracle lmoℬ(∇f(W))=argminX∈ℬ⟨∇f(W),X⟩lmo_B(∇ f(W))=argmin_X ∇ f(W),X over a norm ball, and interpret existing normalized gradient methods as optimizing with LMOs. By understanding the induced norm as how the output is amplified in a linear layer under bounded input, pethick2025training propose the Scion optimizer, whose design principle is to bound the RMS norm of intermediate layers in neural network training. kovalev2025understanding interprets the constraint set ℬB as a trust region, and both works prove convergence guarantees for their respective classes of algorithms. These frameworks cast optimizer choice as an appropriate choice of norm or constraint set, and subsequent work varies that choice along several axes. lau2025polargrad write vector preconditioning and matrix-aware updates as the same subproblem under different norms: the time-dependent quadratic metric gives the variable-metric methods that date back to fletcher1970new and are made data-driven in AdaGrad (duchi2011adaptive), while the spectral norm gives ‖G‖nucmsign⁡(G)\|G\|_nuc\, (G), the solution of bernstein2024old. crawshaw2025exploration aggregate per-layer norms into a product norm, whose oracle decomposes into scaled per-layer oracles; this expresses the standard practice of running Muon and Adam on separate parameter classes. Others vary the geometry itself: the Schatten-p family interpolates between the Frobenius and spectral cases (ma2026schattor; pethick2026use), and preconditioned matrix norms widen the family further (veprikov2025preconditioned). Certain frameworks deviate from the norm-choice template. Lion-K (chen2024lion) replaces the sign map by the subgradient of a convex K, and shows that the iterates converge to KKT points of minX⁡F⁡(X) _XF(X) subject to ∗(X)≤1/λK^*(X)≤ 1/λ for decoupled weight decay λ. Muon is the special case =∥⋅∥trK=\|·\|_tr (chen2025muon). shumaylov2026muon instead exhibit updates of the form (GG⊤)−cG(G )^-cG whose best-performing exponents fall in the quasi-norm regime c>1/2c>1/2, where the update is provably not the LMO of any unitarily invariant norm ball; they establish convergence directly, through alignment and a descent potential rather than through the oracle. These frameworks share one structural feature: a single geometry determines the update map, and the same geometry supplies the analysis. An LMO over a norm ball is aligned with the negative gradient, ⟨G,lmoℬ(G)⟩=−ρ‖G‖∗ G,lmo_B(G) =-ρ\|G\|_*, which yields a descent lemma under the matching norm-smoothness assumption; Lion-K replaces this with a Lyapunov function built from K. Chained LMOs do not in general inherit either structure, because the LMO class is not closed under composition. Even when every constituent map is an LMO, many practical chains are not the LMO of any fixed constraint set (Section 2), nor in general the subgradient of a convex K, so neither the oracle-based nor the Lyapunov-based argument transfers. We therefore analyze chains directly, in the spirit of shumaylov2026muon. For positively aligned chains, an α-approximate LMO framework recovers convergence, at the cost of an unavoidable dimension dependence. Chains that admit negative-alignment directions, in contrast, can diverge even on smooth convex objectives (Section 2.2). Closest to our setting is lau2026symmetry, whose contribution is a design principle: update maps derived from the symmetry group of each weight block, with NorMuon, Muon++, MuonEq and Aurora (li2025normuon; zhang2026muon+; chang2026muoneq; dewulf2026aurora) placed within that perspective. We instead study the effect of chaining itself: how composed updates relate to existing frameworks, whether they converge, and how they can improve linear associative memory learning. Their construction composes row-norm and one-sided spectral maps only, so the elementwise chains of si2025adamuon and liu2026muon fall outside it; our framework places no such restrictions. A.2 Existing Muon variants as chained LMOs Recent Muon variants add normalization at different stages of the orthogonalization pipeline. Post-orthogonalization methods correct the row- or column-wise imbalance left after the matrix-sign step. NorMuon (li2025normuon) does so with a row-wise second-moment state, while Muon++ (zhang2026muon+) applies stateless row, column, column–row, or row–column normalization directly after msign . Pre-orthogonalization methods instead precondition the momentum before Newton–Schulz to improve the spectrum entering the orthogonalization step: MuonEq (chang2026muoneq) uses row/column equilibration, and Muon2 (liu2026muon) uses Adam-style second-moment preconditioning. AdaMuon (si2025adamuon) acts on both sides, applying an elementwise sign before orthogonalization and a second-moment normalization afterward, which combines coordinate-wise adaptivity with orthogonal updates. Aurora (dewulf2026aurora) alternates row normalization and the matrix sign on tall matrices through the chain (msign∘RN)∘2( ) 2, so that the update is approximately balanced in both its spectrum and its row norms. Some of these optimizers are stateful because of exponential moving averages. bernstein2024old recover norm-based steepest-descent methods from Adam (kingma2014adam), Shampoo (gupta2018shampoo), and Prodigy (mishchenko2023prodigy) by switching off their EMAs, and we use the same reduction to expose the underlying normalization maps. For example, Adam with β1=β2=0 _1= _2=0 reduces to SgnSgn; similarly, NorMuon reduces to RN∘msignRN , Muon2 to msign∘Sgn , and AdaMuon to Sgn∘msign∘SgnSgn . We use this reduction only to isolate the normalization chain, not to identify the full stateful optimizer with its stateless counterpart. A.3 Linear Associative Memory and Muon An associative memory stores key–value pairs (ki,vi)i=1N\(k_i,v_i)\_i=1^N, with ki∈ℝdkk_i ^d_k and vi∈ℝdvv_i ^d_v, in a matrix W∈ℝdv×dkW ^d_v× d_k so that presenting kik_i recalls viv_i; for orthonormal keys, W=∑iviki⊤W= _iv_ik_i does so exactly. Attention can be read this way (zhong2025understanding): without the softmax normalization it computes ot=∑i≤texp⁡(qt⊤ki/dk)vi=Stϕ(qt)o_t= _i≤ t (q_t k_i/ d_k)\,v_i=S_tφ(q_t) with St=∑i≤tviϕ(ki)⊤S_t= _i≤ tv_iφ(k_i) , where ϕφ is the infinite-dimensional feature map of the exponential kernel, and FFN layers are similar through ReLU kernels. From an interpretability perspective, factual associations in a transformer can be edited by a rank-one update to an MLP projection (meng2022locating; meng2022mass). A growing line of work relates the success of Muon to associative memory. wang2026muon find empirically that Muon’s advantage over Adam comes mainly from the value/output and FFN weights, precisely the parameters that act as associative memories. To study this theoretically, wang2026muon; kim2026sharp; li2026muon turn to the linear associative memory, where W retrieves pair i correctly when vi⊤Wki>vj⊤Wkiv_i Wk_i>v_j Wk_i for all j≠ij≠ i, and is trained with the cross-entropy loss ℒL on these logits, with pair i occurring with frequency pip_i. When the frequencies are imbalanced, Muon provably retrieves better than gradient descent, because it whitens the frequencies. At initialization, −∇ℒ(0)≈∑ipiviki⊤- (0)≈ _ip_i\,v_ik_i , so gradient descent learns each pair at a rate proportional to pip_i. With orthonormal embeddings, this sum is a singular value decomposition with singular values pip_i, and msign⁡(∑ipiviki⊤)=∑iviki⊤ ( _ip_i\,v_ik_i )= _iv_ik_i gives every pair the same weight. wang2026muon; li2026muon analyze orthonormal embeddings, and kim2026sharp extend the result to Gaussian embeddings. Appendix B Appendix: Proofs in section 2 Sign convention. Throughout the paper, LMO means linear minimization oracle: lmoK⁡(G)∈arg⁡minX∈K⟨G,X⟩.lmo_K(G)∈ _X∈ K G,X . We write N=Nk∘⋯∘N1N=N_k ·s N_1 for a chain of normalization maps and T=−NT=-N for the signed update used in minimization. For the norm balls and odd oracle selections considered here, lmoj=−Njlmo_D_j=-N_j, so T=(−1)k+1(lmok∘⋯∘lmo1)=−N,Xt+1=Xt+ηtT(Mt)=Xt−ηtN(Mt).T=(-1)^k+1 (lmo_D_k ·s _D_1 )=-N, X_t+1=X_t+ _tT(M_t)=X_t- _tN(M_t). For any fixed set C, with −C:=−X:X∈C-C:=\-X:X∈ C\, N(G)∈argmaxZ∈C⟨G,Z⟩⟺−N(G)∈argminX∈−C⟨G,X⟩.N(G)∈ _Z∈ C G,Z -N(G)∈ _X∈-C G,X . Thus the non-LMO proofs below test whether N is a linear maximization oracle; their conclusions concern the signed update T=−NT=-N as a linear minimization oracle. The reflected set is the same set only when C=−C=-C, as for the norm balls used here. We write msign⁡(G)=UV⊤msign(G)=UV for the canonical polar factor of G=UΣV⊤G=U V , with zero singular directions mapped to zero, and SgnSgn for the entrywise sign with Sgn⁡(0)=0Sgn(0)=0. Row and column normalization are denoted by RNRN and CNCN, respectively. The proofs use two necessary conditions for a linear maximization oracle. B.1 Non-LMO results Lemma B.1 (Lemma 1 of the paper). Let N⁡(G)∈arg⁡maxX∈K⟨G,X⟩N(G)∈ _X∈ K G,X for every G, for a fixed set K, and let Y=N⁡(G)Y=N(G). If ‖N⁡(Y)‖F=‖Y‖F\|N(Y)\|_F=\|Y\|_F, then N⁡(Y)=YN(Y)=Y. Proof. Since every oracle output belongs to K, Y∈KY∈ K. Optimality at the input Y gives ⟨Y,N⁡(Y)⟩≥⟨Y,Y⟩=‖Y‖F2. Y,N(Y) ≥ Y,Y =\|Y\|_F^2. By Cauchy–Schwarz and the assumed equality of norms, ⟨Y,N⁡(Y)⟩≤‖Y‖F‖N⁡(Y)‖F=‖Y‖F2. Y,N(Y) ≤\|Y\|_F\|N(Y)\|_F=\|Y\|_F^2. Hence equality holds in Cauchy–Schwarz. The equal norms imply N⁡(Y)=YN(Y)=Y, including the case Y=0Y=0. ∎ Lemma B.2 (Lemma 2 of the paper). If N⁡(G)∈arg⁡maxX∈K⟨G,X⟩N(G)∈ _X∈ K G,X for every G, for a fixed set K, and N⁡(−G)=−N⁡(G)N(-G)=-N(G), then ⟨G,N⁡(G)⟩≥0. G,N(G) ≥ 0. Proof. Because N⁡(−G)=−N⁡(G)N(-G)=-N(G) is itself an output of the oracle, it belongs to K. Therefore optimality at G implies ⟨G,N⁡(G)⟩≥⟨G,−N⁡(G)⟩=−⟨G,N⁡(G)⟩. G,N(G) ≥ G,-N(G) =- G,N(G) . Thus 2⟨G,N⁡(G)⟩≥02 G,N(G) ≥ 0. ∎ Consequently, to rule out a linear maximization-oracle representation for N, it suffices to find G such that N⁡(N⁡(G))≠N⁡(G)N(N(G))≠ N(G) but ‖N⁡(N⁡(G))‖F=‖N⁡(G)‖F\|N(N(G))\|_F=\|N(G)\|_F, or, for an odd map, to find G such that ⟨G,N⁡(G)⟩<0 G,N(G) <0. By the sign convention above, either argument rules out an LMO representation for the signed update −N-N. Lemma B.3 (Tall matrix-sign and row-normalization witnesses). Assume m>nm>n. Each normalization map N listed in parts (i) and (i) of Theorem 1 has an input G for which N⁡(N⁡(G))≠N⁡(G)N(N(G))≠ N(G) and ‖N⁡(N⁡(G))‖F=‖N⁡(G)‖F\|N(N(G))\|_F=\|N(G)\|_F. Proof. Assume throughout this subsection that m>nm>n, and put r=m−n≥1.r=m-n≥ 1. For 0<t<10<t<1, define A⁡(t)∈ℝm×nA(t) ^m× n by the rows A(t)1,: A(t)_1,: =(1−t,t,0,…,0), =( 1-t, t,0,…,0), A(t)2,: A(t)_2,: =(1−t,−t,0,…,0), =( 1-t,- t,0,…,0), followed by r copies of e2⊤e_2 , and then e3⊤,…,en⊤.e_3 ,…,e_n . Thus A⁡(t)A(t) has exactly 2+r+(n−2)=m2+r+(n-2)=m rows. Every row has Euclidean norm one, so RN⁡(A⁡(t))=A⁡(t).RN(A(t))=A(t). Moreover, A(t)⊤A(t)=diag⁡(2(1−t),r+2t, 1,…,1).A(t) A(t)= (2(1-t),\,r+2t,\,1,…,1 ). Define Q(t):=msign(A(t))=A(t)(A(t)⊤A(t))−1/2.Q(t):=msign(A(t))=A(t) (A(t) A(t) )^-1/2. Since the columns of A⁡(t)A(t) are orthogonal, column normalization produces exactly the same matrix: CN⁡(A⁡(t))=Q⁡(t).CN(A(t))=Q(t). Since Q⁡(t)Q(t) has orthonormal columns, msign⁡(Q⁡(t))=Q⁡(t),CN⁡(Q⁡(t))=Q⁡(t).msign(Q(t))=Q(t),CN(Q(t))=Q(t). The only nontrivial transition is row normalization of Q⁡(t)Q(t). For either of its first two rows, the squared second coordinate after normalization is t/(r+2t)12+t/(r+2t)=2tr+4t. t/(r+2t) 12+t/(r+2t)= 2tr+4t. Hence, with ϕ⁡(t):=2tr+4t,φ(t):= 2tr+4t, we have the exact identities RN⁡(A⁡(t)) (A(t)) =A⁡(t), =A(t), (5) CN⁡(A⁡(t)) (A(t)) =msign⁡(A⁡(t))=Q⁡(t), =msign(A(t))=Q(t), CN⁡(Q⁡(t)) (Q(t)) =msign⁡(Q⁡(t))=Q⁡(t), =msign(Q(t))=Q(t), RN⁡(Q⁡(t)) (Q(t)) =A⁡(ϕ⁡(t)). =A(φ(t)). The map ϕφ is strictly increasing: ϕ′(t)=2r(r+4t)2>0.φ (t)= 2r(r+4t)^2>0. At t0=12,t_0= 12, we have ϕ⁡(t0)=1r+2<12=t0.φ(t_0)= 1r+2< 12=t_0. Strict monotonicity therefore gives t0>ϕ⁡(t0)>ϕ2(t0)>⋯t_0>φ(t_0)>φ^2(t_0)>·s for every finite number of iterations. We now apply this one-dimensional dynamics to the maps in parts (i)–(i) of Theorem 1. For part (i), equation 5 gives the following transitions. The map N=msign∘RNN=msign sends A⁡(t)A(t) to Q⁡(t)Q(t), and a second application sends it to Q⁡(ϕ⁡(t))Q(φ(t)). The map RN∘msignRN sends A⁡(t)A(t) to A⁡(ϕ⁡(t))A(φ(t)), then to A(ϕ2(t))A(φ^2(t)). The maps msign∘RN∘CNmsign and CN∘RN∘msignCN both send A⁡(t)A(t) to Q⁡(ϕ⁡(t))Q(φ(t)), then to Q(ϕ2(t))Q(φ^2(t)). Finally, RN∘CN∘msignRN has the same transitions as RN∘msignRN on this family, because CN⁡(Q⁡(t))=Q⁡(t)CN(Q(t))=Q(t). All five maps are therefore nonidempotent at t=t0t=t_0. Along these witness orbits, outputs have Frobenius norm n n when the final operation is msignmsign or CNCN, and m m when it is RNRN. Lemma B.1 excludes a fixed-set linear maximization-oracle representation for each normalization map, and hence an LMO representation for its negative. For part (i), let F=msign∘RNF=msign . Then F⁡(A⁡(t))=Q⁡(t)F(A(t))=Q(t) and F⁡(Q⁡(t))=Q⁡(ϕ⁡(t))F(Q(t))=Q(φ(t)), so for every k≥1k≥ 1, Fk(A⁡(t))=Q⁡(ϕk−1(t))F^k(A(t))=Q(φ^k-1(t)). Applying the same map FkF^k once more replaces ϕk−1(t)φ^k-1(t) by ϕ2k−1(t)φ^2k-1(t), which is distinct at t=t0t=t_0. For the reverse composition H=RN∘msignH=RN , the same argument gives the distinct iterates A⁡(ϕk(t0))A(φ^k(t_0)) and A⁡(ϕ2k(t0))A(φ^2k(t_0)). Again the outputs along these witness orbits have constant Frobenius norm, so Lemma B.1 proves part (i). ∎ Lemma B.4 (Finite mixed row/column-normalization words). Every finite composition of RNRN and CNCN that contains both is not a linear maximization oracle of a fixed set K. Proof. Repeated adjacent copies of the same normalization may be removed because RN2=RN,CN2=CN.RN^2=RN,CN^2=CN. Thus it suffices to consider an alternating sequence of RNRN and CNCN. Put r=m−1,s=n−1,r=m-1,s=n-1, and fix any constant a>1a>1. For t>0t>0, define the positive row-normalized matrix R⁡(t)R(t) by its squared entries: R(t)ij2=atat+s,i=1,j=1,1at+s,i=1,j≥2,t+s,i≥2,j=1,1t+s,i≥2,j≥2.R(t)_ij^2= cases atat+s,&i=1,\ j=1,\\[5.69054pt] 1at+s,&i=1,\ j≥ 2,\\[5.69054pt] tt+s,&i≥ 2,\ j=1,\\[5.69054pt] 1t+s,&i≥ 2,\ j≥ 2. cases Let C1(t)C_1(t) denote the squared norm of the first column and C2(t)C_2(t) the squared norm of any other column: C1(t)=atat+s+rt+s,C_1(t)= atat+s+r tt+s, and C2(t)=1at+s+r1t+s.C_2(t)= 1at+s+r 1t+s. Column normalization multiplies the squared first-to-other coordinate ratio in every row by C2(t)/C1(t)C_2(t)/C_1(t). Row normalization preserves these within-row ratios. Therefore RN⁡(CN⁡(R⁡(t)))=R⁡(F⁡(t)),RN(CN(R(t)))=R(F(t)), where F⁡(t)=tC2(t)C1(t)=(1+ra)t+msamt+s⁡(a+r).F(t)=t C_2(t)C_1(t)= (1+ra)t+msam\,t+s(a+r). This scalar map is strictly increasing. Indeed, F′(t)=rs(a−1)2(amt+s⁡(a+r))2>0.F (t)= rs(a-1)^2 (am\,t+s(a+r) )^2>0. The fixed-point equation F⁡(t)=tF(t)=t is a quadratic with positive leading coefficient and negative constant term: amt2+[s⁡(a+r)−(1+ra)]t−ms=0.am\,t^2+ [s(a+r)-(1+ra) ]t-ms=0. It therefore has exactly one positive root, denoted t⋆t_ . Moreover, F⁡(1)−1=n⁡(1−a)am+s⁡(a+r)<0.F(1)-1= n(1-a)am+s(a+r)<0. Thus 1>t⋆1>t_ . Since F is increasing, t⋆<F⁡(t)<t⁡(t>t⋆).t_ <F(t)<t\,(t>t_ ). Consequently 1>F⁡(1)>F2(1)>⋯>t⋆,1>F(1)>F^2(1)>·s>t_ , and no two finite iterates coincide. Suppose first that a reduced chain of RNRN and CNCN ends in RNRN. Starting from R⁡(t)R(t), every effective CNCN-then-RNRN pair applies one copy of F; any initial RNRN is inert. Hence for some integer k≥1k≥ 1, N⁡(R⁡(t))=R⁡(Fk(t)).N(R(t))=R(F^k(t)). Applying the same word again gives N⁡(N⁡(R⁡(t)))=R⁡(F2k(t)).N(N(R(t)))=R(F^2k(t)). At t=1t=1 these are distinct. If instead the word ends in CNCN, transpose the entire construction. Transposition exchanges RNRN and CNCN, and ⟨G,X⟩=⟨G⊤,X⊤⟩. G,X = G ,X . Thus a linear maximization-oracle representation for the original word would induce one for the transposed RN-ending word, which we just ruled out. Finally, on the positive witness family, an RN-ending word has output norm m m, while a CN-ending word has output norm n n. Lemma B.1 therefore proves that every finite composition containing both RNRN and CNCN is not a linear maximization oracle of any fixed set. ∎ The limiting row/column equilibration map ℰE uses the compatible normalization ∥ℰ(G)i,:∥22=n/m\|E(G)_i,:\|_2^2=n/m and ∥ℰ(G):,j∥22=1\|E(G)_:,j\|_2^2=1, so every equilibrated matrix has squared Frobenius norm n. The following lemma identifies the balanced representative inside a positive diagonal orbit. Lemma B.5 (Uniqueness in a positive diagonal orbit). Let B,CB,C be strictly positive matrices in the same positive diagonal orbit: C=D1BD2,D1,D2≻0.C=D_1BD_2,D_1,D_2 0. If B and C have the same squared row sums and the same squared column sums, then B=CB=C. Proof. Write Sij=Bij2,Uij=Cij2.S_ij=B_ij^2,U_ij=C_ij^2. For suitable scalars αi,βj _i, _j, log⁡UijSij=αi+βj. U_ijS_ij= _i+ _j. Therefore ∑ij(Uij−Sij)log⁡UijSij=∑iαi∑j(Uij−Sij)+∑jβj∑i(Uij−Sij)=0, _ij(U_ij-S_ij) U_ijS_ij= _i _i _j(U_ij-S_ij)+ _j _j _i(U_ij-S_ij)=0, because the row and column marginals agree. But each summand satisfies (u−s)log⁡(u/s)≥0(u-s) (u/s)≥ 0 for u,s>0u,s>0, with equality only when u=su=s. Hence U=SU=S, and positivity gives B=CB=C. ∎ Lemma B.6 (The limiting equilibration map is not a linear maximization oracle). The limiting equilibration map E in Theorem 1 - (i) is not a linear maximization oracle of any fixed set. Proof. If B is already balanced and L,R≻0L,R 0 are diagonal, then ℰ⁡(LBR)=B.E(LBR)=B. (6) We now construct a positive balanced matrix explicitly. First suppose m≥nm≥ n. Set b=12m,a=n+12m,c=1m.b= 1 2m,a= n+12m,c= 1 m. Let the first n rows form the n×n× n matrix (a−b)In+b⊤,(a-b)I_n+b11 , and let each of the remaining m−nm-n rows equal c⊤.c1 . Call the resulting matrix B. For a top row, a2+(n−1)b2=n+12m+n−12m=nm.a^2+(n-1)b^2= n+12m+ n-12m= nm. For a lower row, nc2=nm.nc^2= nm. Each column therefore has squared norm nm+(m−n)c2=nm+m−nm=1. nm+(m-n)c^2= nm+ m-nm=1. Thus B is balanced. When m<nm<n, use Bm,n=nmBn,m⊤,B_m,n= nm\,B_n,m , which has squared row norms n/mn/m and squared column norms one. We prove that E is not a linear maximization oracle by a local perturbation inside the balanced manifold. Let L=diag⁡(2,1,…,1),R=diag⁡(2,1,…,1),L= (2,1,…,1),R= (2,1,…,1), and define G=LBR.G=LBR. By equation 6, E⁡(G)=B.E(G)=B. Let Z be supported on the leading 2×22× 2 block with Z11=Z22=1,Z12=Z21=−1.Z_11=Z_22=1,\,Z_12=Z_21=-1. Its row sums and column sums are zero. For sufficiently small t, define B⁡(t)B(t) entrywise by B(t)ij=Bij2+tZij.B(t)_ij= B_ij^2+tZ_ij. All entries remain positive. Since Z has zero row and column sums, B⁡(t)B(t) has exactly the same squared row and column norms as B, so B⁡(t)B(t) is balanced. Hence E⁡(B⁡(t))=B⁡(t),E(B(t))=B(t), so B⁡(t)B(t) is itself an output of E. On the other hand, dt⟨G,B⁡(t)⟩|t=0=12∑ijLiiRjjZij. . ddt G,B(t) |_t=0= 12 _ijL_iR_jZ_ij. Only the first 2×22× 2 block contributes, giving dt⟨G,B⁡(t)⟩|t=0=12(2−1)(2−1)=12>0. . ddt G,B(t) |_t=0= 12(2-1)(2-1)= 12>0. Thus for all sufficiently small t>0t>0, ⟨G,B⁡(t)⟩>⟨G,B⟩=⟨G,E⁡(G)⟩. G,B(t) > G,B = G,E(G) . If E were a linear maximization oracle of a fixed set K, then every output B⁡(t)B(t) would belong to K, contradicting optimality of E⁡(G)=BE(G)=B at the input G. Therefore the limiting equilibration E is not a linear maximization oracle of any fixed set. ∎ Now we find explicit counterexamples of Lemma 2 for msign∘ℰ \ E and the iteration of sign and msign , proving these chains are not linear maximization oracles. Lemma B.7 (Negative alignment after equilibration). For every shape covered by Theorem 1, the odd map msign∘ℰmsign admits an input G such that ⟨G,(msign∘ℰ)(G)⟩<0. G,(msign )(G) <0. Hence msign∘ℰmsign is not a linear maximization oracle of any fixed set. Proof. Write N=msign∘ℰN=msign . Both constituents are odd, so N is odd. By Lemma B.2, it is therefore enough to find G with ⟨G,msign⁡(ℰ⁡(G))⟩<0. G,msign(E(G)) <0. First suppose m>nm>n, and use the balanced positive matrix B from the previous subsection. Put r=m−n,s=a+(n−1)b,λ=s2+nrc2.r=m-n,\,s=a+(n-1)b,\,λ=s^2+nrc^2. The leading n×n× n block of B is (a−b)In+bJn.(a-b)I_n+bJ_n. A direct calculation gives B⊤B=(a−b)2(I−1nJn)+λ1nJn.B B=(a-b)^2 (I- 1nJ_n )+λ\, 1nJ_n. Hence, writing Q=msign⁡(B),Q=msign(B), the leading n×n× n block of Q is Qtop=I−1nJn+sλ1nJn.Q_ top=I- 1nJ_n+ s λ 1nJ_n. (7) Since r>0r>0, λ>s2,λ>s^2, and therefore every off-diagonal entry of the block equation 7 is strictly negative: (Qtop)ij=1n(sλ−1)<0,i≠j.(Q_ top)_ij= 1n ( s λ-1 )<0,\,i≠ j. Choose, say, the entry (1,2)(1,2), and set δ=−B12Q12>0.δ=-B_12Q_12>0. For a scalar M≥1M≥ 1, let LM=diag⁡(M,1,…,1),RM=diag⁡(1,M,1,…,1),L_M= (M,1,…,1),\,R_M= (1,M,1,…,1), and define GM=LMBRM.G_M=L_MBR_M. By diagonal invariance, E⁡(GM)=B,E(G_M)=B, so N⁡(GM)=Q.N(G_M)=Q. The selected (1,2)(1,2) contribution to ⟨GM,Q⟩ G_M,Q is −δM2-δ M^2. Every other entry receives weight at most M. Moreover ∑ij|BijQij|≤‖B‖F‖Q‖F=n, _ij|B_ijQ_ij|≤\|B\|_F\|Q\|_F=n, because B has squared Frobenius norm n and Q has n orthonormal columns. Therefore ⟨GM,Q⟩≤−δM2+nM. G_M,Q ≤-δ M^2+nM. Taking M>nδM> nδ gives ⟨GM,N⁡(GM)⟩<0. G_M,N(G_M) <0. For m<nm<n, transpose the corresponding n×mn× m witness. Compatible equilibration obeys the transpose relation up to the positive scalar n/m n/m, and the polar factor removes this scalar. Thus the same negative inner product is obtained. Now suppose m=n=q≥3.m=n=q≥ 3. Consider first the 3×33× 3 positive circulant matrix B3=114(123312231).B_3= 1 14 pmatrix1&2&3\\ 3&1&2\\ 2&3&1 pmatrix. Every row and every column has norm one. Furthermore, B3⊤B3=314I3+1114J3,B_3 B_3= 314I_3+ 1114J_3, and direct evaluation of the polar factor gives Q3:=msign⁡(B3)=13(1−311+31+31−3111+31−3).Q_3:=msign(B_3)= 13 pmatrix1- 3&1&1+ 3\\ 1+ 3&1- 3&1\\ 1&1+ 3&1- 3 pmatrix. In particular, (Q3)11=1−33<0.(Q_3)_11= 1- 33<0. For q>3q>3, let B0=B3⊕Iq−3.B_0=B_3 I_q-3. This matrix is balanced and invertible, but has zero entries. To obtain a strictly positive witness, define the squared-entry matrix Sε=(1−ε)(B0∘B0)+εJqq,S_ =(1- )(B_0 B_0)+ J_qq, where ∘ denotes the Hadamard product, and set (Bε)ij=(Sε)ij.(B_ )_ij= (S_ )_ij. Every row and column of SεS_ sums to one, so BεB_ is balanced. It is strictly positive. Since Bε→B0B_ → B_0 and B0B_0 is invertible, BεB_ is invertible for all sufficiently small ε>0 >0. Continuity of the polar factor at invertible matrices gives msign⁡(Bε)⟶Q3⊕Iq−3.msign(B_ ) Q_3 I_q-3. Hence, for sufficiently small fixed ε , (Bε)11>0,[msign⁡(Bε)]11<0.(B_ )_11>0,\,[msign(B_ )]_11<0. We can therefore repeat the preceding diagonal-amplification argument, now scaling the first row and first column. For sufficiently large M, GM=LMBεLMG_M=L_MB_ L_M satisfies E⁡(GM)=BεE(G_M)=B_ and ⟨GM,msign⁡(Bε)⟩<0. G_M,msign(B_ ) <0. Thus msign∘Emsign E has a negative-alignment witness in every shape m,n≥5m,n≥ 5. Since it is odd, Lemma B.2 proves that it is not a linear maximization oracle of any fixed set. ∎ Lemma B.8 (Mixed matrix-sign and entrywise-sign words). Every finite composition of msignmsign and SgnSgn that contains both admits an input G with negative alignment, ⟨G,N⁡(G)⟩<0. G,N(G) <0. Hence no such composition is a linear maximization oracle of any fixed set. Proof. Write P=msignP=msign and S=SgnS=Sgn. Since both maps are idempotent, every finite chain containing both reduces to an alternating chain of length at least two. Both maps are odd, so by Lemma B.2 it suffices to construct a negative-alignment witness for every alternating chain. We first consider square matrices of size q≥5q≥ 5. Suppose the first applied operation is S. Put k=q−1k=q-1, J=⊤J=11 , B=J−2IB=J-2I, and A=B+2e1e1⊤A=B+2e_1e_1 . Thus A has +1+1 in its (1,1)(1,1) entry, −1-1 on the remaining diagonal, and +1+1 off the diagonal. Let Q=P⁡(A)Q=P(A). To determine its signs, let u=k−1/2(0,1,…,1)⊤u=k^-1/2(0,1,…,1) and decompose ℝq=span⁡e1,u⊕e1,u⟂R^q=span\e_1,u\ \e_1,u\ . On the orthogonal complement, A=−2IA=-2I, while on span⁡e1,uspan\e_1,u\ it is represented by HA=(1k−2).H_A= pmatrix1& k\\ k&k-2 pmatrix. Since detHA=−2 H_A=-2, its eigenvalues have opposite signs. Setting D=k2−2k+9D= k^2-2k+9, its matrix sign is msign⁡(HA)=1D(3−k2k2k−3).msign(H_A)= 1D pmatrix3-k&2 k\\ 2 k&k-3 pmatrix. It follows that Q has negative diagonal entries and positive off-diagonal entries for k≥4k≥ 4. Hence S⁡(Q)=BS(Q)=B. The eigenvalues of B are q−2q-2 along 1 and −2-2 on ⟂1 , so P⁡(B)=2J/q−I=:RP(B)=2J/q-I=:R. For q≥5q≥ 5, R again has negative diagonal and positive off-diagonal entries, hence S⁡(R)=BS(R)=B. Thus every alternating chain beginning with S outputs Q, B, or R. Now take Gt=B+(t+1)e1e1⊤G_t=B+(t+1)e_1e_1 . Since S⁡(Gt)=AS(G_t)=A, the possible chain outputs remain Q,B,RQ,B,R, independently of t. Their (1,1)(1,1) entries are all strictly negative. Therefore ⟨Gt,X⟩=⟨B,X⟩+(t+1)X11→−∞ G_t,X = B,X +(t+1)X_11→-∞ for every X∈Q,B,RX∈\Q,B,R\. Taking t sufficiently large gives a negative-alignment witness for every chain beginning with S. Now suppose the first applied operation is P. Let w=(8q,8q,1,…,1)⊤w=(8q,8q,1,…,1) , let dw=‖w‖22d_w=\|w\|_2^2, and define the Householder reflection Q=I−2ww⊤/dwQ=I-2w /d_w. Its diagonal entries are positive and its off-diagonal entries are negative, so S⁡(Q)=2I−J=:BS(Q)=2I-J=:B. Since B has eigenvalue 2−q2-q along 1 and eigenvalue 22 on ⟂1 , we have P⁡(B)=I−2J/q=:RP(B)=I-2J/q=:R, and again S⁡(R)=BS(R)=B. Thus every such alternating chain eventually alternates between B and R. It remains to choose an input with polar factor Q and negative alignment with both B and R. Define z=(I−ww⊤/dw)z=(I-w /d_w)1. Then z⟂wz w, so Qz=zQz=z. Moreover, if r=‖z‖22r=\|z\|_2^2, then r=q−(17q−2)2128q2+q−2>q2r=q- (17q-2)^2128q^2+q-2> q2 for q≥5q≥ 5. Consequently z⊤Bz=2r−r2<0z Bz=2r-r^2<0 and z⊤Rz=r−2r2/q<0z Rz=r-2r^2/q<0. For γ>0γ>0, set Gγ=Q⁡(I+γzz⊤)G_γ=Q(I+γ z ). The second factor is positive definite, so P⁡(Gγ)=QP(G_γ)=Q. Since Qz=zQz=z, we also have Gγ=Q+γzz⊤G_γ=Q+γ z . Therefore both ⟨Gγ,B⟩ G_γ,B and ⟨Gγ,R⟩ G_γ,R tend to −∞-∞ as γ→∞γ→∞. Taking γ sufficiently large gives a negative-alignment witness for every chain beginning with P. Finally, for rectangular matrices, take q=min⁡m,nq= \m,n\ and zero-pad the corresponding q×q× q witness. Both P and S commute with zero padding, and Frobenius inner products are unchanged. Lemma B.2 therefore rules out a linear maximization-oracle representation for every finite composition containing both msignmsign and SgnSgn. ∎ Theorem B.9 (Theorem 1 of the paper). Let m,n≥5m,n≥ 5, and let normalizations map zero rows and columns to zero. For each map N:ℝm×n→ℝm×nN:R^m× n ^m× n listed below, there is no fixed set K⊆ℝm×nK ^m× n such that N⁡(G)∈arg⁡maxX∈K⟨G,X⟩N(G)∈ _X∈ K G,X for every G: (i) msign∘RNmsign , msign∘RN∘CNmsign , RN∘msignRN , CN∘RN∘msignCN , and RN∘CN∘msignRN for m>nm>n. (i) Alternation with matrix sign: (msign∘RN)∘k(msign ) k and (RN∘msign)∘k(RN ) k for every k≥1k≥ 1 and m>nm>n. (i) Every finite composition of RNRN and CNCN that contains both, in particular (RN∘CN)∘k(RN ) k; its limit ℰE; and msign∘ℰmsign , with the last two considered on matrices with no zero entries. (iv) Every finite composition of msignmsign and SgnSgn that contains both. Equivalently, for every listed normalization map N, the signed update −N-N is not an LMO of any fixed set. Proof. Parts (i) and (i) follow from Lemma B.3 and Lemma B.1. The finite-word claim in part (i) follows from Lemma B.4 and Lemma B.1; the limit E is excluded by Lemma B.6, and msign∘Emsign E is excluded by Lemmas B.7 and B.2. Part (iv) follows from Lemmas B.8 and B.2. ∎ B.2 Convergence and non-convergence guarantees Proposition B.10 (Proposition 1 of the paper). Let T be positively scale invariant with ⟨G,T⁡(G)⟩>0 G,T(G) >0 for some G. For the convex, 11-smooth quadratic f⁡(X)=⟨G,X⟩22‖G‖F2f(X)= G,X ^22\|G\|_F^2 with X0=G/‖G‖F2X_0=G/\|G\|_F^2, the iterates Xt+1=Xt+ηtT(∇f(Xt))X_t+1=X_t+ _tT(∇ f(X_t)) satisfy f⁡(Xt)→∞f(X_t)→∞ for every nonnegative step-size sequence with ∑tηt=∞ _t _t=∞. Proof. Let a:=⟨G,T⁡(G)⟩>0a:= G,T(G) >0. The Hessian of f is the rank-one orthogonal projector H⁡[Z]=G⁡⟨G,Z⟩/‖G‖F2H[Z]=G G,Z /\|G\|_F^2, so f is convex, nonnegative, and 11-smooth. Put ct=⟨G,Xt⟩c_t= G,X_t . Since c0=1c_0=1, whenever ct>0c_t>0, ∇f(Xt)=(ct/‖G‖F2)G∇ f(X_t)=(c_t/\|G\|_F^2)G is a positive multiple of G. Positive scale invariance gives T⁡(∇f(Xt))=T⁡(G)T(∇ f(X_t))=T(G). Thus ct+1=ct+ηtac_t+1=c_t+ _ta, and induction yields ct=1+a∑j<tηj>0,Xt=X0+(∑j<tηj)T(G).c_t=1+a _j<t _j>0, X_t=X_0+ ( _j<t _j )T(G). Consequently, f⁡(Xt)=(1+a∑j<tηj)22‖G‖F2⟶∞.f(X_t)= (1+a _j<t _j )^22\|G\|_F^2 ∞. ∎ The witnesses in Lemmas B.8 and B.7 have ⟨G,N⁡(G)⟩<0 G,N(G) <0 for the normalization map N. Their signed updates T=−NT=-N therefore satisfy ⟨G,T⁡(G)⟩>0 G,T(G) >0, as required by Proposition B.10. Equivalently, the stateless iteration Xt+1=Xt−ηtN(∇f(Xt))X_t+1=X_t- _tN(∇ f(X_t)) diverges. The Aurora claim follows by the same sign conversion from the next lemma. Lemma B.11 (Negative alignment for two Aurora cycles). If n=2n=2 and m≥258m≥ 258, then (msign∘RN)∘2(msign ) 2 has a full-column-rank input G with nonzero entries such that ⟨G,(msign∘RN)∘2(G)⟩<0. G,(msign ) 2(G) <0. Proof. Let m=r+2m=r+2 with r≥256r≥ 256, put s=5+1/rs= 5+1/r, and define q=1s(r−1,…,r−1⏟r entries,−2,1)⊤,Q=[m−1/2q].q= 1s ( r^-1,…,r^-1_r entries,-2,1 ) ,\,Q= bmatrixm^-1/21&q bmatrix. Then Tq=01^Tq=0 and ‖q‖=1\|q\|=1, so QTQ=I2Q^TQ=I_2. Write c=m/sc= m/s and a=c/ra=c/r. The rows of RN⁡(Q)RN(Q) have slopes a,−2c,ca,-2c,c, with the first slope repeated r times. Let C=RN(Q)TRN(Q)C=RN(Q)^TRN(Q). Direct summation gives C22<3C_22<3, C11≤mC_11≤ m, and cC12≥c21+a2−12=(r+2)r25r2+2r+2−12≥r6−12.cC_12≥ c^21+a^2- 12= (r+2)r^25r^2+2r+2- 12≥ r6- 12. For h=detCh= C, the inverse square root of a positive definite 2×22× 2 matrix is a positive scalar multiple of (C22+h−C12−C12C11+h) ( smallmatrixC_22+h&-C_12\\ -C_12&C_11+h smallmatrix ). Hence the first coordinate of the last row of V:=msign⁡(RN⁡(Q))V:=msign(RN(Q)) has the sign of C22+h−cC12C_22+h-cC_12. Since h≤3mh≤ 3m and, for r≥256r≥ 256, r/6−1/2>3+3(r+2)r/6-1/2>3+ 3(r+2), this coordinate is strictly negative. Now let Aε=RN⁡([,εq])A_ =RN([1, q]). A QR expansion gives msign⁡(Aε)→Qmsign(A_ )→ Q as ε↓0 0: the second QR diagonal is ε+O⁡(ε3) +O( ^3), the off-diagonal term is O⁡(ε3)O( ^3), and the polar factor of the resulting upper-triangular 2×22× 2 matrix tends to the identity. Therefore (msign∘RN)2([,εq])⟶V.(msign )^2([1, q]) V. For sufficiently small fixed ε>0 >0, the inner product between the last row of AεA_ and the last row of this output is negative. Multiplying only the last input row by a sufficiently large positive scalar leaves the first row-normalization, and hence the chained output, unchanged, while making the total alignment negative. The input remains full column rank with nonzero entries. This proves the claimed negative-alignment example for (msign∘RN)∘2(msign ) 2 when n=2n=2 and m≥258m≥ 258. ∎ Definition B.12 (Definition 1 of the paper). Let K be a compact set containing the origin and define mK(G):=minX∈K⁡⟨G,X⟩.m_K(G):= _X∈ K G,X . For ν∈(0,1]ν∈(0,1], a map U is a ν-approximate LMO of K if, for every G, U⁡(G)∈K,⟨G,U⁡(G)⟩≤νmK(G).U(G)∈ K, G,U(G) ≤ν m_K(G). Write hK(G):=maxX∈K⁡⟨G,X⟩h_K(G):= _X∈ K G,X . The nonnegative minimization gap is −mK(G)=hK(−G)=h−K(G)-m_K(G)=h_K(-G)=h_-K(G). For an origin-symmetric set K=−K=-K, this equals hK(G)h_K(G), the notation used for the norm-ball bounds in the main text. Proposition B.13 (Proposition 2 of the paper). Let f have an L-Lipschitz gradient, let Δ=f⁡(X0)−inff<∞ =f(X_0)- f<∞, and let U be a ν-approximate LMO of a compact set K containing the origin. Write R:=maxX∈K⁡‖X‖FR:= _X∈ K\|X\|_F. (i) Gradient descent. With Xt+1=Xt+ηU(∇f(Xt))X_t+1=X_t+η U(∇ f(X_t)) and η=2ΔLR2T,1T∑t=0T−1[−mK(∇f(Xt))]≤Rν2LΔT.η= 2 LR^2T, 1T _t=0^T-1[-m_K(∇ f(X_t))]≤ Rν 2L T. (i) SGD with momentum. With unbiased stochastic gradients of conditional variance at most σ2σ^2, Mt=βMt−1+(1−β)g^t,Xt+1=Xt+ηU(Mt),M_t=β M_t-1+(1-β) g_t, X_t+1=X_t+η U(M_t), 1−β=T−1/21-β=T^-1/2, η=cT−3/4η=cT^-3/4, and fixed initialization with finite initial tracking error, the same average in expectation is O(ν−1[Δ/c+R(LRc+σ)]T−1/4).O\! (ν^-1 [ /c+R(LRc+σ) ]T^-1/4 ). For K=−K=-K, both statements can equivalently be written with hK(∇f(Xt))h_K(∇ f(X_t)) in place of −mK(∇f(Xt))-m_K(∇ f(X_t)). Proof. Feasibility implies ‖U⁡(v)‖F≤R\|U(v)\|_F≤ R. The minimum-value function is R-Lipschitz: |mK(g)−mK(v)|≤R‖g−v‖F m_K(g)-m_K(v) ≤ R\|g-v\|_F. The robust descent-alignment inequality is ⟨g,U⁡(v)⟩≤νmK(g)+(1+ν)R‖g−v‖F. g,U(v) ≤ν m_K(g)+(1+ν)R\|g-v\|_F. (8) Indeed, ⟨g,U⁡(v)⟩=⟨v,U⁡(v)⟩+⟨g−v,U⁡(v)⟩≤νmK(v)+R‖g−v‖F, g,U(v) = v,U(v) + g-v,U(v) ≤ν m_K(v)+R\|g-v\|_F, and the Lipschitz bound for mKm_K gives equation 8. For part (i), let gt=∇f(Xt)g_t=∇ f(X_t). Smoothness and approximate optimality give f⁡(Xt+1)≤f⁡(Xt)+νηmK(gt)+LR2η22.f(X_t+1)≤ f(X_t)+νη m_K(g_t)+ LR^2η^22. Summing and using f⁡(X0)−f⁡(XT)≤Δf(X_0)-f(X_T)≤ yields 1T∑t=0T−1[−mK(gt)]≤ΔνηT+LR2η2ν. 1T _t=0^T-1[-m_K(g_t)]≤ νη T+ LR^2η2ν. The stated choice of η balances the two terms and proves part (i). The degenerate cases R=0R=0 or Δ=0 =0 have zero minimization gap along the corresponding constant iterates and can be treated separately. For part (i), write g^t=gt+ζt g_t=g_t+ _t, where ⁡[ζt∣ℱt]=0E[ _t _t]=0 and ⁡[‖ζt‖F2∣ℱt]≤σ2E[\| _t\|_F^2 _t]≤σ^2. Let et=Mt−gte_t=M_t-g_t. Applying equation 8 with v=Mtv=M_t and then smoothness gives νη[−mK(gt)]≤f⁡(Xt)−f⁡(Xt+1)+(1+ν)Rη‖et‖F+LR2η22.νη[-m_K(g_t)]≤ f(X_t)-f(X_t+1)+(1+ν)Rη\|e_t\|_F+ LR^2η^22. Thus it remains to bound T−1∑t‖et‖FT^-1 _tE\|e_t\|_F. For t≥1t≥ 1, the tracking recursion is et=βet−1+β⁡(gt−1−gt)+(1−β)ζt.e_t=β e_t-1+β(g_t-1-g_t)+(1-β) _t. Since ‖Xt−Xt−1‖F≤Rη\|X_t-X_t-1\|_F≤ Rη, smoothness gives ‖gt−gt−1‖F≤LRη\|g_t-g_t-1\|_F≤ LRη. Unrolling the recursion, bounding the drift by the triangle inequality, and using martingale orthogonality for the noise yields ‖et‖F≤βt+1E0+β1−βLRη+σ1−β1+β,E\|e_t\|_F≤β^t+1E_0+ β1-βLRη+σ 1-β1+β, (9) where E0:=‖M−1−g0‖F<∞E_0:=E\|M_-1-g_0\|_F<∞. Averaging and summing the descent inequality gives 1T∑t<T⁡[−mK(gt)]≤1ν[ΔηT+CR(LRη1−β+σ1−β+E0(1−β)T)+LR2η2] 1T _t<TE[-m_K(g_t)]≤ 1ν [ η T+CR ( LRη1-β+σ 1-β+ E_0(1-β)T )+ LR^2η2 ] for a universal constant C. Setting 1−β=T−1/21-β=T^-1/2 and η=cT−3/4η=cT^-3/4 gives 1T∑t<T[−mK(gt)]=O(ν−1[Δ/c+R(LRc+σ)]T−1/4). 1T _t<TE[-m_K(g_t)]=O\! (ν^-1[ /c+R(LRc+σ)]T^-1/4 ). For fixed initialization, the initial-error and direct-smoothness terms have lower orders T−1/2T^-1/2 and T−3/4T^-3/4, respectively. If rF⊆KrB_F K for some r>0r>0, then −mK(g)≥r‖g‖F-m_K(g)≥ r\|g\|_F, so these bounds also control the gradient norm. Here FB_F is the Frobenius unit ball. The proof uses momentum before applying U; an EMA of a statistic internal to U is not covered. ∎ Proposition B.14 (proposition 3 of the paper). Let Krow:=X:maxi∥Xi,:∥2≤1.K_ row:=\X: _i\|X_i,:\|_2≤ 1\. Then −(RN∘msign)-(RN ) is a 1/m1/ m-approximate LMO of KrowK_ row, while −(msign∘RN)-(msign ) is a 1/m1/ m-approximate LMO of the spectral-norm unit ball. By transposition, the same holds for CNCN with 1/n1/ n. Proof. We prove the two row-normalized statements; the column-normalized versions follow by transposition. The zero input is immediate, so assume G≠0G≠ 0. First let Q=msign⁡(G)Q=msign(G) and write the polar factorization G=QHG=QH, where H=(GTG)1/2⪰0H=(G^TG)^1/2 0. Let D row-normalize the nonzero rows of Q, and set Dii=1D_i=1 on zero rows, so that RN⁡(Q)=DQRN(Q)=DQ. Since QQTQQ^T is an orthogonal projection, every row of Q has norm at most one and hence D⪰ID I. Also HQTQ=HHQ^TQ=H. Therefore ⟨G,RN⁡(Q)⟩=tr⁡(HQTDQ)=tr⁡H+tr⁡(HQT(D−I)Q)≥tr⁡H=‖G‖∗. G,RN(Q) =tr(HQ^TDQ)=trH+tr\! (HQ^T(D-I)Q ) =\|G\|_*. The support function of Krow=X:maxi∥Xi,:∥2≤1K_ row=\X: _i\|X_i,:\|_2≤ 1\ is hKrow(G)=∑i∥Gi,:∥2h_K_ row(G)= _i\|G_i,:\|_2. Since ∑i∥Gi,:∥2≤m∥G∥F≤m∥G∥∗ _i\|G_i,:\|_2≤ m\|G\|_F≤ m\|G\|_*, we obtain ⟨G,(RN∘msign)(G)⟩≥1mhKrow(G). G,(RN )(G) ≥ 1 mh_K_ row(G). Because Krow=−KrowK_ row=-K_ row, the signed output −RN⁡(Q)-RN(Q) is feasible and ⟨G,−RN⁡(Q)⟩≤−1mhKrow(G)=1mmKrow(G). G,-RN(Q) ≤- 1 mh_K_ row(G)= 1 mm_K_ row(G). This is the required approximate minimization inequality. For the reverse order, let A=RN⁡(G)A=RN(G) and write G=diag⁡(d1,…,dm)AG= (d_1,…,d_m)A, where di=∥Gi,:∥2d_i=\|G_i,:\|_2. For each nonzero row, ∥Ai,:∥2=1\|A_i,:\|_2=1. Put M=ATAM=A^TA. On range⁡(M)range(M), (M†)1/2⪰‖A‖op−1I(M )^1/2 \|A\|_op^-1I, and each Ai,:TA_i,:^T belongs to this range. Hence, with P=msign⁡(A)=A(M†)1/2P=msign(A)=A(M )^1/2, ⟨Ai,:,Pi,:⟩≥1‖A‖op≥1m A_i,:,P_i,: ≥ 1\|A\|_op≥ 1 m for every nonzero row. Zero rows have di=0d_i=0 and contribute nothing. Thus ⟨G,P⟩≥1m∑idi≥1m‖G‖∗. G,P ≥ 1 m _id_i≥ 1 m\|G\|_*. The last inequality is the nuclear-norm triangle inequality applied to the single-row decomposition of G. Write Ksp:=X:‖X‖op≤1K_ sp:=\X:\|X\|_op≤ 1\. Since ‖P‖op≤1\|P\|_op≤ 1, the signed output −P-P is feasible and ⟨G,−P⟩≤−1m‖G‖∗=1mmKsp(G). G,-P ≤- 1 m\|G\|_*= 1 mm_K_ sp(G). Thus −(msign∘RN)-(msign ) is a 1/m1/ m-approximate LMO. ∎ Proposition B.15 (Proposition 4 of the paper). Let n=2n=2. For every closed convex set K containing the image of the signed map in question, the best uniform approximation factor satisfies ν≤3(m−1)−1/4ν≤ 3(m-1)^-1/4 for −(RN∘msign)-(RN ) when m is odd, and ν≤3m−1ν≤ 3 m-1 for −(msign∘RN)-(msign ). Proof. For an unbounded K, interpret mKm_K as an infimum; if it is −∞-∞ in a tested direction, no positive approximation factor is possible. Set n=2n=2. First let N=RN∘msignN=RN and U=−NU=-N. Put r=m−1r=m-1, and assume r is even. Define p,q∈ℝmp,q ^m by p1=1−r−1/2,pi=r−3/4(i≥2),p_1= 1-r^-1/2, p_i=r^-3/4 (i≥ 2), q1=0,qi=r−1/2,2≤i≤r/2+1,−r−1/2,r/2+2≤i≤m.q_1=0, q_i= casesr^-1/2,&2≤ i≤ r/2+1,\\ -r^-1/2,&r/2+2≤ i≤ m. cases Then ‖p‖=‖q‖=1\|p\|=\|q\|=1 and p⊤q=0p q=0, so Q=[p,q]Q=[p,q] has orthonormal columns. For G=Qdiag⁡(1,ε)G=Q (1, ), ε>0 >0, we have msign⁡(G)=Qmsign(G)=Q. With a=r−3/4a=r^-3/4 and b=r−1/2b=r^-1/2, direct rowwise normalization gives ⟨G,N⁡(G)⟩=1−r−1/2+ra2+εb2a2+b2≤3 G,N(G) = 1-r^-1/2+r a^2+ b^2 a^2+b^2≤ 3 for sufficiently small ε , for example 0<ε≤r−10< ≤ r^-1. Define Q′Q by Q′1,:=e2⊤Q _1,:=e_2 and Q′i,:=r−1/2e1⊤Q _i,:=r^-1/2e_1 for i≥2i≥ 2. Then (Q′)⊤Q′=I2(Q ) Q =I_2, so N⁡(Q′)=RN⁡(Q′)N(Q )=RN(Q ), whose first row is e2⊤e_2 and remaining rows are e1⊤e_1 . Since −N⁡(Q′)∈K-N(Q )∈ K, mK(G)≤⟨G,−N⁡(Q′)⟩=−ra=−r1/4.m_K(G)≤ G,-N(Q ) =-ra=-r^1/4. If U is a ν-approximate LMO of K, then −3≤⟨G,U⁡(G)⟩≤νmK(G)≤−νr1/4,-3≤ G,U(G) ≤ν m_K(G)≤-ν r^1/4, which proves ν≤3(m−1)−1/4ν≤ 3(m-1)^-1/4. For the reverse chain, let N=msign∘RNN=msign and U=−NU=-N. Set ε=m−2 =m^-2 and define G∈ℝm×2G ^m× 2 by G1,:=e1⊤,Gi,:=εe1⊤(2≤i≤m−1),Gm,:=εe2⊤.G_1,:=e_1 , G_i,:= e_1 \ (2≤ i≤ m-1), G_m,:= e_2 . Then A=RN⁡(G)A=RN(G) has m−1m-1 copies of e1⊤e_1 and one copy of e2⊤e_2 , so A⊤A=diag⁡(m−1,1)A A= (m-1,1). Hence ⟨G,N⁡(G)⟩=1+(m−2)εm−1+ε≤3m−1. G,N(G) = 1+(m-2) m-1+ ≤ 3 m-1. Take Q′Q with first row e1⊤e_1 and all remaining rows e2⊤/m−1e_2 / m-1. It has orthonormal columns and satisfies N⁡(Q′)=Q′N(Q )=Q . Thus −Q′∈K-Q ∈ K, and mK(G)≤⟨G,−Q′⟩=−1−εm−1≤−1.m_K(G)≤ G,-Q =-1- m-1≤-1. Approximate minimization gives −3m−1≤⟨G,U⁡(G)⟩≤νmK(G)≤−ν,- 3 m-1≤ G,U(G) ≤ν m_K(G)≤-ν, so ν≤3/m−1ν≤ 3/ m-1. Both bounds vanish as m→∞m→∞, proving Proposition 4. ∎ Appendix C Appendix: Proofs in section 3 We use the model and scaling assumptions of Section 3. For clarity, we first specify the notation used throughout this appendix. Let X=[x1,…,xN]∈ℝd×N,Y=[y1,…,yN]∈ℝn×N,X=[x_1,…,x_N] ^d× N, Y=[y_1,…,y_N] ^n× N, where the entries of X and Y are independent centered Gaussians with variances 1/d1/d and 1/n1/n, respectively, and X is independent of Y. Write pi=i−α/∑j=1Nj−αp_i=i^-α/ _j=1^Nj^-α, and define PN=IN−1N⊤,s=Yp,K=PNdiag⁡(p)Y⊤∈ℝN×n.P_N=I_N- 1N11 , s=Yp, K=P_N (p)Y ^N× n. Differentiating the cross-entropy loss at zero gives the exact identity G=−∇ℒ(0)=∑i=1Npi(Σxi−1N∑j=1NΣxj)yi⊤=ΣXK.G=- (0)= _i=1^Np_i ( x_i- 1N _j=1^N x_j )y_i = XK. In particular, K⊤ei=piyi−sN,K⊤K=∑i=1Npi2yiyi⊤−1Nss⊤.K e_i=p_iy_i- sN, K K= _i=1^Np_i^2y_iy_i - 1Ns . We abbreviate the three updates by WG=G,WM=pol⁡(G),WC=pol⁡(RN⁡(G)).W_G=G, W_M=pol(G), W_C=pol(RN(G)). Here polpol is the exact canonical polar factor, equivalently msignmsign in the main text. All inverse matrix powers below are taken on the full-column-rank event, which has probability one for all sufficiently large n. The statements concern these exact updates, not a fixed number of Newton–Schulz iterations. Set h=(n/N)α,L=log⁡(enNd),s1=∑a=1dσa,s2=∑a=1dσa2.h=(n/N)^α, L= (enNd), s_1= _a=1^d _a, s_2= _a=1^d _a^2. The assumptions imply L≍log⁡nL n, N/n→∞N/n→∞, d/n→∞d/n→∞, and n−1/2+ϵ≤h≤n−(1−β+ϵ),s1≍d1−β≥cn,c≤s2≤C, n^-1/2+ε≤ h≤ n^-(1-β+ε), s_1 d^1-β≥ cn, c≤ s_2≤ C, h→0,hn1−β≤n−ϵ,L/nh→0,nNh→0, h→ 0, hn^1-β≤ n^-ε, L/nh→ 0, nNh→ 0, Lαnd→0,ns2Lα+2s1h≤Cn−ϵLα+2→0. L^α nd→ 0, ns_2\,L^α+2s_1h≤ Cn^-εL^α+2→ 0. Constants c,Cc,C may change between displays and depend only on the fixed parameters in the assumptions. A subscript b also permits dependence on a fixed b∈(0,1)b∈(0,1). The eigenvalues of a positive semidefinite matrix are indexed in nonincreasing order. Conditional OℙO_P bounds below are uniform over the specified query event; bounds for an index i≥bNi≥ bN are uniform over deterministic indices in this range. The proof is organized by conditioning: first on the query realization Y, then on the observed gradient G, and finally on the residual score randomness left after conditioning on (Y,G)(Y,G). C.1 Query-side event and population geometry Lemma C.1 (High-probability regularity of the queries). There are constants c,C>0c,C>0 such that the event ℰYE_Y defined by maxi≤N⁡‖yi‖2≤2,‖1NYY⊤‖op≤Cn, _i≤ N\|y_i\|_2≤ 2, \| 1NY \|_op≤ Cn, (10) cn−2αI⪯K⊤K⪯CI,c≤tr⁡(K⊤K)≤C, cn^-2αI K K CI, c (K K)≤ C, (11) tr((K⊤K)−1/2)≥cnα+1,∥(K⊤K)−1/2s∥2≤Cn, \! ((K K)^-1/2 )≥ cn^α+1, \|(K K)^-1/2s\|_2≤ C n, (12) λℓ(K⊤K)≥cℓ−2α(1≤ℓ≤⌊A∗L⌋), _ (K K)≥ c ^-2α (1≤ ≤ A_*L ), (13) and supB≠01N#i:‖Byi‖2≤η‖B‖Fn≤C⁡(η+nN)(0<η<1) _B≠ 0 1N\# \i:\|By_i\|_2≤η \|B\|_F n \≤ C (η+ nN ) (0<η<1) (14) satisfies ℙY(ℰYc)→0P_Y(E_Y^c)→ 0. Here A∗>0A_*>0 is a fixed sufficiently large constant, and the supremum in equation 14 is over all nonzero real matrices with n columns and any finite number of rows. The inequalities hold simultaneously for all B and η. For a deterministic index i, define s−i=∑j≠ipjyj,K−i=PN(diag⁡(p)−pieiei⊤)Y⊤.s_-i= _j≠ ip_jy_j, K_-i=P_N ( (p)-p_ie_ie_i )Y . Let ℰY,−iE_Y,-i consist of equation 11–equation 13 with K,sK,s replaced by K−i,s−iK_-i,s_-i. Thus ℰY,−iE_Y,-i is measurable with respect to Y−i=(yj)j≠iY_-i=(y_j)_j≠ i; it does not include conditions involving yiy_i. For every fixed b∈(0,1)b∈(0,1), supi≥bNℙY((ℰY∩ℰY,−i)c)→0. _i≥ bNP_Y\! ((E_Y _Y,-i)^c )→ 0. Proof. Gaussian norm and singular-value bounds give equation 10 with probability tending to one, since log⁡N=o⁡(n) N=o(n) and N/n→∞N/n→∞. Write S=∑j=1Npj2yjyj⊤,K⊤K=S−1Nss⊤.S= _j=1^Np_j^2y_jy_j , K K=S- 1Ns . To control the centering term, represent yj=εjy~jy_j= _j y_j, where the εj _j are independent Rademacher variables, independent of the Gaussian vectors y~j y_j. Conditional on Y~ Y, let A=[pjy~j]jA=[p_j y_j]_j. Then v=s⊤S−1s=ε⊤A⊤(AA⊤)−1Aε.v=s S^-1s= A (A )^-1A . The matrix between the signs is a rank-n orthogonal projection. Consequently εv=n,Varε⁡(v)=2∑j≠k[A⊤(AA⊤)−1A]jk 2≤2n.E_ v=n, _ (v)=2 _j≠ k [A (A )^-1A ]_jk^\,2≤ 2n. Thus v≤2nv≤ 2n with probability at least 1−2/n1-2/n. On this event, for all sufficiently large n, (1−v/N)S⪯K⊤K⪯S,s⊤(K⊤K)−1s=v1−v/N≤Cn.(1-v/N)S K K S, s (K K)^-1s= v1-v/N≤ Cn. The equality is the rank-one inverse formula. This proves the second bound in equation 12 and allows spectral bounds for S to be transferred to K⊤K K with fixed constants. The first 4n4n queries satisfy ∑j≤4nyjyj⊤⪰cI _j≤ 4ny_jy_j cI except with probability e−cne^-cn. Since pj≥cn−αp_j≥ cn^-α on this block, S⪰cn−2αIS cn^-2αI. The upper trace bound follows from maxj⁡‖yj‖≤2 _j\|y_j\|≤ 2 and ∑jpj2≤C _jp_j^2≤ C. The first term p12‖y1‖2p_1^2\|y_1\|^2 gives a positive lower trace bound with probability tending to one. For q=⌊A∗L⌋=o⁡(n)q= A_*L =o(n), the n×qn× q matrix [y1,…,yq][y_1,…,y_q] has all singular values in [c,C][c,C] except with probability e−cne^-cn. The nonzero eigenvalues of its weighted Gram matrix therefore give λℓ(S)≥cpℓ2≥cℓ−2α _ (S)≥ cp_ ^2≥ c ^-2α for ℓ≤q ≤ q. This proves equation 11 and equation 13. For the inverse trace, put r=⌊n/4⌋r= n/4 and S>r=∑j>rpj2yjyj⊤S_>r= _j>rp_j^2y_jy_j . For a fixed unit vector u, nu⊤S>runu S_>ru is a weighted sum of independent squared standard Gaussians. Its weights obey ∑j>rpj2≤Cn1−2α,∑j>rpj4≤Cn1−4α,maxj>r⁡pj2≤Cn−2α. _j>rp_j^2≤ Cn^1-2α, _j>rp_j^4≤ Cn^1-4α, _j>rp_j^2≤ Cn^-2α. The Gaussian moment-generating-function bound, at a tail parameter of order n, followed by a 1/41/4-net of the sphere, yields ‖S>r‖op≤Cn−2α\|S_>r\|_op≤ Cn^-2α except with probability e−cne^-cn. The first r terms of S have rank at most r. Hence at least n−rn-r eigenvalues of K⊤K⪯SK K S are at most Cn−2αCn^-2α, and tr((K⊤K)−1/2)≥cnα+1.tr\! ((K K)^-1/2 )≥ cn^α+1. We next prove equation 14, including its uniformity. Write zi=nyiz_i= n\,y_i and an=n/Na_n= n/N. On equation 10, ‖N−1∑izizi⊤‖op≤C0\|N^-1 _iz_iz_i \|_op≤ C_0 for a fixed C0C_0. For a fixed unit vector v and 0<t≤10<t≤ 1, ℙ⁡(|v⊤zi|≤2t)≤C1t.P(|v z_i|≤ 2t)≤ C_1t. A Bernoulli Chernoff bound therefore gives ℙ1N#i:|v⊤zi|≤2t>C2t≤e−cNtP\! \ 1N\#\i:|v z_i|≤ 2t\>C_2t \≤ e^-cNt for a sufficiently large fixed C2C_2; when C2t≥1C_2t≥ 1 the event is empty. For each dyadic t∈[an,1]t∈[a_n,1], also including the smallest dyadic t not below ana_n, take a sphere net of mesh δ=c0t3/2δ=c_0t^3/2. Its cardinality is at most (3/δ)n(3/δ)^n. A union bound over these nets has failure probability at most Clog⁡(e/an)expCnlog⁡(e/an)−cNan=o⁡(1),C (e/a_n) \! \Cn (e/a_n)-cNa_n \=o(1), because N/n→∞N/n→∞ and Nan/n=N/n≫log⁡(eN/n)Na_n/n= N/n (eN/n). For an arbitrary unit vector u, choose a net point v within δ. If |u⊤zi|≤t|u z_i|≤ t but |v⊤zi|>2t|v z_i|>2t, then |(u−v)⊤zi|>t|(u-v) z_i|>t. Thus 1N#i:|(u−v)⊤zi|>t≤1Nt2∑i((u−v)⊤zi)2≤C0δ2t2≤Ct. 1N\#\i:|(u-v) z_i|>t\≤ 1Nt^2 _i((u-v) z_i)^2≤ C_0δ^2t^2≤ Ct. Combining these bounds and using monotonicity between dyadic thresholds and below ana_n proves, simultaneously for all unit u and all t>0t>0, 1N#i:|u⊤zi|≤t≤C⁡(t+an). 1N\#\i:|u z_i|≤ t\≤ C(t+a_n). For any nonzero B, diagonalize B⊤B‖B‖F2=∑kwkukuk⊤,wk≥0,∑kwk=1. B B\|B\|_F^2= _kw_ku_ku_k , w_k≥ 0, _kw_k=1. If ‖Bzi‖≤η‖B‖F\|Bz_i\|≤η\|B\|_F, at least half of the weight is on indices k for which |uk⊤zi|≤2η|u_k z_i|≤ 2η. Therefore 1N#i:‖Bzi‖≤η‖B‖F≤2∑kwk1N#i:|uk⊤zi|≤2η≤C⁡(η+an). split 1N\#\i:\|Bz_i\|≤η\|B\|_F\&≤ 2 _kw_k 1N\#\i:|u_k z_i|≤ 2η\\\ &≤ C(η+a_n). split This is equation 14; it is pathwise and allows B to depend on the data. Finally, for a fixed b>0b>0, a deterministic index i≥bNi≥ bN lies outside the first 4n4n and first ⌊A∗L⌋ A_*L queries once n is large. Deleting its weighted term leaves these blocks intact and only decreases the positive semidefinite tail matrix. The same centering argument applies with its weight set to zero; the projection still has rank n. All failure bounds are independent of the particular deterministic i≥bNi≥ bN. This proves the final assertion. ∎ Lemma C.2 (Bulk leverage on ℰYE_Y). On ℰYE_Y, tj:=(K⊤K)−1/2K⊤ejt_j:=(K K)^-1/2K e_j satisfies ‖tj‖2≤1\|t_j\|_2≤ 1 for every j. For every fixed b>0b>0, maxi≥bN⁡‖ti‖2≤Cbh. _i≥ bN\|t_i\|_2≤ C_bh. Proof. The matrix K(K⊤K)−1/2K(K K)^-1/2 has orthonormal columns, so every row has norm at most one. For i≥bNi≥ bN, ∥ti∥2≤pi∥(K⊤K)−1/2∥op∥yi∥2+1N∥(K⊤K)−1/2s∥2.\|t_i\|_2≤ p_i\|(K K)^-1/2\|_op\|y_i\|_2+ 1N\|(K K)^-1/2s\|_2. On ℰYE_Y, the first term is at most CbN−αnα=CbhC_bN^-αn^α=C_bh, while the second is Cn/N=o⁡(h)C n/N=o(h). Absorbing the latter term into the constant proves the result. ∎ On ℰYE_Y, define the response matrices for a fresh ξ∼N⁡(0,In)ξ N(0,I_n) by Sξ=ξ⊤K⊤Kξ,CY=ξξ⊤Sξ,HY=ξξ⊤Sξ,JY=HYCY−1/2.S_ξ=ξ K K\,ξ, C_Y=E ξ S_ξ, H_Y=E ξ S_ξ, J_Y=H_YC_Y^-1/2. Once Y is fixed, these are deterministic matrices. Lemma C.3 (Response moments on ℰYE_Y). On ℰYE_Y, cI⪯CY⪯CI,cI⪯JY⪯I.cI C_Y CI, cI J_Y I. The matrices CY,HY,JYC_Y,H_Y,J_Y commute with K⊤K K. For every fixed positive integer r, for all sufficiently large n, supY∈ℰY⁡[(ξ⊤K⊤Kξ)−r]≤Cr<∞. _Y _YE [(ξ K Kξ)^-r ]≤ C_r<∞. If (1−e)K⊤K⪯K~⊤K~⪯(1+e)K⊤K,0≤e≤12,(1-e)K K K K (1+e)K K, 0≤ e≤ 12, and J~ J is defined from K~⊤K~ K K in the same way, then ‖J~−JY‖op≤Ce.\| J-J_Y\|_op≤ Ce. The same conclusions apply to K−iK_-i on ℰY,−iE_Y,-i. Proof. For a fixed r, choose a fixed integer q>2rq>2r. In an eigenbasis of K⊤K K, equation 13 gives Sξ≥cq−2α∑ℓ=1qξℓ2.S_ξ≥ cq^-2α _ =1^q _ ^2. The negative r-th moment of χq2 _q^2 is finite; more explicitly it is 2−rΓ(q/2−r)/Γ⁡(q/2)2^-r (q/2-r)/ (q/2). This proves the inverse-moment bound, uniformly on ℰYE_Y. Sign symmetry in this eigenbasis makes CYC_Y and HYH_Y diagonal. Put cj=ξj2Sξ,hj=ξj2Sξ.c_j=E _j^2S_ξ, h_j=E _j^2 S_ξ. Cauchy–Schwarz gives hj2≤cjh_j^2≤ c_j. Also 1≤hj[ξj2Sξ],⁡[ξj2Sξ]=tr⁡(K⊤K)+2λj(K⊤K)≤C.1≤ h_j\,E[ _j^2 S_ξ], [ _j^2S_ξ]=tr(K K)+2 _j(K K)≤ C. Since ⁡[ξj2Sξ]≤⁡[ξj2Sξ]E[ _j^2 S_ξ]≤ E[ _j^2S_ξ], it follows that hj≥ch_j≥ c. A further application of Cauchy–Schwarz and the inverse second-moment bound gives cj≤(ξj4)1/2(Sξ−2)1/2≤Cc_j≤(E _j^4)^1/2(ES_ξ^-2)^1/2≤ C. Thus c≤cj≤Cc≤ c_j≤ C and c≤hj/cj≤1c≤ h_j/ c_j≤ 1, proving the spectral bounds and commutation statements. The relative Loewner comparison gives |S~ξ−1−Sξ−1|≤CeSξ−1| S_ξ^-1-S_ξ^-1|≤ CeS_ξ^-1 and |S~ξ−1/2−Sξ−1/2|≤CeSξ−1/2| S_ξ^-1/2-S_ξ^-1/2|≤ CeS_ξ^-1/2 pointwise. Testing the resulting matrix expectations against arbitrary vectors shows ‖C~−CY‖op+‖H~−HY‖op≤Ce\| C-C_Y\|_op+\| H-H_Y\|_op≤ Ce. Both covariance matrices are bounded below by a fixed positive multiple of the identity. The integral formula D−1/2=1π∫0∞t−1/2(D+tI)−1dtD^-1/2= 1π _0^∞t^-1/2(D+tI)^-1\,dt and the resolvent identity therefore give ∥C~−1/2−CY−1/2∥op≤Ce\| C^-1/2-C_Y^-1/2\|_op≤ Ce. Taking the difference of the products H~C~−1/2 H C^-1/2 and HYCY−1/2H_YC_Y^-1/2 proves the last bound. Only the Gram-matrix conditions were used, so the same proof applies on ℰY,−iE_Y,-i. ∎ Lemma C.4 (Population retrieval margin on ℰYE_Y). Fix b∈(0,1)b∈(0,1). There is cb>0c_b>0 such that for every fixed ε>0 >0, supi≥bNℙY(ti⊤JYyi<cbhormaxj≠i⁡|tj⊤JYyi|>εh|ℰY)→0. _i≥ bNP_Y\! (t_i J_Yy_i<c_bh\ \ or\ \ _j≠ i|t_j J_Yy_i|> h\, |\,E_Y )→ 0. Proof. We first record the relative square-root estimate needed below. For positive definite A0,A1A_0,A_1, suppose δ=∥A0−1/2(A1−A0)A0−1/2∥F≤12.δ=\|A_0^-1/2(A_1-A_0)A_0^-1/2\|_F≤ 12. Then ∥A1−1/2A01/2−I∥op≤Cδ.\|A_1^-1/2A_0^1/2-I\|_op≤ Cδ. Indeed, let A⁡(t)=A0+t⁡(A1−A0)A(t)=A_0+t(A_1-A_0) and H⁡(t)=A(t)1/2H(t)=A(t)^1/2. Differentiating H(t)2=A⁡(t)H(t)^2=A(t) gives HH˙+H˙H=A1−A0H H+ HH=A_1-A_0. If F⁡(t)=H(t)−1(A1−A0)H(t)−1F(t)=H(t)^-1(A_1-A_0)H(t)^-1, then in an eigenbasis of H⁡(t)H(t) with eigenvalues λj>0 _j>0, [H−1H˙]jk=λkλj+λkF(t)jk,‖H−1H˙‖F≤‖F⁡(t)‖F.[H^-1 H]_jk= _k _j+ _kF(t)_jk, \|H^-1 H\|_F≤\|F(t)\|_F. Since A⁡(t)⪰(1−tδ)A0A(t) (1-tδ)A_0, ∥F(t)∥F≤δ1−tδ,∥H(t)−1A01/2∥op≤(1−tδ)−1/2.\|F(t)\|_F≤ δ1-tδ, \|H(t)^-1A_0^1/2\|_op≤(1-tδ)^-1/2. It follows that ‖dt(H(t)−1A01/2)‖F≤δ(1−tδ)3/2. \| ddt (H(t)^-1A_0^1/2 ) \|_F≤ δ(1-tδ)^3/2. Integration from zero to one proves the claim. In particular, this estimate uses the Frobenius norm of the relative perturbation, not a condition-number-dependent additive perturbation bound. Fix a deterministic i≥bNi≥ bN, and work on ℰY∩ℰY,−iE_Y _Y,-i. Let J−iJ_-i denote the response operator obtained from K−i⊤K−iK_-i K_-i. Direct expansion gives K⊤K−K−i⊤K−i=pi2(1−1/N)yiyi⊤−piN(yis−i⊤+s−iyi⊤).K K-K_-i K_-i=p_i^2(1-1/N)y_iy_i - p_iN(y_is_-i +s_-iy_i ). Set ei:=‖(K−i⊤K−i)−1/2(K⊤K−K−i⊤K−i)(K−i⊤K−i)−1/2‖F.e_i:= \|(K_-i K_-i)^-1/2 (K K-K_-i K_-i )(K_-i K_-i)^-1/2 \|_F. Since pinα≍bhp_in^α _bh, ei≤Cb(h2+hnN)=Ob(h2).e_i≤ C_b (h^2+ h nN )=O_b(h^2). Hence (1−ei)K−i⊤K−i⪯K⊤K⪯(1+ei)K−i⊤K−i.(1-e_i)K_-i K_-i K K (1+e_i)K_-i K_-i. Lemma C.3 and the relative square-root estimate just proved (with A0=K−i⊤K−iA_0=K_-i K_-i and A1=K⊤KA_1=K K) give ‖JY−J−i‖op \|J_Y-J_-i\|_op ≤Cbh2, ≤ C_bh^2, (15) ‖(K⊤K)−1/2(K−i⊤K−i)1/2−I‖op \|(K K)^-1/2(K_-i K_-i)^1/2-I \|_op ≤Cbh2, ≤ C_bh^2, (16) ‖(K⊤K)−1/2JY−(K−i⊤K−i)−1/2J−i‖op \|(K K)^-1/2J_Y-(K_-i K_-i)^-1/2J_-i \|_op ≤Cbnαh2. ≤ C_bn^αh^2. (17) For j≠ij≠ i, define tj−i:=(K−i⊤K−i)−1/2(pjyj−s−iN),t_j^-i:=(K_-i K_-i)^-1/2 (p_jy_j- s_-iN ), and ti−i:=−(K−i⊤K−i)−1/2s−i/Nt_i^-i:=-(K_-i K_-i)^-1/2s_-i/N. These are rows of K−i(K−i⊤K−i)−1/2K_-i(K_-i K_-i)^-1/2, so ‖tj−i‖2≤1\|t_j^-i\|_2≤ 1. Also, tj=(K⊤K)−1/2(K−i⊤K−i)1/2tj−i−piN(K⊤K)−1/2yi,t_j=(K K)^-1/2(K_-i K_-i)^1/2t_j^-i- p_iN(K K)^-1/2y_i, hence maxj≠i⁡‖tj−tj−i‖2≤Cb(h2+h/N). _j≠ i\|t_j-t_j^-i\|_2≤ C_b(h^2+h/N). (18) Conditional on a realization of Y−iY_-i satisfying ℰY,−iE_Y,-i, the deleted objects are fixed and yi∼N⁡(0,In/n)y_i N(0,I_n/n) is independent. At this step we do not condition on ℰYE_Y, which also involves yiy_i. Each wrong deleted response is centered Gaussian with variance at most 1/n1/n. A Gaussian tail bound and a union bound give maxj≠i⁡|(tj−i)⊤J−iyi|=Oℙ(Ln)=oℙ(h). _j≠ i|(t_j^-i) J_-iy_i|=O_P\! ( Ln )=o_P(h). The replacement error from (tj−i,J−i)(t_j^-i,J_-i) to (tj,JY)(t_j,J_Y) is o⁡(h)o(h) by equation 15, equation 18, and ‖yi‖≤2\|y_i\|≤ 2. For the correct coordinate, put Ai:=(K−i⊤K−i)−1/2J−i.A_i:=(K_-i K_-i)^-1/2J_-i. This matrix is symmetric positive definite and independent of yiy_i, with tr⁡Ai≥cnα+1,‖Ai‖op≤Cnα,‖Ai‖F≤Cnα+1/2.trA_i≥ cn^α+1,\,\|A_i\|_op≤ Cn^α,\,\|A_i\|_F≤ Cn^α+1/2. Therefore ⁡[yi⊤Aiyi∣Y−i]=1ntr⁡Ai≥cnα,E[y_i A_iy_i Y_-i]= 1ntrA_i≥ cn^α, Var⁡(yi⊤Aiyi∣Y−i) (y_i A_iy_i Y_-i) =2n2‖Ai‖F2 = 2n^2\|A_i\|_F^2 ≤Cn2α−1. ≤ Cn^2α-1. Chebyshev gives yi⊤Aiyi≥cnα/2y_i A_iy_i≥ cn^α/2 with failure O⁡(1/n)O(1/n). Finally, ti⊤JYyi=piyi⊤(K⊤K)−1/2JYyi−1Ns⊤(K⊤K)−1/2JYyi.t_i J_Yy_i=p_iy_i (K K)^-1/2J_Yy_i- 1Ns (K K)^-1/2J_Yy_i. The matrix replacement error is Ob(h3)O_b(h^3) after multiplying by pip_i, while 1N|s⊤(K⊤K)−1/2JYyi|≤CnN=o(h). 1N |s (K K)^-1/2J_Yy_i |≤ C nN=o(h). Since pinα≍bhp_in^α _bh, the correct coordinate is at least cbh−o⁡(h)c_bh-o(h). The preceding bounds show that the margin failure event, intersected with ℰY∩ℰY,−iE_Y _Y,-i, has probability tending to zero, with a bound depending on i only through i≥bNi≥ bN. Since supi≥bNℙY((ℰY∩ℰY,−i)c)→0 _i≥ bNP_Y((E_Y _Y,-i)^c)→ 0 and ℙY(ℰY)→1P_Y(E_Y)→ 1, conditioning on ℰYE_Y yields the stated supremum. ∎ C.2 Condition on Y: randomness revealed through G Fix Y∈ℰY _Y. For the a-th row of X, define ξa=d(K⊤K)−1/2K⊤Xa,:⊤. _a= d\,(K K)^-1/2K X_a,: . Conditional on Y, the vectors ξ1,…,ξd _1,…, _d are independent N⁡(0,In)N(0,I_n), and Ga,:=σadξa⊤(K⊤K)1/2.G_a,:= _a d _a (K K)^1/2. The conditional law of the ξa _a does not depend on Y, so the family ξaa=1d\ _a\_a=1^d is independent of Y. Conversely, (Y,G)(Y,G) determines all ξa _a, because every σa _a is positive. Define the observed-gradient matrices ℋ =d(K⊤K)−1/2G⊤G(K⊤K)−1/2, =d(K K)^-1/2G G(K K)^-1/2, =d(K⊤K)−1/2G⊤Σ2G(K⊤K)−1/2. =d(K K)^-1/2G ^2G(K K)^-1/2. Conditional on Y, ℋ=∑a=1dσa2ξaξa⊤,=∑a=1dσa4ξaξa⊤.H= _a=1^d _a^2 _a _a , = _a=1^d _a^4 _a _a . Lemma C.5 (Observed-gradient estimates). There is C<∞C<∞ such that supY∈ℰYℙξ(‖ℋ−1ℋ‖op>Cn∣Y)→0, _Y _YP_ξ\! (\|HJ^-1H\|_op>Cn Y )→ 0, and supY∈ℰYξ[d∥G(K⊤K)−1/2∥∗∣Y]≤Cn3/2−β. _Y _YE_ξ\! [ d\,\|G(K K)^-1/2\|_* Y ]≤ Cn^3/2-β. For every fixed δ>0δ>0, supY∈ℰYℙξ(|1ntr⁡ℋ−∑a≥1σa2|>δor|1ntr⁡−∑a≥1σa4|>δ|Y)→0. _Y _YP_ξ\! ( | 1ntrH- _a≥ 1 _a^2 |>δ\ or\ | 1ntrJ- _a≥ 1 _a^4 |>δ\ |\ Y )→ 0. Proof. Split the first 2n2n rows from the rest: ℋ=ℋ0+ℋ1H=H_0+H_1 and =0+1J=J_0+J_1. Let Ξ0 _0 have rows ξa⊤ _a , a≤2na≤ 2n, and let B have rows σa2ξa⊤ _a^2 _a on the same block. Then 0=B⊤BJ_0=B B and ℋ0=Ξ0⊤B=B⊤Ξ0H_0= _0 B=B _0. Since B(B⊤B)−1B⊤B(B B)^-1B is an orthogonal projection, ℋ0−1ℋ0⪯ℋ00−1ℋ0⪯Ξ0⊤Ξ0⪯CnIH_0J^-1H_0 _0J_0^-1H_0 _0 _0 CnI except with probability e−cne^-cn. The same block gives ⪰cn−4βΞ0⊤Ξ0⪰cn1−4βIJ cn^-4β _0 _0 cn^1-4βI. On each dyadic block of rows with indices of order 2kn2^kn, the maximum weight σa2 _a^2 is at most C(2kn)−2βC(2^kn)^-2β and the unweighted Gaussian Gram matrix has operator norm at most C2knC2^kn, with failure e−c2kne^-c2^kn. Summing over these blocks, using 2β>12β>1, gives ‖ℋ1‖op≤C∑k≥1(2kn)1−2β≤Cn1−2β\|H_1\|_op≤ C _k≥ 1(2^kn)^1-2β≤ Cn^1-2β except with probability Ce−cnCe^-cn. Therefore ℋ1−1ℋ1⪯CnIH_1J^-1H_1 CnI. Finally, (ℋ0+ℋ1)−1(ℋ0+ℋ1)⪯2ℋ0−1ℋ0+2ℋ1−1ℋ1,(H_0+H_1)J^-1(H_0+H_1) 2H_0J^-1H_0+2H_1J^-1H_1, which proves the first assertion. The matrix dG(K⊤K)−1/2 d\,G(K K)^-1/2 has rows σaξa⊤ _a _a . Bound the nuclear norm of the first n rows by the sum of their row norms, and bound that of the remaining rows by n n times their Frobenius norm. Jensen’s inequality gives ξ[d∥G(K⊤K)−1/2∥∗∣Y]≤n∑a≤nσa+n(n∑a>nσa2)1/2≤Cn3/2−β. splitE_ξ\! [ d\,\|G(K K)^-1/2\|_* Y ]&≤ n _a≤ n _a+ n (n _a>n _a^2 )^1/2\\ &≤ Cn^3/2-β. split For the trace statements, ξtr⁡ℋn=s2,Varξ⁡(tr⁡ℋn)=2n∑a≤dσa4≤Cn.E_ξ trHn=s_2, _ξ\! ( trHn )= 2n _a≤ d _a^4≤ Cn. The analogous variance for tr⁡/ntrJ/n is 2∑a≤dσa8/n≤C/n2 _a≤ d _a^8/n≤ C/n. The corresponding finite sums converge to the stated infinite sums, so Chebyshev’s inequality completes the proof. All these laws and bounds are uniform in Y. ∎ For any update W∈ℝd×nW ^d× n, write mi(W)=K(K⊤K)−1G⊤Wyi,mji(W)=ej⊤mi(W).m_i(W)=K(K K)^-1G Wy_i, m_ji(W)=e_j m_i(W). When ‖ΣWyi‖2>0\| Wy_i\|_2>0, define a(i)(W)=dmi(W)‖ΣWyi‖2.a^(i)(W)= d\,m_i(W)\| Wy_i\|_2. If Wyi=0Wy_i=0, set a(i)(W)=0a^(i)(W)=0; all scores then vanish and strict retrieval fails. This convention covers zero denominators in the statements below. For the two baselines, nonzero queries have nonzero denominators almost surely. Corollary C.6 (Conditional mean-energy bound). For any update W with col⁡(W)⊆col⁡(G)col(W) (G), on the event ‖ℋ−1ℋ‖op≤Cn\|HJ^-1H\|_op≤ Cn, ‖a(i)(W)‖22≤Cnfor every i.\|a^(i)(W)\|_2^2≤ Cn every i. Proof. The case Wyi=0Wy_i=0 follows from the convention above. Otherwise, write Wyi=GqWy_i=Gq and q~=(K⊤K)1/2q q=(K K)^1/2q. Direct substitution gives a(i)(W)=K(K⊤K)−1/2ℋq~q~⊤q~.a^(i)(W)=K(K K)^-1/2 H q q J q. The matrix K(K⊤K)−1/2K(K K)^-1/2 has orthonormal columns. Hence ∥a(i)(W)∥22≤∥ℋ−1/2∥op2=∥ℋ−1ℋ∥op≤Cn.\|a^(i)(W)\|_2^2≤\|HJ^-1/2\|_op^2=\|HJ^-1H\|_op≤ Cn. ∎ Lemma C.7 (Correct conditional mean for gradient descent). Fix b∈(0,1)b∈(0,1) and u>0u>0. Then supY∈ℰYsupi≥bNℙξ(|ai(i)(WG)|>u∣Y)→0. _Y _Y _i≥ bNP_ξ\! (|a_i^(i)(W_G)|>u Y )→ 0. Proof. For WG=GW_G=G, set q=(K⊤K)1/2yiq=(K K)^1/2y_i and t=tit=t_i. The case q=0q=0 is immediate. Otherwise, ai(i)(WG)=t⊤ℋq⊤q.a_i^(i)(W_G)= t Hq q Jq. Write e=q/‖q‖e=q/\|q\| and t=⟨t,e⟩e+t⟂t= t,e e+t_ . The Gaussian coordinates along e and t⟂t_ are independent. Conditioning on those along e gives the distributional identity ai(i)(WG)=⟨t,e⟩∑a≤dσa2ga2∑a≤dσa4ga2+‖t⟂‖2g0,a_i^(i)(W_G) d= t,e _a≤ d _a^2g_a^2 _a≤ d _a^4g_a^2+\|t_ \|_2g_0, where g0,g1,g2,…g_0,g_1,g_2,… are independent standard Gaussians. Under a common coupling across d, the numerator is bounded by ∑a≥1σa2ga2<∞ _a≥ 1 _a^2g_a^2<∞ almost surely, and the denominator is at least σ12|g1|>0 _1^2|g_1|>0 almost surely. Thus the displayed ratio is uniformly tight, independently of Y and i. Lemma C.2 gives ‖ti‖≤Cbh\|t_i\|≤ C_bh for i≥bNi≥ bN, and h→0h→ 0. This proves the assertion. ∎ Lemma C.8 (Correct conditional mean for Muon). For every fixed u>0u>0 and δ>0δ>0, supY∈ℰYℙξ[1N#i:|ai(i)(WM)|>u>δ|Y]→0. _Y _YP_ξ\! [ 1N\#\i:|a_i^(i)(W_M)|>u\>δ\ |\ Y ]→ 0. Proof. Put W=WMW=W_M, and define νi=ti⊤(K⊤K)−1/2G⊤Wyi,vi=∥ΣWyi∥2. _i=t_i (K K)^-1/2G Wy_i, v_i=\| Wy_i\|_2. For vi>0v_i>0, ai(i)(W)=dνi/via_i^(i)(W)= d\, _i/v_i. The correct expansion is obtained with A=(K⊤K)−1G⊤W,z=A⊤s=W⊤G(K⊤K)−1s:νi=piyi⊤Ayi−1Nz⊤yi.A=(K K)^-1G W, z=A s=W G(K K)^-1s: _i=p_iy_i Ay_i- 1Nz y_i. On ℰYE_Y, pathwise even when A depends on the data, 1N∑i|yi⊤Ayi|≤tr⁡(|sym⁡A|YY⊤N)≤Cn‖A‖∗, 1N _i|y_i Ay_i| \! (|symA| Y N )≤ Cn\|A\|_*, where sym⁡A=(A+A⊤)/2symA=(A+A )/2 and ‖sym⁡A‖∗≤‖A‖∗\|symA\|_*≤\|A\|_*. Since ‖W‖op≤1\|W\|_op≤ 1, Lemma C.5 gives ξ[d‖A‖∗∣Y]≤∥(K⊤K)−1/2∥opξ[d∥G(K⊤K)−1/2∥∗∣Y]≤Cnα+3/2−β. splitE_ξ[ d\,\|A\|_* Y]&≤\|(K K)^-1/2\|_op\,E_ξ[ d\,\|G(K K)^-1/2\|_* Y]\\ &≤ Cn^α+3/2-β. split Also, although W depends on G, its operator norm bound implies ξ[d‖z‖22∣Y]≤ξ[d‖G(K⊤K)−1s‖22∣Y]=s2∥(K⊤K)−1/2s∥22≤Cn. splitE_ξ[d\|z\|_2^2 Y]& _ξ[d\|G(K K)^-1s\|_2^2 Y]\\ &=s_2\|(K K)^-1/2s\|_2^2≤ Cn. split Consequently d‖A‖∗=Oℙ(nα+3/2−β) d\,\|A\|_*=O_P(n^α+3/2-β) and d‖z‖2=Oℙ(n) d\,\|z\|_2=O_P( n), uniformly in Y. For fixed b>0b>0, use pi≤CbN−αp_i≤ C_bN^-α on i≥bNi≥ bN and 1N∑i|z⊤yi|≤(z⊤YY⊤Nz)1/2≤Cn‖z‖2 1N _i|z y_i|≤ (z Y Nz )^1/2≤ C n\|z\|_2 to obtain 1N∑i≥bN|dνi|≤CbN−αnd‖A‖∗+CNnd‖z‖2=Oℙ,b(hn1/2−β+N−1). split 1N _i≥ bN| d\, _i|&≤ C_b N^-αn d\,\|A\|_*+ CN n d\,\|z\|_2\\ &=O_P,b\! (hn^1/2-β+N^-1 ). split (19) Since WMW_M is the polar factor of the full-column-rank matrix G, ‖ΣWM‖F2=tr⁡[(G⊤G)−1G⊤Σ2G]=tr⁡(ℋ−1)≥tr⁡‖ℋ‖op≥tr⁡tr⁡ℋ. split\| W_M\|_F^2&=tr\! [(G G)^-1G ^2G ]=tr(H^-1J)\\ &≥ trJ\|H\|_op≥ trJtrH. split The trace limits in Lemma C.5 have positive finite values. Hence ‖ΣWM‖F≥c0>0\| W_M\|_F≥ c_0>0 with conditional probability tending to one, uniformly in Y. On this event, apply equation 14 pathwise with B=ΣWMB= W_M. For each fixed η∈(0,1)η∈(0,1), all but at most CN(η+n/N)CN(η+ n/N) indices satisfy vi≥c0η/n.v_i≥ c_0η/ n. No independence of B and Y is required. Counting the exceptional indices and then using Markov’s inequality for the remaining deterministic list of numerators gives 1N#i:|ai(i)(WM)|>u≤b+1N+Cη+CnN+nuc0ηN∑i≥bN|dνi|. split 1N\#\i:|a_i^(i)(W_M)|>u\≤ b+ 1N+Cη+C nN+ nuc_0η N _i≥ bN| d\, _i|. split By equation 19, the last term is Oℙ,b(hn1−β+n/Nuη)=oℙ(1).O_P,b\! ( hn^1-β+ n/Nuη )=o_P(1). For a given δ>0δ>0, choose b,ηb,η small and fixed so that b+Cη<δ/2b+Cη<δ/2, and then send n→∞n→∞. The exceptional observed-gradient event also has probability tending to zero uniformly in Y, which proves the assertion. ∎ For the chained update, the answer-coordinate scales cancel exactly: RN⁡(G)=RN⁡(XK)RN(G)=RN(XK), since every σa _a is positive. Define W~a,:=ξa⊤CY−1/2dξa⊤K⊤Kξa. W_a,:= _a C_Y^-1/2 d\, _a K K _a. Lemma C.9 (Fixed-Y approximation for the chain). There is C<∞C<∞ such that supY∈ℰYℙξ(‖WC−W~‖op>CLα+1nd|Y)→0, _Y _YP_ξ\! (\|W_C- W\|_op>CL^α+1 nd\ |\ Y )→ 0, and supY∈ℰYℙξ(maxa≤d∥(WC)a,:∥2>CLα+1nd|Y)→0. _Y _YP_ξ\! ( _a≤ d\|(W_C)_a,:\|_2>CL^α+1 nd\ |\ Y )→ 0. Moreover, supY∈ℰYsupi≤Nℙξ(‖mi(WC)−s1dK(K⊤K)−1/2JYyi‖∞>Cns2dLα+2|Y)→0. _Y _Y _i≤ NP_ξ\! ( \|m_i(W_C)- s_1dK(K K)^-1/2J_Yy_i \|_∞>C ns_2dL^α+2\ |\ Y )→ 0. Proof. We give the concentration, polar perturbation, and conditional-mean steps separately. Throughout the proof Y∈ℰY _Y is fixed. Concentration of the proxy Gram matrix. Write Za=CY−1/2ξaξa⊤CY−1/2ξa⊤K⊤Kξa,W~⊤W~=1d∑a=1dZa,ξZa=I.Z_a= C_Y^-1/2 _a _a C_Y^-1/2 _a K K _a, W W= 1d _a=1^dZ_a, _ξZ_a=I. Choose q=⌈AL⌉q= AL with A a sufficiently large fixed constant; in Lemma C.1, take A∗>2A_*>2A. In an eigenbasis of K⊤K K, ξa⊤K⊤Kξa≥cq−2α∑ℓ=1qξa,ℓ2. _a K K _a≥ cq^-2α _ =1^q _a, ^2. Let aT_a be the event ‖ξa‖2≤Cn\| _a\|^2≤ Cn and ∑ℓ≤qξa,ℓ2≥q/2 _ ≤ q _a, ^2≥ q/2. By choosing the fixed constants sufficiently large, Gaussian and chi-square tails give supY∈ℰYℙξ(ac∣Y)≤e−20L _Y _YP_ξ(T_a^c Y)≤ e^-20L for all sufficiently large n. On aT_a, ξa⊤K⊤Kξa≥cL1−2α,0⪯Za⪯CnL2α−1I. _a K K _a≥ cL^1-2α, 0 Z_a CnL^2α-1I. For Z^a=Zaa Z_a=Z_a1_T_a, positivity gives Z^a⪯IE Z_a I and Z^a2⪯CnL2α−1Z^a⪯CnL2α−1I.E Z_a^2 CnL^2α-1E Z_a CnL^2α-1I. Thus the centered truncated sum has variance parameter at most CdnL2α−1CdnL^2α-1 and summand norm at most CnL2α−1CnL^2α-1. The inverse moments in Lemma C.3 also give ‖Za‖op2≤C[‖ξa‖4(ξa⊤K⊤Kξa)2]≤Cn2,E\|Z_a\|_op^2≤ CE\! [ \| _a\|^4( _a K K _a)^2 ]≤ Cn^2, so Cauchy–Schwarz controls the truncation bias: ‖Z^a−I‖op≤Cne−10L.\|E Z_a-I\|_op≤ Cne^-10L. Apply matrix Bernstein to both signs of the centered truncated sum (tropp2012user, Theorem 6.1). Together with ℙ⁡(⋃a≤dac∣Y)≤de−20LP( _a≤ dT_a^c Y)≤ de^-20L, it gives ‖W~⊤W~−I‖op≤C⁡(Lαnd+L2αnd)≤CLαnd=:en\| W W-I\|_op≤ C (L^α nd+L^2α nd )≤ CL^α nd=:e_n with conditional probability tending to one uniformly in Y. The last inequality uses Lαn/d→0L^α n/d→ 0; the truncation bias is negligible on this scale. A relative polar estimate. For positive definite matrices M,PM,P with δ=‖P−I‖op≤1/2δ=\|P-I\|_op≤ 1/2, we claim ‖pol⁡(PM)−I‖op≤C⁡(1+log⁡cond⁡(M))δ,\|pol(PM)-I\|_op≤ C (1+ (M) )δ, (20) where cond⁡(M)=λmax(M)/λmin(M)cond(M)= _ (M)/ _ (M). Here is a direct proof of the spectral-norm estimate. Set E=P−IE=P-I and write the polar decomposition (I+tE)M=U⁡(t)H(t),0≤t≤1,(I+tE)M=U(t)H(t), 0≤ t≤ 1, where U⁡(t)U(t) is orthogonal and H⁡(t)H(t) is positive definite. The factors are differentiable along this nonsingular path. Put Ω⁡(t)=U(t)⊤U˙(t) (t)=U(t) U(t) and D⁡(t)=U(t)⊤E(I+tE)−1U(t).D(t)=U(t) E(I+tE)^-1U(t). The matrix D⁡(t)D(t) is symmetric, because E commutes with I+tEI+tE, and ‖D⁡(t)‖op≤2δ\|D(t)\|_op≤ 2δ. Differentiating the polar decomposition and subtracting its transpose gives the Sylvester equation ΩH+HΩ=DH−HD. H+H =DH-HD. Its integral solution is Ω=∫0∞e−sH(DH−HD)e−sHs. = _0^∞e^-sH(DH-HD)e^-sH\,ds. If μ=λmin(H)μ= _ (H) and Λ=λmax(H) = _ (H), then ‖Ω‖op≤2|D|∫0∞op⁡‖e−sH‖op‖He−sH‖ops≤C‖D‖op(1+log⁡(Λ/μ)).\| \|_op≤ 2\|D\|_op _0^∞\|e^-sH\|_op\|He^-sH\|_op\,ds≤ C\|D\|_op (1+ ( /μ) ). For the last inequality, split the integral at 1/Λ1/ and 1/μ1/μ: on the first interval use ‖He−sH‖op≤Λ\|He^-sH\|_op≤ ; on the middle interval use ‖He−sH‖op≤1/(es)\|He^-sH\|_op≤ 1/(es); on the last interval the product is at most μe−2sμ e^-2sμ. Moreover, Λμ=cond⁡((I+tE)M)≤1+δ1−δcond⁡(M)≤3cond⁡(M). μ=cond((I+tE)M)≤ 1+δ1-δcond(M)≤ 3cond(M). Since U⁡(0)=IU(0)=I and ‖U˙(t)‖op=‖Ω⁡(t)‖op\| U(t)\|_op=\| (t)\|_op, integration over t∈[0,1]t∈[0,1] proves equation 20. Application to WCW_C and the row bound. The factorization RN⁡(G)d=W~CY1/2(K⊤K)1/2 RN(G) d= W\,C_Y^1/2(K K)^1/2 is exact. On the proxy-Gram event, set P=(W~⊤W~)1/2,V=pol⁡(W~),M=CY1/2(K⊤K)1/2.P=( W W)^1/2, V=pol( W), M=C_Y^1/2(K K)^1/2. Then W~=VP W=VP and V⊤V=IV V=I. By Lemma C.3, the two factors in M commute, so M is symmetric positive definite. Its eigenvalues lie between cn−αcn^-α and C, and thus 1+log⁡cond⁡(M)≤CL1+ (M)≤ CL. Left multiplication by an isometry preserves the polar factorization, so WC=Vpol⁡(PM),‖P−I‖op≤Cen.W_C=Vpol(PM), \|P-I\|_op≤ Ce_n. Applying equation 20 yields ‖WC−W~‖op=‖pol⁡(PM)−P‖op≤‖pol⁡(PM)−I‖op+‖P−I‖op≤CLen≤CLα+1nd. split\|W_C- W\|_op&=\|pol(PM)-P\|_op\\ &≤\|pol(PM)-I\|_op+\|P-I\|_op\\ &≤ CLe_n≤ CL^α+1 nd. split For the row bound, all aT_a hold with conditional probability at least 1−de−20L1-de^-20L. On these events, maxa≤d∥W~a,:∥2≤CLα−1/2nd. _a≤ d\| W_a,:\|_2≤ CL^α-1/2 nd. The square matrix pol⁡(PM)pol(PM) is orthogonal and V=W~P−1V= WP^-1, with ‖P−1‖op≤2\|P^-1\|_op≤ 2 for large n. Therefore ∥(WC)a,:∥2=∥Va,:∥2≤2∥W~a,:∥2≤CLα−1/2nd,\|(W_C)_a,:\|_2=\|V_a,:\|_2≤ 2\| W_a,:\|_2≤ CL^α-1/2 nd, which is stronger than the stated row bound. Conditional-mean concentration. Fix a deterministic i and write b=CY−1/2yib=C_Y^-1/2y_i, so ‖b‖≤C\|b\|≤ C. The contribution from row a to coordinate j of mi(W~)m_i( W) is σadfaj,faj=(tj⊤ξa)(b⊤ξa)ξa⊤K⊤Kξa. _adf_aj, f_aj= (t_j _a)(b _a) _a K K _a. Its expectation is faj=tj⊤HYCY−1/2yi=tj⊤JYyiEf_aj=t_j H_YC_Y^-1/2y_i=t_j J_Yy_i. Gaussian moments, ‖tj‖≤1\|t_j\|≤ 1, and the inverse-moment bound imply faj2≤CEf_aj^2≤ C uniformly in a,j,i,Ya,j,i,Y. For each row, refine aT_a by also requiring maxj≤N⁡|tj⊤ξa|≤CL,|b⊤ξa|≤CL. _j≤ N|t_j _a|≤ C L, |b _a|≤ C L. Gaussian tails and a union bound over j allow the complement of this refined event to have probability at most e−20Le^-20L, after increasing fixed constants. On the refined event, |faj|≤CLα+1/2|f_aj|≤ CL^α+1/2. The truncated variables remain independent across a; their total weighted variance is at most Cs2/d2Cs_2/d^2, and their centered summand bound is CmaxaσaLα+1/2/dC _a _aL^α+1/2/d. Their total bias is at most Cs1e−10L/dCs_1e^-10L/d, which is negligible since s1≤ds2s_1≤ ds_2. Scalar Bernstein and a union bound over j≤Nj≤ N now give maxj≤N⁡|ej⊤mi(W~)−s1dtj⊤JYyi|≤Cd(s2L+maxa≤d⁡σaLα+3/2)≤Cs2dLα+3/2 split _j≤ N |e_j m_i( W)- s_1dt_j J_Yy_i |&≤ Cd ( s_2L+ _a≤ d _aL^α+3/2 )\\ &≤ C s_2dL^α+3/2 split with conditional failure probability tending to zero. The constants depend on i only through ‖yi‖≤2\|y_i\|≤ 2, so the same failure bound holds uniformly over deterministic i≤Ni≤ N. Finally, the replacement error obeys |ej⊤(mi(WC)−mi(W~))|≤∥G(K⊤K)−1/2tj∥2∥WC−W~∥op∥yi∥2.|e_j (m_i(W_C)-m_i( W))|≤\|G(K K)^-1/2t_j\|_2\,\|W_C- W\|_op\,\|y_i\|_2. Here ∥G(K⊤K)−1/2tj∥22=1d∑a≤dσa2(ξa⊤tj)2.\|G(K K)^-1/2t_j\|_2^2= 1d _a≤ d _a^2( _a t_j)^2. The weighted chi-square moment-generating-function bound, together with ‖tj‖≤1\|t_j\|≤ 1, ∑aσa4≤s22 _a _a^4≤ s_2^2, and a union bound over j, gives maxj≤N∥G(K⊤K)−1/2tj∥2≤Cs2Ld _j≤ N\|G(K K)^-1/2t_j\|_2≤ C s_2Ld with conditional probability tending to one uniformly in Y. Combining this with the operator-norm approximation and the proxy-mean bound yields ‖mi(WC)−s1dK(K⊤K)−1/2JYyi‖∞≤Cns2dLα+3/2≤Cns2dLα+2, \|m_i(W_C)- s_1dK(K K)^-1/2J_Yy_i \|_∞≤ C ns_2dL^α+3/2≤ C ns_2dL^α+2, which proves the last assertion. ∎ C.3 Condition on (Y,G)(Y,G): residual score randomness Once (Y,G)(Y,G) is fixed, any update W measurable with respect to (Y,G)(Y,G), its conditional mean, and all query-side quantities are fixed. The remaining randomness is the component of X not revealed by G=ΣXKG= XK. Define the retrieval score Sji(W)=xj⊤ΣWyi.S_ji(W)=x_j Wy_i. Lemma C.10 (Exact conditional score law). Let W be any measurable function of (Y,G)(Y,G). For a fixed query yiy_i, write i(W)=(S1i(W),…,SNi(W))⊤S_i(W)=(S_1i(W),…,S_Ni(W)) . Conditional on (Y,G)(Y,G), i(W)∼N⁡(mi(W),‖ΣWyi‖22d[IN−K(K⊤K)−1K⊤]),S_i(W) N\! (m_i(W), \| Wy_i\|_2^2d [I_N-K(K K)^-1K ] ), where mi(W)=K(K⊤K)−1G⊤Wyim_i(W)=K(K K)^-1G Wy_i. Proof. Conditional on Y, put U=K(K⊤K)−1/2U=K(K K)^-1/2. For each answer-coordinate row, Xa,:⊤=1dUξa+xa,⟂,xa,⟂=(IN−U⊤)Xa,:⊤.X_a,: = 1 dU _a+x_a, , x_a, =(I_N-U )X_a,: . These two Gaussian components are independent, and the residual rows are mutually independent. Given (Y,G)(Y,G), the first components are fixed because ξa=dσa−1(K⊤K)−1/2Ga,:⊤. _a= d\, _a^-1(K K)^-1/2G_a,: . The residual rows retain mean zero and covariance d−1(IN−UU⊤)d^-1(I_N-U ). Substituting their decomposition into i(W)=X⊤ΣWyiS_i(W)=X Wy_i gives the stated covariance. For the mean, the observed component satisfies 1d∑a=1dσaUξa(Wyi)a=K(K⊤K)−1G⊤Wyi, 1 d _a=1^d _aU _a(Wy_i)_a=K(K K)^-1G Wy_i, which is the claimed formula. ∎ Lemma C.11 (Wrong-score maximum). Suppose n→∞n→∞ and N/n→∞N/n→∞. Let P be any rank-n orthogonal projection in ℝNR^N, let a∈ℝNa ^N satisfy ‖a‖2≤Cn\|a\|_2≤ C n for a fixed C, and let z∼N⁡(0,IN−P)z N(0,I_N-P). Put τn=18log⁡(N/n) _n= 18 (N/n). There is a deterministic sequence εn→0 _n→ 0, depending only on n,N,Cn,N,C, such that for every P,a,iP,a,i, ℙ(ai+zi>maxj≠i(aj+zj))≤|ai|>τn+εn.P\! (a_i+z_i> _j≠ i(a_j+z_j) ) 1\|a_i|> _n\+ _n. Proof. Write z=(IN−P)gz=(I_N-P)g, where g∼N⁡(0,IN)g N(0,I_N), and set w=a−Pgw=a-Pg. Then a+z=g+wa+z=g+w and, since ‖Pg‖2∼χn2\|Pg\|^2 _n^2, ‖w‖2≤C′n\|w\|^2≤ C n except with probability e−cne^-cn. Let r=log⁡(N/n)r= (N/n) and BN=#⁡j:gj>rB_N=\#\j:g_j>r\. For all sufficiently large n, r≥1r≥ 1 and BN=Nℙ(g1>r)≥cNnr.EB_N=NP(g_1>r)≥ c Nnr. A Chernoff bound gives BN≥BN/2B_N _N/2 except with probability exp(−cNn/r) (-c Nn/r). On this event, if every coordinate of g+wg+w were below r/2r/2, then ‖w‖2≥r24BN≥crNn>C′n\|w\|^2≥ r^24B_N≥ cr Nn>C n for all sufficiently large n, a contradiction. Hence maxj⁡(aj+zj)≥r/2 _j(a_j+z_j)≥ r/2 except on the two events just bounded. If |ai|≤r/8|a_i|≤ r/8, then ai+zi≤r/4a_i+z_i≤ r/4 except with probability at most e−r2/128e^-r^2/128, because Var⁡(zi)≤1Var(z_i)≤ 1. Thus coordinate i cannot be the maximum on the intersection of these events. One may take, for sufficiently large n, εn=e−cn+exp(−cNn/r)+e−r2/128→0. _n=e^-cn+ (-c Nn/r)+e^-r^2/128→ 0. The argument does not assume independence of the coordinates of z, or independence between g and w. ∎ Theorem C.12 (Theorem 2 of the paper). Fix α>1α>1, 1/2<β<11/2<β<1, and a fixed deterministic positive sequence satisfying cσa−β≤σa≤Cσa−β(a≥1)c_σa^-β≤ _a≤ C_σa^-β\,(a≥ 1) for constants 0<cσ≤Cσ<∞0<c_σ≤ C_σ<∞. Fix 0<ϵ<(β−1/2)/20<ε<(β-1/2)/2 and Γ>1/(1−β) >1/(1-β). Let n→∞n→∞, and let N=N⁡(n)N=N(n) and d=d⁡(n)d=d(n) be integer sequences satisfying, eventually, n 1+(1−β+ϵ)/α≤N≤n 1+(1/2−ϵ)/α,n 1/(1−β)≤d≤nΓ.n^\,1+(1-β+ε)/α≤ N≤ n^\,1+(1/2-ε)/α,\,n^\,1/(1-β)≤ d≤ n . Then, on the same realized dataset, (Acc⁡(WC),Acc⁡(WM),Acc⁡(WG))→ℙ(1,0,0). (Acc(W_C),Acc(W_M),Acc(W_G) ) P(1,0,0). Proof. Let ℰXE_X denote ∥ℋ−1ℋ∥op≤Cn\\|HJ^-1H\|_op≤ Cn\. Lemma C.1 and Lemma C.5 give ℙ⁡(ℰYc)→0,supY∈ℰYℙξ(ℰXc∣Y)→0.P(E_Y^c)→ 0, _Y _YP_ξ(E_X^c Y)→ 0. Both WGW_G and WMW_M have their columns in col⁡(G)col(G). On ℰY∩ℰXE_Y _X, the exact score law and Corollary C.6 show that, for a nonzero query, the scores normalized by ‖ΣWyi‖/d\| Wy_i\|/ d have the law a(i)(W)+z,‖a(i)(W)‖2≤Cn,z∼N⁡(0,IN−K(K⊤K)−1K⊤).a^(i)(W)+z, \|a^(i)(W)\|_2≤ C n, z N\! (0,I_N-K(K K)^-1K ). The query is zero with probability zero; under the zero-denominator convention, such a query is not recovered in any case. The projection K(K⊤K)−1K⊤K(K K)^-1K has rank n. For GD, Lemma C.11 bounds conditional recovery probability by |ai(i)(WG)|>τn+εn1\|a_i^(i)(W_G)|> _n\+ _n on the good events. For any fixed u>0u>0, eventually τn≥u _n≥ u. Lemma C.7 therefore implies, for every fixed b∈(0,1)b∈(0,1), supi≥bNℙ⁡i is recovered by WG→0. _i≥ bNP\i is recovered by W_G\→ 0. Consequently ⁡[Acc⁡(WG)]≤b+1N+supi≥bNℙ⁡i is recovered by WG.E[Acc(W_G)]≤ b+ 1N+ _i≥ bNP\i is recovered by W_G\. Send n→∞n→∞ and then b↓0b 0. Markov’s inequality gives Acc⁡(WG)→ℙ0Acc(W_G) P0. For Muon, for any fixed u>0u>0 and eventually τn≥u _n≥ u, 1N#i:|ai(i)(WM)|>τn≤1N#i:|ai(i)(WM)|>u→ℙ0 1N\#\i:|a_i^(i)(W_M)|> _n\≤ 1N\#\i:|a_i^(i)(W_M)|>u\ P0 by Lemma C.8 and ℙ⁡(ℰYc)→0P(E_Y^c)→ 0. This fraction lies in [0,1][0,1], so its expectation also tends to zero. Averaging Lemma C.11 over i on ℰY∩ℰXE_Y _X, and adding the probabilities of the complementary events, yields ⁡[Acc⁡(WM)]→0,Acc⁡(WM)→ℙ0.E[Acc(W_M)]→ 0, (W_M) P0. For the chain, fix b∈(0,1)b∈(0,1) and a deterministic i≥bNi≥ bN. Lemma C.4 implies that, conditional on ℰYE_Y, with probability tending to one uniformly over such i, ti⊤JYyi≥cbh,maxj≠i⁡|tj⊤JYyi|≤cb4h.t_i J_Yy_i≥ c_bh, _j≠ i|t_j J_Yy_i|≤ c_b4h. Fix a query realization satisfying this margin. Lemma C.9 then implies mii(WC)−maxj≠i⁡mji(WC)≥cb2s1dhm_i(W_C)- _j≠ im_ji(W_C)≥ c_b2 s_1dh except on an observed-gradient event whose conditional probability tends to zero uniformly in Y and deterministic i. Indeed, the ratio of its mean-approximation error to (s1/d)h(s_1/d)h is ns2Lα+2s1h≤Cn−ϵLα+2→0. ns_2\,L^α+2s_1h≤ Cn^-εL^α+2→ 0. Intersect this event with the row-norm event in the same lemma, whose conditional failure probability also tends to zero uniformly in Y. On that intersection, ‖ΣWCyi‖2d≤‖yi‖2d(∑a≤dσa2∥(WC)a,:∥22)1/2≤Cns2dLα+1. \| W_Cy_i\|_2 d≤ \|y_i\|_2 d ( _a≤ d _a^2\|(W_C)_a,:\|_2^2 )^1/2≤ C ns_2dL^α+1. Now condition on (Y,G)(Y,G). Lemma C.10 shows that every centered score is Gaussian with variance bounded by the square of the last display. A union bound over j≤Nj≤ N gives ℙ(maxj≤N|Sji(WC)−mji(WC)|>Cns2dLα+3/2|Y,G)→0P\! ( _j≤ N|S_ji(W_C)-m_ji(W_C)|>C ns_2dL^α+3/2\ |\ Y,G )→ 0 uniformly over the stated good events. The fluctuation threshold is o⁡((s1/d)h)o((s_1/d)h) by the same scaling calculation. Therefore the correct score is strictly greater than all wrong scores with probability tending to one. Combining the query, observed-gradient, and residual-score bounds yields supi≥bNℙ⁡i is not recovered by WC→0. _i≥ bNP\i is not recovered by W_C\→ 0. This is a bound uniform over deterministic i, not a claim that a single union bound recovers every query. It suffices because ⁡[1−Acc⁡(WC)]≤b+1N+supi≥bNℙ⁡i is not recovered by WC.E[1-Acc(W_C)]≤ b+ 1N+ _i≥ bNP\i is not recovered by W_C\. Send n→∞n→∞ and then b↓0b 0, and use Markov’s inequality to obtain Acc⁡(WC)→ℙ1Acc(W_C) P1. All three updates are defined on the same realized dataset. A union bound over their three accuracy deviation events proves (Acc⁡(WC),Acc⁡(WM),Acc⁡(WG))→ℙ(1,0,0). (Acc(W_C),Acc(W_M),Acc(W_G) ) P(1,0,0). Equivalently, for every fixed δ∈(0,1)δ∈(0,1), ℙAcc(WC)≥1−δ,maxAcc(WM),Acc(WG)≤δ→1.P\! \Acc(W_C)≥ 1-δ, \Acc(W_M),Acc(W_G)\≤δ \→ 1. ∎ Appendix D Appendix: Proofs in section 4 Proof of Proposition 5 of the paper. Write =∘k∘NSslice∘k,=ANlayer∘ANcol∘ANrow.T=A k _slice k, =AN_layer _col _row. We first verify the chained-LMO representation. For an axis a∈layer,row,cola∈\layer,row,col\, let ℱaF_a denote the fibers along that axis and define a=Z:‖Zf‖2≤1for every f∈ℱa.K_a= \Z:\|Z_f\|_2≤ 1\ for every f _a \. Since the Frobenius inner product decomposes over fibers, maxZ∈a⁡⟨M,Z⟩=∑f∈ℱa‖Mf‖2, _Z _a M,Z = _f _a\|M_f\|_2, and a maximizer is obtained by normalizing every nonzero fiber independently. Hence ANaAN_a is, up to the conventional sign, an exact LMO over aK_a. Similarly, with slice=Z:∥Zℓ,:,:∥op≤1for every slice ℓ,K_slice= \Z:\|Z_ ,:,:\|_op≤ 1\ for every slice \, we have maxZ∈slice⟨M,Z⟩=∑ℓ∥Mℓ,:,:∥∗, _Z _slice M,Z = _ \|M_ ,:,:\|_*, whose maximizer is obtained by applying the matrix sign independently to each slice. Thus every constituent of T is an exact LMO, and T is a chained LMO. It remains to show that the composition is not itself an LMO. For G∈ℝm×nG ^m× n, let ι(G)1,:,:=G,ι(G)ℓ,:,:=0(ℓ>1). (G)_1,:,:=G, (G)_ ,:,:=0 ( >1). Row- and column-axis normalization preserve the signs and zero pattern of the nonzero slice. The subsequent layer-axis normalization acts on fibers of the form (c,0,…,0)(c,0,…,0), and therefore ⁡(ι⁡(G))=ι⁡(Sgn⁡(G)).A( (G))= (Sgn(G)). Consequently, for every k≥1k≥ 1, ∘k(ι⁡(G))=ι⁡(Sgn⁡(G)),A k( (G))= (Sgn(G)), and hence ⁡(ι⁡(G))=ι⁡(Sgn∘msign∘Sgn⁡(G)).T( (G))= \! (Sgn (G) ). By Lemma B.8, for every m,n≥5m,n≥ 5 there exists G∈ℝm×nG ^m× n such that ⟨G,Sgn∘msign∘Sgn⁡(G)⟩<0. G,Sgn (G) <0. Hence ⟨ι⁡(G),⁡(ι⁡(G))⟩<0. (G),T( (G)) <0. Thus TensorChain admits a negative-alignment input, and by Lemma 2, it is not a LMO of any fixed set. A smooth counterexample where TensorChain optimization diverges can be found by using Proposition 1. ∎ Appendix E Appendix: Replacing NS with Cheaper Chains (a) FFN replacement ablation. (b) Layerwise ablation. Figure 2: Replacing Newton–Schulz with cheaper axis-wise normalizations on Qwen3-0.6B after 8.8B training tokens. (a) NS5 is retained on Q,K,V,OQ,K,V,O while the FFN uses the indicated replacement. (b) A single row–column equilibration sweep is applied to the indicated parameter classes, with NS5 retained elsewhere. Experimental setup. We study whether inexpensive axis-wise normalizations can replace the Newton–Schulz (NS) step of Muon on the olmo2norm Qwen3-0.6B architecture. We use QK normalization, unfused QKV projections, Kaiming initialization, tied input/output embeddings, variable-length attention, selective activation checkpointing, and compilation. Training and validation use fineweb_edu_350bt and fineweb_edu_350bt_eval, respectively. We train for 8.8B tokens (1×1× Chinchilla) with sequence length 4096, global batch size 256 (1.05M tokens), and 8401 optimization steps. All runs use one node with 8 GPUs under replicated data parallelism and seed 42; validation is performed every 200 steps on one batch. Checkpoints are retained and training metrics are logged to TensorBoard. We use 10% learning-rate warmup followed by linear decay to zero, matrix weight decay 0.10.1, momentum coefficient β=0.95β=0.95, and AdamW learning rate 8×10−38× 10^-3 for non-matrix parameters. All matrix updates use the same s_mu_scale=max scaling, 0.2max⁡(fan_in,fan_out)0.2 (fan\_in,fan\_out). Muon uses five Newton–Schulz iterations (NS5). For each variant, we sweep the matrix learning rate over 3,5,7.5,10×10−3\3,5,7.5,10\× 10^-3 and report the best validation loss. For short-axis normalization, whose best point lies at the lower boundary of the sweep, we additionally evaluate 1×10−31× 10^-3 and 1.5×10−31.5× 10^-3. Replacing NS on the FFN. Our primary motivation is computational: the large FFN matrices account for the largest share of the NS computation, making them a natural target for a cheaper replacement. We therefore retain NS5 on the attention projections Q,K,V,OQ,K,V,O and replace NS on the FFN matrices W1,W2,W3W_1,W_2,W_3. We compare row and column ℓ2 _2 normalization, normalization along the long or short matrix axis, and alternating row–column equilibration. For the latter, we test both a single sweep (k=1k=1) and five repeated sweeps (k=5k=5). We additionally evaluate row normalization and repeated row–column equilibration as replacements for NS on all matrix parameters. Among the FFN replacements, row–column equilibration performs best, while increasing from one to five sweeps yields only a small improvement. We therefore fix k=1k=1 for the subsequent weight-wise ablation. Which weights can use the cheaper replacement? Having identified a single row–column sweep as the strongest inexpensive replacement, we next vary which parameter classes use it while retaining NS5 on all remaining matrices. We consider the FFN alone, Q,KQ,K alone, V,OV,O alone, and combinations of these groups. The strongest configuration replaces NS only on Q,KQ,K, reaching validation loss 2.59452.5945 compared with 2.59052.5905 for Muon. Replacing additional parameter classes consistently degrades performance, and no configuration outperforms Muon. Table 5: Selected Newton–Schulz replacement results on Qwen3-0.6B after 8.8B training tokens. Each entry reports the best result over the matrix learning-rate sweep. Δ denotes the difference from Muon. Variant Val. loss Δ Muon (NS5 on all matrices) 2.5905 – NS5 on QKVO, replacement on FFN Row–column, k=5k=5 2.6022 +0.012 Row–column, k=1k=1 2.6059 +0.015 Long-axis ℓ2 _2 2.6134 +0.023 Column ℓ2 _2 2.6136 +0.023 Row ℓ2 _2 2.6255 +0.035 Short-axis ℓ2 _2 2.6509 +0.060 Row–column (k=1k=1) on listed weights, NS5 elsewhere Q,KQ,K 2.5945 +0.004 Q,KQ,K, FFN 2.6026 +0.012 FFN 2.6059 +0.015 V,OV,O 2.6209 +0.030 V,OV,O, FFN 2.6420 +0.052 All matrices (row–column, k=5k=5) 2.6463 +0.056 Overall, cheap axis-wise normalization can approach NS when used selectively, but none of the tested replacements matches Muon. The small difference between one and five row–column sweeps also suggests that simply iterating equilibration further does not close the gap to spectral normalization. Appendix F Appendix: Raw Experiment Results This section provides the complete experimental record underlying the comparisons in the main text. We report the full Qwen3-0.6B baseline sweeps under WSD and cosine schedules, the Muon learning-rate sweeps used for tuning, the full train/validation trajectories and endpoint losses, the TensorChain grouping ablation, and the cosine token-efficiency analysis. Unless otherwise stated, the experimental configuration follows Section 5.1. We retain the cosine results in addition to the main WSD experiments to document schedule-dependent differences in optimizer ranking and learning-rate sensitivity. F.1 Full Qwen3-0.6B baseline comparison The two tables below report the complete common learning-rate grids for the Qwen3-0.6B, 1M-token-batch, 1×1× Chinchilla baseline comparison. TensorChain grouping variants are reported separately in Section F.5. For AdaMuon, the learning rate grid was insufficient so we ran two additional experiments, lr=27e−3lr=27e-3 and 40e−340e-3 and could find that the optimal learning rate is 27e−327e-3 with val loss 2.8405, which is strictly worse than TensorChain’s optimal val loss 2.8379. Table 6: Full Qwen3-0.6B, 1M-token-batch, 1×1× Chinchilla baseline sweep under WSD. Entries are final validation losses; methods are ordered by their best validation loss, and bold marks the best learning rate for each method. Method 1×10−31× 10^-3 2×10−32× 10^-3 3×10−33× 10^-3 5×10−35× 10^-3 8×10−38× 10^-3 TensorChain 2.9089 2.9002 2.9000 2.9040 2.9131 Muon2 2.9171 2.9106 2.9104 2.9157 2.9260 Aurora-LN 2.9204 2.9128 2.9116 2.9130 2.9225 AdaMuon 2.9290 2.9156 2.9142 2.9188 2.9294 MuonEq 2.9235 2.9163 2.9151 2.9186 2.9278 MuonEq-diag 2.9258 2.9184 2.9167 2.9211 2.9305 Muon+ 2.9279 2.9203 2.9180 2.9234 2.9334 Aurora 2.9280 2.9196 2.9184 2.9238 2.9335 Muon 2.9299 2.9224 2.9201 2.9251 2.9354 NorMuon 2.9295 2.9214 2.9201 2.9256 2.9357 Table 7: Full Qwen3-0.6B, 1M-token-batch, 1×1× Chinchilla baseline sweep under cosine decay. Entries are final validation losses; methods are ordered by their best validation loss, and bold marks the best learning rate for each method. Method 5×10−35× 10^-3 8×10−38× 10^-3 12×10−312× 10^-3 16×10−316× 10^-3 20×10−320× 10^-3 TensorChain 2.8501 2.8409 2.8379 2.8383 2.8399 AdaMuon 2.8674 2.8609 2.8568 2.8519 2.8458 Aurora-LN 2.8614 2.8520 2.8474 2.8476 2.8482 Muon2 2.8626 2.8543 2.8497 2.8498 2.8515 MuonEq 2.8675 2.8578 2.8540 2.8548 2.8563 MuonEq-diag 2.8708 2.8624 2.8572 2.8562 2.8567 Aurora 2.8739 2.8658 2.8624 2.8627 2.8646 Muon+ 2.8739 2.8663 2.8631 2.8630 2.8684 Muon 2.8756 2.8680 2.8639 2.8645 2.8659 NorMuon 2.8767 2.8685 2.8650 2.8652 2.8657 (a) WSD learning rate tuning loss. (b) COS learning rate tuning Figure 3: Full Qwen3-0.6B baseline comparison. F.2 Muon learning-rate sweeps We report each tuning setup separately for Muon and use it for all baselines. The results are below. Table 8: Muon learning-rate sweep for Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla under WSD. Matrix LR 1×10−31× 10^-3 2×10−32× 10^-3 3×10−33× 10^-3 5×10−35× 10^-3 8×10−38× 10^-3 Final val. loss 2.9299 2.9224 2.9201 2.9251 2.9354 Table 9: Muon learning-rate sweep for Qwen3-1.7B, 4M-token batch, 1×1× Chinchilla under WSD. Matrix LR 0.5×10−30.5× 10^-3 1×10−31× 10^-3 2×10−32× 10^-3 4×10−34× 10^-3 8×10−38× 10^-3 Final val. loss 2.6878 2.6672 2.6523 2.6493 2.6589 Table 10: Muon learning-rate sweep for Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla under cosine decay. Matrix LR 1×10−31× 10^-3 2×10−32× 10^-3 3×10−33× 10^-3 5×10−35× 10^-3 8×10−38× 10^-3 12×10−312× 10^-3 16×10−316× 10^-3 20×10−320× 10^-3 Final val. loss 2.9201 2.8993 2.8857 2.8756 2.8680 2.8639 2.8645 2.8659 Table 11: Muon learning-rate sweep for Qwen3-0.6B, 4M-token batch, 1×1× Chinchilla under cosine decay. Matrix LR 8×10−38× 10^-3 12×10−312× 10^-3 18×10−318× 10^-3 27×10−327× 10^-3 40×10−340× 10^-3 60×10−360× 10^-3 Final val. loss 2.9650 2.9568 2.9516 2.9462 2.9459 2.9468 Table 12: Muon learning-rate sweep for Qwen3-1.7B, 1M-token batch, 1×1× Chinchilla under cosine decay. Matrix LR 4×10−34× 10^-3 6×10−36× 10^-3 9×10−39× 10^-3 13.5×10−313.5× 10^-3 20×10−320× 10^-3 Final val. loss 2.6019 2.5974 2.5987 2.6042 2.6139 Table 13: Muon learning-rate sweep for Qwen3-1.7B, 4M-token batch, 1×1× Chinchilla under cosine decay. Matrix LR 5×10−35× 10^-3 8×10−38× 10^-3 10×10−310× 10^-3 15×10−315× 10^-3 20×10−320× 10^-3 Final val. loss 2.6228 2.6120 2.6091 2.6047 2.6049 F.3 WSD training and validation trajectories The following two tables aggregate the final losses across WSD setups. Since all methods within a setup use the same selected learning rate, the learning rate is omitted here. Table 14: Final training loss for the WSD trajectory experiments. Bold marks the best method in each row. Setup Muon Aurora MuonEq Muon2 Aurora-LN TensorChain 0.6B, 1M, 1×1× 2.9446 2.9437 2.9397 2.9350 2.9362 2.9249 1.7B, 4M, 4×4× 2.5254 2.5244 2.5226 2.5208 2.5188 2.5181 Table 15: Final validation loss for the WSD trajectory experiments. Bold marks the best method in each row. Setup Muon Aurora MuonEq Muon2 Aurora-LN TensorChain 0.6B, 1M, 1×1× 2.9201 2.9184 2.9151 2.9104 2.9116 2.8990 1.7B, 4M, 4×4× 2.5310 2.5301 2.5279 2.5264 2.5242 2.5236 (a) Training loss. (b) Validation loss. Figure 4: WSD trajectories for Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla. (a) Training loss. (b) Validation loss. Figure 5: WSD trajectories for Qwen3-1.7B, 4M-token batch, 4×4× Chinchilla. F.4 Cosine training and validation trajectories We analogously aggregate the four cosine settings. Again, learning rates are omitted because all methods within a given setup use the same selected learning rate. Table 16: Final training loss for the cosine trajectory experiments. Bold marks the best method in each row. Setup Muon Aurora MuonEq Muon2 Aurora-LN TensorChain 0.6B, 1M, 1×1× 2.8899 2.8889 2.8806 2.8771 2.8750 2.8636 0.6B, 4M, 4×4× 2.7934 2.7897 2.7845 2.7851 2.7804 2.7748 1.7B, 1M, 1×1× 2.6036 2.6017 2.5978 2.5941 2.5942 2.5922 1.7B, 4M, 4×4× 2.3859 2.3850 2.3814 2.3780 2.3751 2.3778 Table 17: Final validation loss for the cosine trajectory experiments. Bold marks the best method in each row. Setup Muon Aurora MuonEq Muon2 Aurora-LN TensorChain 0.6B, 1M, 1×1× 2.8639 2.8624 2.8540 2.8497 2.8474 2.8379 0.6B, 4M, 4×4× 2.7527 2.7487 2.7434 2.7443 2.7391 2.7344 1.7B, 1M, 1×1× 2.5974 2.5964 2.5924 2.5885 2.5882 2.5857 1.7B, 4M, 4×4× 2.4823 2.4810 2.4784 2.4744 2.4734 2.4747 (a) Training loss. (b) Validation loss. Figure 6: Cosine trajectories for Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla. (a) Training loss. (b) Validation loss. Figure 7: Cosine trajectories for Qwen3-0.6B, 4M-token batch, 4×4× Chinchilla. (a) Training loss. (b) Validation loss. Figure 8: Cosine trajectories for Qwen3-1.7B, 1M-token batch, 1×1× Chinchilla. (a) Training loss. (b) Validation loss. Figure 9: Cosine trajectories for Qwen3-1.7B, 4M-token batch, 4×4× Chinchilla. (a) Training loss. (b) Validation loss. Figure 10: Cosine trajectories for Qwen3-1.7B, 4M-token batch, 4×4× Chinchilla, distinct tokens. F.5 TensorChain grouping ablation TensorChain requires choosing which compatible matrices are grouped along the layer axis. For a 28-layer model, consecutive forms four contiguous groups, e.g., 1:7,8:14,15:21,22:28\1:7\,\8:14\,\15:21\,\22:28\, whereas stride interleaves depth, e.g., 1,5,9,…,25,2,6,10,…,26,3,7,11,…,27,4,8,12,…,28\1,5,9,…,25\,\2,6,10,…,26\,\3,7,11,…,27\,\4,8,12,…,28\. Greedy assigns compatible matrices to parameter-balanced groups, while all jointly groups all compatible matrices 1,2,⋯28\1,2,·s 28\. The results are extracted from the same Qwen3-0.6B full sweeps above. Table 18: TensorChain grouping ablation on Qwen3-0.6B, 1M-token batch, 1×1× Chinchilla. Each entry reports the best point from the corresponding full sweep. WSD Cosine Grouping Best LR Val. loss Best LR Val. loss Consecutive 3×10−33× 10^-3 2.9070 12×10−312× 10^-3 2.8474 Stride 3×10−33× 10^-3 2.8990 12×10−312× 10^-3 2.8357 Greedy 3×10−33× 10^-3 2.8994 16×10−316× 10^-3 2.8370 All 3×10−33× 10^-3 2.9000 12×10−312× 10^-3 2.8379 Consecutive grouping performs worse near its optimum than stride, greedy, and all-layer grouping. Since the latter three mix matrices across more distant depths despite using different assignment rules, this suggests that cross-depth mixing matters more than the precise partition. F.6 Cosine token efficiency The cosine schedule yields smooth validation trajectories throughout training, so we additionally compare methods by the number of training tokens required to reach the same validation loss. We do not use this metric for WSD because its terminal decay concentrates much of the loss reduction near the end of training, making within-run crossing times less informative. For TensorChain and a baseline B, define the common target ℓ⋆(TC,B)=max⁡ℓTC(Dmax),ℓB(Dmax), (TC,B)= \! \ _TC(D_ ), _B(D_ ) \, the worse of the two final validation losses. Both methods therefore reach the target within the observed training horizon, so no extrapolation is required. For a method M, let DM(ℓ)=infD:ℓ^M(D)≤ℓ,D_M( )= \! \D: _M(D)≤ \, where ℓ^M(D) _M(D) is obtained by linear interpolation between the two adjacent validation checkpoints bracketing the target. We define TE⁡(TC,B)=100(1−DTC(ℓ⋆)DB(ℓ⋆))%.TE(TC;B)=100 (1- D_TC( )D_B( ) )\%. Positive values indicate that TensorChain reaches the same validation loss using fewer training tokens, while negative values favor the baseline. Since each pairwise comparison fixes the batch size, the same percentage is equivalently the reduction in optimization steps. Table 19: Cosine-schedule token-efficiency gain of TensorChain. Each entry is the percentage reduction in training tokens required to reach the same validation-loss target. Positive values favor TensorChain. Setting Aurora-LN Muon2 MuonEq Aurora Muon 0.6B, 1M, 1×1× 8.2 9.0 10.3 12.4 12.8 0.6B, 4M, 4×4× 5.4 7.6 7.4 8.9 10.0 1.7B, 1M, 1×1× 3.9 4.0 6.2 7.5 7.8 1.7B, 4M, 4×4× -2.0 -0.5 6.7 7.4 7.8 F.7 Comparison with the Teon optimizer Since we could not find an open-sourced implementation of Teon (zhang2026teon), we implemented our own version of Teon. We concatenated the two consecutive Q, K, V weights and applied Newton–Schulz, following the suggestion of (zhang2026teon). We tuned the learning rate of Teon in Qwen3 0.6B pretraining setup, with 1M batch size and 1x chinchilla tokens, with WSD schedule with 10%\% warmup and decay. To our surprise, the Teon optimizer performed worse than Muon for our setup; The best validation loss for Teon was 2.9222, whereas Muon was 2.9201. Hence we decided to not include Teon in our baselines.

Instant research discovery