Muon is Not That Special: Random or Inverted Spectra Work Just as Well

# Random or Inverted Spectra Work Just as Well
Source: [https://arxiv.org/html/2605.11181](https://arxiv.org/html/2605.11181)
## Muon is Not That Special: Random or Inverted Spectra Work Just as Well

Zakhar Shumaylov1Nathaël Da Costa2Peter Zaika1Bálint Mucsányi2Alex Massucco1 Yoav Gelberg3Carola\-Bibiane Schönlieb1Yarin Gal3Philipp Hennig2

1University of Cambridge2University of Tübingen3University of Oxford

###### Abstract

The recent empirical success of theMuonoptimizer has renewed interest in non\-Euclidean optimization, typically justified by similarities with second\-order methods, and linear minimization oracle \(LMO\) theory\. In this paper, we challenge this geometric narrative through three contributions, demonstrating that precise geometric structure is not the key factor affecting optimization performance\. First, we introduceFreon, a family of optimizers based on Schatten \(quasi\-\)norms, powered by a novel, provably optimal QDWH\-based iterative approximation\.Freonnaturally interpolates betweenSGDandMuon, while smoothly extrapolating into the quasi\-norm regime\. Empirically, the best\-performing Schatten parameters for GPT\-2 lie strictly within the quasi\-norm regime, and thuscannotbe represented by any unitarily invariant LMO\. Second, noting thatFreonperforms well across a wide range of exponents, we introduceKaon, an absurd optimizer that replaces singular values with random noise\. Despite lacking any coherent geometric structure,KaonmatchesMuon’s performance and retains classical convergence guarantees, proving that strict adherence to a precise geometry is practically irrelevant\. Third, having shown that geometry is not the primary driver of performance, we demonstrate it is instead controlled by two local quantities:alignmentanddescent potential\. Ultimately, each optimizer must tune its step size around these two quantities\. While their dynamics are difficult to predict a priori, evaluating them within a stochastic random feature model yields a precise insight:Muonsucceeds not by tracking an ideal global geometry, but by guaranteeing step\-size optimality\.

## 1Introduction

First\-order optimization algorithms have become increasingly central to modern machine learning, fueled by their routine use in training models with billions of parameters on trillions of tokens\. While adaptive gradient methods likeAdamW\(Kingma and Ba,[2017](https://arxiv.org/html/2605.11181#bib.bib23); Duchiet al\.,[2011](https://arxiv.org/html/2605.11181#bib.bib22); Loshchilov and Hutter,[2019](https://arxiv.org/html/2605.11181#bib.bib21)\)have long served as the dominant baselines, recent years have seen a surge of interest in matrix\-based and spectral optimizers\(Guptaet al\.,[2018](https://arxiv.org/html/2605.11181#bib.bib20);[Vyaset al\.,](https://arxiv.org/html/2605.11181#bib.bib19); Jordanet al\.,[2024](https://arxiv.org/html/2605.11181#bib.bib58)\)\.

Algorithms such asShampooandMuonhave demonstrated remarkable success at scale\(Liuet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib54); Teamet al\.,[2025a](https://arxiv.org/html/2605.11181#bib.bib45),[b](https://arxiv.org/html/2605.11181#bib.bib55); DeepSeek\-AI,[2026](https://arxiv.org/html/2605.11181#bib.bib18)\), sparking a renaissance in the study of non\-Euclidean descent methods\([Bernstein and Newhouse,](https://arxiv.org/html/2605.11181#bib.bib17)\)and an increasing suite of proposed modifications\(Du and Su,[2026](https://arxiv.org/html/2605.11181#bib.bib16); Ahnet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib15);[Riabininet al\.,](https://arxiv.org/html/2605.11181#bib.bib14); Siet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib13); Amselet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib12); Gonget al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib44)\)\. Theoretically, this success is almost universally justified through the lens of Linear Minimization Oracles \(LMOs\) and strict geometric preconditioning\(Pethicket al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib41); Kovalev,[2025](https://arxiv.org/html/2605.11181#bib.bib40); Fanet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib39)\)\. In this prevailing narrative,Muon’s performance stems from its exact adherence to a specific target geometry via complete spectrum whitening\. Analysis of simple models\(Kimet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib27); Maet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib28)\)commonly confirms this narrative, especially under distributional heavy tails\(Yuet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib29); Wanget al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib25)\)\.

However, recent benchmarking studies have begun to reveal cracks in this geometric facade, noting that advantages of spectral optimizers diminish under proper baseline tuning and are highly sensitive to batch size\(Wenet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib42); Semenovet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib43)\)\. Concurrently,Su \([2025](https://arxiv.org/html/2605.11181#bib.bib53)\)showed with an isotropic curvature model that optimal updates need only preserve the ordering of singular values, whereasMuon’s whitening imposes much stronger constraints\. Motivated by this, we initially set off to determine an optimal spectral geometry\. Yet, much like analyses ofMuonon quadratics\(Gononet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib24)\), we ultimately find no simple geometric recipe reliably improves performance\. Instead, we are forced to ask: does the exact geometric formulation of spectral optimizers matter, or is the LMO framework masking a more fundamental driver of optimization success?

In this work, we argue that while viewing optimization through the lens of non\-Euclidean geometry provides an elegant theoretical framework, it fails to capture the core phenomena underlying the performance of spectral optimizers in deep learning\. By systematically relaxing the assumptions of geometric preconditioning, we make the following specific contributions:

![Refer to caption](https://arxiv.org/html/2605.11181v1/x1.png)Figure 1:Transformations of the gradient singular value spectrum under various spectral optimizers for a GPT\-2 layer\. From left to right:SGDshows standard gradient spectrum\.Truncated SGDexplicitly zeros out the largest5%5\\%of singular values\.Muonmaps all singular values uniformly to1\.01\.0\.Kaonsets the values randomly\. Finally,Freonapplies normalized updates of the form\(G​G⊤\)−c​G\(GG^\{\\top\}\)^\{\-c\}G, smoothly interpolating the spectral decay between standard SGD \(c=0c=0\), Muon \(c=1/2c=1/2\), and the pseudoinverse\-like endpoint \(c=1c=1\)\. Note the logarithmic scale on the y\-axis\.1. 1\.To explore the broader space of order\-preserving updates, we introduceFreon, a family of Schatten \(quasi\-\)norm updates of the form\(G​G⊤\)−c​G\(GG^\{\\top\}\)^\{\-c\}G\.Freonnaturally interpolates between SGD \(c=0c=0\) and Muon \(c=1/2c=1/2\)\. Through a symmetry analysis analogous toYenet al\.\([2025](https://arxiv.org/html/2605.11181#bib.bib57)\), we argue that extrapolation into thec≥1/2c\\geq 1/2regime is of immense theoretical interest \([Theorem˜2\.6](https://arxiv.org/html/2605.11181#S2.Thmtheorem6)\), despite plunging into a quasi\-norm regime that fundamentally breaks standard LMO theory \([Theorem˜2\.1](https://arxiv.org/html/2605.11181#S2.Thmtheorem1)\)\. To compute these updates stably, we developed an optimal \([Theorems˜2\.7](https://arxiv.org/html/2605.11181#S2.Thmtheorem7)and[E\.5](https://arxiv.org/html/2605.11181#A5.Thmtheorem5)\) QDWH\-based iteration utilizing rational approximations \([Algorithm˜3](https://arxiv.org/html/2605.11181#algorithm3)\) similar toNakatsukasaet al\.\([2010](https://arxiv.org/html/2605.11181#bib.bib9)\); Nakatsukasa and Freund \([2016](https://arxiv.org/html/2605.11181#bib.bib8)\), extending the polar express theory ofAmselet al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib12)\)\. Empirically, unlike simple quadratics, we find that optimal exponents for a GPT\-2 model concentrate in the quasi\-norm regionc∈\[1/2,1\]c\\in\[1/2,1\]\([Figures˜I\.7](https://arxiv.org/html/2605.11181#A9.F7)and[6](https://arxiv.org/html/2605.11181#S3.F6)\), but varies significantly between batches\. As this is strictly outside the range of any unitarily invariant norm, it demonstrates conclusively that standard LMO theory cannot be enough to explain performance gains\.
2. 2\.Our empirical sweeps in[Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[4](https://arxiv.org/html/2605.11181#S2.F4)revealed a profound conundrum: almost allFreonupdates withc\>0c\>0performed remarkably similarly, provided large singular values were sufficiently suppressed\. This led us to hypothesize that the exact geometric structure of the update is largely irrelevant\. To test this, we introduceKaon, an optimizer that computes no LMO, targets no specific Schatten norm, and simply replaces the gradient’s singular values with noise \([Algorithm˜2](https://arxiv.org/html/2605.11181#algorithm2)\)\. Despite possessing zero coherent geometry,Kaonclosely matchesMuon’s performance at scale \([Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[3](https://arxiv.org/html/2605.11181#S1.F3)\) and curiously retains classical convergence guarantees \([Theorem˜2\.5](https://arxiv.org/html/2605.11181#S2.Thmtheorem5)\)\. This shows that precise target spectra are largely irrelevant for optimization performance\.
3. 3\.If LMOs are discarded, how can we mathematically understand the gap between all these updates? Extending the framework ofDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\), we decompose the exact local Taylor expansion into two structurally revealing quantities:batch gradient alignment\(γk\\gamma\_\{k\}\) anddirectional descent potential\(Φk\\Phi\_\{k\}\)\. While this is only mechanistic in quadratic models \([Section˜3\.2](https://arxiv.org/html/2605.11181#S3.SS2)\), it cleanly exposes the fundamental tradeoff \([Section˜3\.3](https://arxiv.org/html/2605.11181#S3.SS3)\): different optimizers implicitly choose how much alignment to sacrifice for increased descent potential, but actually realizing the benefit of a largeΦk\\Phi\_\{k\}hinges on using step sizes that are tuned to exploit it\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/x2.png)Figure 2:NanoGPT learning rate sensitivity\.Final validation loss as the tuned learning rates are jointly scaled, averaged over three seeds with±2\\pm 2std error bars\. Black outlines mark the best learning rate per optimiser\.
![Refer to caption](https://arxiv.org/html/2605.11181v1/x3.png)Figure 3:NanoGPT validation curves\.Validation\-loss curves for Muon, Kaon, and Freon atc=2/3c=2/3andc=3/4c=3/4\. Lines are averages over three seeds, with±2\\pm 2std bands\. The y\-axis is clipped to\[3\.3,6\]\[3\.3,6\]for legibility\.

SetX0=G/‖G‖FX\_\{0\}=G/\\\|G\\\|\_\{F\}

for*t=1,2,…,Tt=1,2,\\ldots,T*do

At=Xt−1​Xt−1⊤A\_\{t\}=X\_\{t\-1\}X\_\{t\-1\}^\{\\top\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
Bt=b​At\+c​At2B\_\{t\}=bA\_\{t\}\+cA\_\{t\}^\{2\}\\phantom\{B^\{a\}\\,X\_\{t\-1\}\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
Xt=a​Xt−1\+Bt​Xt−1X\_\{t\}=aX\_\{t\-1\}\+B\_\{t\}X\_\{t\-1\}\\vphantom\{B^\{b\}\\,X\_\{t\-1\}\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
end for

return

XTX\_\{T\}

Algorithm 1Muon
SetX0=G/‖G‖FX\_\{0\}=G/\\\|G\\\|\_\{F\}

for*t=1,2,…,Tt=1,2,\\ldots,T*do

At=Xt−1​Xt−1⊤A\_\{t\}=X\_\{t\-1\}X\_\{t\-1\}^\{\\top\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
Bt=\(I−At\)2B\_\{t\}=\(I\-A\_\{t\}\)^\{2\}\\phantom\{B^\{a\}\\,X\_\{t\-1\}\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
Xt=4\.1⋅Bt​Xt−1X\_\{t\}=4\.1\\cdot B\_\{t\}X\_\{t\-1\}\\phantom\{B^\{b\}\\,X\_\{t\-1\}\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
end for

return

XT/1\.175X\_\{T\}/1\.175\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}

Algorithm 2Kaon
SetX0=G/‖G‖FX\_\{0\}=G/\\\|G\\\|\_\{F\}

Set

A0=X0​X0⊤A\_\{0\}=X\_\{0\}X\_\{0\}^\{\\top\}
for*t=1,2,…,Tt=1,2,\\ldots,T*do

B=Rt​\(At−1\)B\\phantom\{\{\}\_\{t\}\}=R\_\{t\}\(A\_\{t\-1\}\)\\phantom\{X\_\{t\-1\}X\_\{t\-1\}^\{\\top\}\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
Xt=Ba​Xt−1X\_\{t\}=B^\{a\}\\,X\_\{t\-1\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
At=Bb​At−1A\_\{t\}=B^\{b\}\\,A\_\{t\-1\}\\vphantom\{\(A\_\{t\}\)^\{p\}\_\{t\}\}
end for

μ=\(1n​⟨XT,G⟩\)\(a\+2​b\)/\(2​a−2​b\)\\mu=\\bigl\(\\frac\{1\}\{n\}\\langle X\_\{T\},G\\rangle\\bigr\)^\{\(a\+2b\)/\(2a\-2b\)\}

return

XT/μX\_\{T\}/\\mu

Algorithm 3Freon

Algorithms: Simplified pseudocode illustrating the main differences between the optimizers considered\. Full algorithm for Freon is shown in[Algorithm˜4](https://arxiv.org/html/2605.11181#algorithm4), utilizing block\-QR for the rational map\.
## 2Background and Methods

We begin by considering the classical approach to introducingMuon, based on the idea of regularized steepest descent\. Given a norm∥⋅∥\\\|\\cdot\\\|on an inner product space \(spectral norm in the case ofMuon\), we define the Linear Minimization Oracle \(LMO\) and the corresponding dual norm as:

lmo∥⋅∥⁡\(V\)=arg⁡max‖U‖≤1​⟨U,V⟩‖V‖∗=max‖U‖=1​⟨U,V⟩=⟨V,lmo∥⋅∥⁡\(V\)⟩\.\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(\{V\}\)=\\underset\{\\\|\{U\}\\\|\\leq 1\}\{\\arg\\max\\;\}\\langle\{U\},\{V\}\\rangle\\qquad\\\|\{V\}\\\|\_\{\*\}=\\underset\{\\\|\{U\}\\\|=1\}\{\\max\}\\langle U,V\\rangle=\\langle V,\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(V\)\\rangle\.\(1\)Then, for a differentiable loss functionf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\rightarrow\\mathbb\{R\}, withGk=∇f​\(Xk\)G\_\{k\}=\\nabla f\(X\_\{k\}\), regularized steepest descent updates by minimizing a first\-order approximation of the loss under a quadratic penalty:

Xk\+1=arg⁡min𝑋​\{f​\(Xk\)\+⟨Gk,X−Xk⟩\+12​η​‖X−Xk‖2\}=Xk−η​‖Gk‖∗​lmo∥⋅∥⁡\(Gk\)\.\{X\}\_\{k\+1\}=\\underset\{\{X\}\}\{\\arg\\min\}\\left\\\{f\(\{X\}\_\{k\}\)\+\\left\\langle\{G\}\_\{k\},\{X\}\-\{X\}\_\{k\}\\right\\rangle\+\\frac\{1\}\{2\\eta\}\\left\\\|\{X\}\-\{X\}\_\{k\}\\right\\\|^\{2\}\\right\\\}=\{X\}\_\{k\}\-\\eta\\left\\\|\{G\}\_\{k\}\\right\\\|\_\{\*\}\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\\left\(\{G\}\_\{k\}\\right\)\.For example, forMuon,lmo∥⋅∥⁡\(W\)=\(W​W⊤\)−1/2​W\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(W\)=\(WW^\{\\top\}\)^\{\-1/2\}W, where the inverse is taken in the sense of a Moore\-Penrose pseudo\-inverse\. Unless otherwise specified, the proofs of theorems in this section can be found in[Appendix˜D](https://arxiv.org/html/2605.11181#A4)\.

### 2\.1Limitations of the LMO Framework

For simplicity, in this section we restrict our attention to a single layer, with a matrix domainℝm×n\\mathbb\{R\}^\{m\\times n\}\. We consider a differentiable loss functionf:ℝm×n→ℝf:\\mathbb\{R\}^\{m\\times n\}\\to\\mathbb\{R\}\. Letr=min⁡\(m,n\)r=\\min\(m,n\)denote the maximal rank, and write the singular value decompositionGk=∇f​\(Xk\)=Uk​diag\(σk\)⁡Vk⊤G\_\{k\}=\\nabla f\(X\_\{k\}\)=U\_\{k\}\\operatorname\*\{diag\}\(\\sigma\_\{k\}\)V\_\{k\}^\{\\top\}with, for simplicity, strictly positive singular valuesσk∈ℝ\>0r\\sigma\_\{k\}\\in\\mathbb\{R\}\_\{\>0\}^\{r\}\. We define Preconditioned Spectral Descent updates as updates of the formXk\+1=Xk−αk​⟨Gk,Dk⟩​DkX\_\{k\+1\}=X\_\{k\}\-\\alpha\_\{k\}\\langle G\_\{k\},D\_\{k\}\\rangle D\_\{k\}withDk=Uk​diag⁡\(pk​\(σk\)\)​Vk⊤D\_\{k\}=U\_\{k\}\\operatorname\{diag\}\(p\_\{k\}\(\\sigma\_\{k\}\)\)V\_\{k\}^\{\\top\}for some sequence of vector valued mappingspk:ℝ\>0r→ℝ\>0rp\_\{k\}:\\mathbb\{R\}\_\{\>0\}^\{r\}\\to\\mathbb\{R\}\_\{\>0\}^\{r\}\.

The two simplest members of this family are gradient descent \(GD\) andMuon\(spectral descent\)\. GD setspk​\(S\)=S/‖S‖2p\_\{k\}\(S\)=S/\\\|S\\\|\_\{2\}, whileMuonsetspk​\(S\)=1p\_\{k\}\(S\)=\{1\}, practically achieved using Newton–Schulz iterations \([Algorithm˜1](https://arxiv.org/html/2605.11181#algorithm1)\)\. However, not all preconditioned spectral descent updates are steepest descent updates:Freonforc\>1/2c\>1/2,TruncatedSGDandKaonare not, based on the following:

###### Theorem 2\.1\(Preconditioned Spectral Descent vs\. Regularized Steepest Descent\)\.

Fix the functionffand the pointXkX\_\{k\}\. A regularized steepest descent update w\.r\.t\. a unitarily invariant matrix norm fromXkX\_\{k\}is a preconditioned spectral descent update\. On the other hand, a preconditioned spectral descent update fromXkX\_\{k\}is a steepest descent update w\.r\.t\. a unitarily invariant matrix norm only ifpkp\_\{k\}preserves the order≤\\leqof the entries ofσk\\sigma\_\{k\}\.

###### Note 2\.2\.

This theorem exposes a fundamental limitation of LMOs; maximizing descent cannot remove ‘bad’ large or intermediate singular directions from a batched gradient\.

But, as we show, LMOs are not necessary for convergence of preconditioned spectral descent\.

###### Theorem 2\.3\(Convergence of Preconditioned Spectral Descent\)\.

Let∥⋅∥\\\|\\cdot\\\|be a unitarily invariant matrix norm, and let∥⋅∥∗\\\|\\cdot\\\|\_\{\*\}be its dual norm\. Letf:ℝm×n→ℝf:\\mathbb\{R\}^\{m\\times n\}\\to\\mathbb\{R\}be continuously differentiable, bounded below byfminf\_\{\\min\}, withL𝒳L\_\{\\mathcal\{X\}\}\-Lipschitz continuous gradients with respect to∥⋅∥\\\|\\cdot\\\|\.

Consider preconditioned spectral descent updates as in[Section˜2\.1](https://arxiv.org/html/2605.11181#S2.SS1), and assume there exist constantsmk≥0m\_\{k\}\\geq 0andMk\>0M\_\{k\}\>0satisfying the divergence condition∑k=0∞\(mk/Mk\)2=∞\\sum\_\{k=0\}^\{\\infty\}\\left\(\{m\_\{k\}\}/\{M\_\{k\}\}\\right\)^\{2\}=\\infty, such that for allσ∈ℝ\>0r\\sigma\\in\\mathbb\{R\}\_\{\>0\}^\{r\}, the mappingpkp\_\{k\}satisfies the following sufficient descent and boundedness conditions:

⟨σ,pk​\(σ\)⟩≥mk​‖σ‖∗and‖pk​\(σ\)‖≤Mk\.\\langle\\sigma,p\_\{k\}\(\\sigma\)\\rangle\\geq m\_\{k\}\\\|\\sigma\\\|\_\{\*\}\\qquad\\text\{and\}\\qquad\\\|p\_\{k\}\(\\sigma\)\\\|\\leq M\_\{k\}\.Then, choosingαk=η​1L𝒳​Mk2\\alpha\_\{k\}=\\eta\\frac\{1\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}for anyη∈\(0,2\)\\eta\\in\(0,2\), the following convergence rate holds forK≥1K\\geq 1:

min0≤k<K⁡‖Gk‖∗2≤L𝒳​\(f​\(X0\)−fmin\)η​\(1−η2\)​∑k=0K−1\(mkMk\)2andlim infk→∞‖Gk‖∗=0\.\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}\\leq\\frac\{L\_\{\\mathcal\{X\}\}\(f\(X\_\{0\}\)\-f\_\{\\min\}\)\}\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\}\\qquad\\text\{and\}\\qquad\\liminf\_\{k\\to\\infty\}\\\|G\_\{k\}\\\|\_\{\*\}=0\.

###### Note 2\.4\.

The proof relies on standard convergence arguments, requiring that updates do not become too small, do not blow up and are mostly in the right direction\. Whenpkp\_\{k\}is derived from anlmo∥⋅∥\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}, these constants evaluate tomk=Mk=1m\_\{k\}=M\_\{k\}=1\. Crucially, this theorem accommodates non\-norm\-descent methods likeKaonandFreon\. While this result does not yield a tight convergence rate, it provides an essential sanity check: such methods will work in practice\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/x4.png)Figure 4:Final validation loss versus learning rate for all optimizer families on WikiText\-2 \(118M tokens\)\.\(a\)TruncatedSGDimproves overSGDbut fails to close the gap toMuon– suppressing large singular values is necessary but not sufficient\.\(b\)Kaonclosely matchesMuondespite replacing singular values with random noise\.\(c\)Freonwithc≈2/3c\\approx 2/3closely matchesMuon; the optimalcclies strictly outside the range of any proper Schatten norm \(c\>1/2c\>1/2\)\.
### 2\.2TruncatedSGD

As a first step in relaxing the strict geometric assumptions of the LMO framework, we isolate the role of the largest singular values in optimization stability\. We studyTruncatedSGD, an optimizer that zeros out the toppp% of singular values while leaving the rest unchanged\. Although a full SVD is impractical,TruncatedSGDserves a pedagogical purpose: it tests whether standardSGDfails primarily due to instability induced by excessively large singular values\.

Concurrently with our work,Jianget al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib37)\)propose a Newton\-Schulz based approximation to fixed\-value clipping, similarly targeting large singular values under the assumption that batch noise concentrates in this part of the spectrum\. In line with this motivation, we observe in[Figure˜4](https://arxiv.org/html/2605.11181#S2.F4)thatTruncatedSGDimproves overSGD, confirming that suppressing large singular values is*necessary*for stability\. However, it still falls short ofMuon, showing that suppression alone is*insufficient*\. The view that clipping primarily removes noise implicitly treats the underlying gradient direction as intrinsically meaningful; our results in[Section˜3\.3](https://arxiv.org/html/2605.11181#S3.SS3)ultimately challenge this premise, suggesting instead that the gradient direction itself can be poorly aligned with useful progress\.

### 2\.3Kaon

Motivated byTruncatedSGD’s gains overSGDand observation of[Section˜2\.4](https://arxiv.org/html/2605.11181#S2.SS4)that only the broad suppression/amplification of singular values matters, we ask: does the target geometry matter at all? To test this, we replaced the gradient’s singular values with uniform random noise\. This stochastic spectrum redistribution still performed surprisingly well, closely matchingMuon\.

To make this ‘absurd optimizer’ efficient, we draw on ideas for pseudorandom number generation using chaotic iterations\(Phatak and Rao,[1995](https://arxiv.org/html/2605.11181#bib.bib30)\)\. We use the generalized third\-order logistic recurrencext\+1=λ​xt​\(1−xt2\)2,x\_\{t\+1\}=\\lambda x\_\{t\}\(1\-x\_\{t\}^\{2\}\)^\{2\},lifted to matrix iterations \([Algorithm˜2](https://arxiv.org/html/2605.11181#algorithm2)\)\. Settingλ=4\.1\\lambda=4\.1pushes the recurrence deep into the chaotic regime\. We visualize an example in[Figure˜1](https://arxiv.org/html/2605.11181#S1.F1), and illustrate the stationary distribution of these iterations in[Section˜F\.1](https://arxiv.org/html/2605.11181#A6.SS1),[Figure˜1\(a\)](https://arxiv.org/html/2605.11181#A6.F1.sf1)\.

We emphasize thatKaonis primarily a pedagogical construction\. Its strong performance \([Figure˜4](https://arxiv.org/html/2605.11181#S2.F4)\) is a direct counter\-argument to LMO\-based explanations of spectral optimizers: state\-of\-the\-art performance does not require adherence to a precise target geometry, and can even be achieved with random spectral reshaping, while retaining standard convergence guarantees:

###### Theorem 2\.5\(Convergence of Random Spectral Descent\)\.

Let the setting be identical to[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)\. Suppose at each iterationkk, the stochastic mapping is defined specifically aspk​\(σ\)=Ek‖Ek‖p\_\{k\}\(\\sigma\)=\\frac\{E\_\{k\}\}\{\\\|E\_\{k\}\\\|\}, where the entries ofEkE\_\{k\}are i\.i\.d\. samples of a probability distribution onℝ\>0\\mathbb\{R\}\_\{\>0\}independent ofkk\.

Ifαk=ηL𝒳\\alpha\_\{k\}=\\frac\{\\eta\}\{L\_\{\\mathcal\{X\}\}\}for any constantη∈\(0,2\)\\eta\\in\(0,2\), then almost surely, the sequencemin0≤k<K⁡‖Gk‖∗2\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}achieves a convergence rate of𝒪​\(1/K\)\\mathcal\{O\}\(1/K\)\.

### 2\.4Freon

Following the setup of[Section˜2](https://arxiv.org/html/2605.11181#S2), we consider steepest descent in Schattenpp\-norms, defined as‖X‖p=\(∑σi​\{X\}p\)1/p\\\|X\\\|\_\{p\}=\\left\(\\sum\\sigma\_\{i\}\\\{X\\\}^\{p\}\\right\)^\{1/p\}\. By Hölder’s inequality, the dual of thepp\-Schatten norm is theqq\-Schatten norm\(Bhatia,[2013](https://arxiv.org/html/2605.11181#bib.bib4)\), where1p\+1q=1\\frac\{1\}\{p\}\+\\frac\{1\}\{q\}=1forp∈\[1,∞\]p\\in\[1,\\infty\]\. As before, let the singular value decomposition of the gradient beG=U​diag⁡\(σ\)​V⊤G=U\\operatorname\{diag\}\(\\sigma\)V^\{\\top\}\. The steepest descent direction forp∈\[1,∞\]p\\in\[1,\\infty\]is:

lmo∥⋅∥p⁡\(G\)=U​diag⁡\(\(σ‖G‖q\)q−1\)​V⊤or equivalentlyD=\(Gn​Gn⊤\)−c​Gn,\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{p\}\}\(G\)=U\\operatorname\{diag\}\\left\(\\left\(\\frac\{\\sigma\}\{\\\|G\\\|\_\{q\}\}\\right\)^\{q\-1\}\\right\)\{V\}^\{\\top\}\\quad\\text\{or equivalently\}\\quad D=\\left\(G\_\{n\}G\_\{n\}^\{\\top\}\\right\)^\{\-c\}G\_\{n\},\(2\)where we denote the dual\-normalized gradient asGn=G/‖G‖qG\_\{n\}=G/\\\|G\\\|\_\{q\}andc=1−q/2c=1\-q/2\. This update structure, parametrized by the dual exponentqq, forms the basis of the generalFreonoptimizer\.

##### A note onp∈\(0,1\)p\\in\(0,1\):

While the derived update rule extends analytically top∈\(0,1\)p\\in\(0,1\)\(whereq<0q<0\), thepp\-Schatten ceases to be a proper norm in this regime, becoming aquasi\-normdue to the violation of the triangle inequality\. Strictly speaking, the dual of app\-quasi\-norm forp<1p<1is the∞\\infty\-norm \(spectral norm\)\. Consequently, exact steepest descent under a truepp\-quasi\-norm constraint does not yield the continuous deformation above \([Theorem˜2\.1](https://arxiv.org/html/2605.11181#S2.Thmtheorem1)prevents that\); rather, the LMO collapses, activating only the singular vector of the largest singular value\. Therefore, the generalFreonupdate withq<0q<0does not correspond to a steepest descent under the quasi\-norm either\. It is better viewed as a preconditioned spectral descent \(c\.f\.[Section˜2\.1](https://arxiv.org/html/2605.11181#S2.SS1)\)\.

#### 2\.4\.1Norm Hyperparameter Transfer

If one were to execute the exact update derived above with the same step size across different values ofpp, the optimization would fail dramatically\. This is because restricting the update to the standard unitpp\-ball\(‖Δ​W‖p≤1\)\(\\\|\\Delta W\\\|\_\{p\}\\leq 1\)causes the maximum achievable step size in the spectral space to shrink at a rate dependent onpp, as the rankr=min⁡\(m,n\)r=\\min\(m,n\)increases\. To anchor the updates to a consistent scale regardless ofpp, we must tie the updates to a single overarching invariant ball: the spectral ball\. We achieve this by scaling the update constraint precisely byr1/pr^\{1/p\}, enforcing‖Δ​W‖p=r1/p\.\\\|\\Delta W\\\|\_\{p\}=r^\{1/p\}\.This means that for anyp\>0p\>0,1r1/p​‖Δ​W‖p≤‖Δ​W‖∞\\frac\{1\}\{r^\{1/p\}\}\\\|\\Delta W\\\|\_\{p\}\\leq\\\|\\Delta W\\\|\_\{\\infty\}\. By shifting to this mean Schatten norm, defined by∥⋅∥p,mean=1r1/p∥⋅∥p\\\|\\cdot\\\|\_\{p,\\texttt\{mean\}\}=\\frac\{1\}\{r^\{1/p\}\}\\\|\\cdot\\\|\_\{p\}, we geometrically inflate thepp\-Schatten ball, guaranteeing that the baseline maximum step size is always𝒪​\(1\)\\mathcal\{O\}\(1\)independent ofpp, as domain of thelmo∥⋅∥\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}:\{‖D‖p,mean≤1\}⊆\{‖D‖∞≤1\}\\\{\\\|D\\\|\_\{p,\\texttt\{mean\}\}\\leq 1\\\}\\subseteq\\\{\\\|D\\\|\_\{\\infty\}\\leq 1\\\}for allp≥0p\\geq 0\. We also note, that concurrently with our work, a similar normalization was introduced byXuet al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib31)\), restricted to the context of vector\-induced row\-wise and column\-wisep,qp,qnorms, based on arguments of dimension independence\.

#### 2\.4\.2Particular case ofp=0p=0\(c=1c=1\)

The mean\-norm scaling introduced in[Section˜2\.4\.1](https://arxiv.org/html/2605.11181#S2.SS4.SSS1)does more than just stabilize hyperparameters; it mathematically enables us to explore the limiting case wherep→0p\\to 0\(andq→0q\\to 0\)\. When combined with the normalization from[Section˜2\.4\.1](https://arxiv.org/html/2605.11181#S2.SS4.SSS1)the scaling in[Equation˜2](https://arxiv.org/html/2605.11181#S2.E2)converges:∥G∥q/r1/q→det\(GG⊤\)1/2​n\\\|G\\\|\_\{q\}/r^\{1/q\}\\xrightarrow\[\]\{\}\\operatorname\{det\}\(GG^\{\\top\}\)^\{1/2n\}asq→0q\\to 0\. The beauty of thep=0p=0limit extends far beyond this simple closed form; up to scaling, this specific update turns out to be invariant under layerwise symmetries\.

###### Theorem 2\.6\(Equivariance ofFreon\(c=1\)\(c=1\)\. Informal, formal statement in[Theorem˜D\.2](https://arxiv.org/html/2605.11181#A4.Thmtheorem2)\)\.

Consider an arbitrary neural networkϕW\\phi\_\{W\}with layered parametersW=\(W1,…,WL\)W=\(W\_\{1\},\\dots,W\_\{L\}\)\. We define a symmetry ofϕW\\phi\_\{W\}to be any transformationW↦g⋅WW\\mapsto g\\cdot Wsuch thatϕW=ϕg⋅W\\phi\_\{W\}=\\phi\_\{g\\cdot W\}\. Then updates of the formWl−η​\(Gl​Gl⊤\)−1​GlW\_\{l\}\-\\eta\(G\_\{l\}G\_\{l\}^\{\\top\}\)^\{\-1\}G\_\{l\}for each layerWlW\_\{l\}are equivariant under scaling and orthogonal layerwise symmetries ofϕW\\phi\_\{W\}, i\.e\.g⋅Wk=\[g⋅W\]kg\\cdot W^\{k\}=\[g\\cdot W\]^\{k\}, whereWkW^\{k\}represents the parameters afterkkupdates\. Therefore, the optimization trajectory of the neural networkϕW\\phi\_\{W\}is invariant, i\.e\.ϕ\[g⋅W\]k=ϕWk\\phi\_\{\[g\\cdot W\]^\{k\}\}=\\phi\_\{W^\{k\}\}\.

#### 2\.4\.3Efficient Computation: From Polar to Hölder Express

To implementFreonpractically, we require an efficient method to compute terms of the form\(G​G⊤\)−ab​G\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}Gfrom[Equation˜2](https://arxiv.org/html/2605.11181#S2.E2), whereq=2​\(1−ab\)q=2\(1\-\\frac\{a\}\{b\}\)andab∈ℚ\\frac\{a\}\{b\}\\in\\mathbb\{Q\}\.

##### On polynomial iterations:

A natural approach is a polynomial scheme akin to\(Jordanet al\.,[2024](https://arxiv.org/html/2605.11181#bib.bib58); Amselet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib12)\)\. Extending the Polar Express theory \([Theorem˜E\.5](https://arxiv.org/html/2605.11181#A5.Thmtheorem5)\), we find optimal polynomialsPkP\_\{k\}such thatx​Pk​\(xb\)→1xP\_\{k\}\(x^\{b\}\)\\to 1under composition \(Appendix[E\.1](https://arxiv.org/html/2605.11181#A5.SS1), Remark[E\.6](https://arxiv.org/html/2605.11181#A5.Thmtheorem6)\), which can be applied iteratively via updates of the formOk\+1=Pk​\(\(Ok​Ok⊤\)b/2​\(G​G⊤\)a−b/2\)​OkO\_\{k\+1\}=P\_\{k\}\\left\(\(O\_\{k\}O\_\{k\}^\{\\top\}\)^\{b/2\}\(GG^\{\\top\}\)^\{a\-b/2\}\\right\)O\_\{k\}\. While theoretically convergent \(Appendix[E\.1](https://arxiv.org/html/2605.11181#A5.SS1), Remark[E\.9](https://arxiv.org/html/2605.11181#A5.Thmtheorem9)\), these iterations are highly unstable in practice fora/b≠1/2a/b\\neq 1/2\([Section˜E\.3](https://arxiv.org/html/2605.11181#A5.SS3)\)\. Requiring a computationally prohibitive∼𝒪​\(b​log⁡σmin\)\\sim\\mathcal\{O\}\(b\\log\{\\sigma\_\{\\min\}\}\)steps, they diverge rapidly, especially in lower precision\. This instability is well\-known for polynomial iterations outside ofa/b=1/2a/b=1/2\(Higham,[2008](https://arxiv.org/html/2605.11181#bib.bib36)\)\. In line withZhanget al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib35)\), we similarly identify two primary causes: the explicit formation of poorly conditioned matrices \(which introduces spurious negative eigenvalues\) and the rapid accumulation of floating\-point errors\. A robust method must therefore avoid explicitly forming such matrices and converge in as few steps as possible\.

##### On rational iterations:

Generalized rational functions, rooted in Zolotarev’s best rational approximants to the sign function\(Achieser,[1992](https://arxiv.org/html/2605.11181#bib.bib51); Nakatsukasa and Freund,[2016](https://arxiv.org/html/2605.11181#bib.bib8)\), offer significantly greater stability\. For instance, the QR\-based Dynamically Weighted Halley \(QDWH\) method\(Nakatsukasaet al\.,[2010](https://arxiv.org/html/2605.11181#bib.bib9); Nakatsukasa and Freund,[2016](https://arxiv.org/html/2605.11181#bib.bib8)\)uses these to approximate the polar map\. Drawing on this classical theory\(Achieser,[1992](https://arxiv.org/html/2605.11181#bib.bib51); Cheney and Loeb,[1964](https://arxiv.org/html/2605.11181#bib.bib49)\)and use in matrix computations\(Nakatsukasa and Freund,[2016](https://arxiv.org/html/2605.11181#bib.bib8); Gawlik and Nakatsukasa,[2021](https://arxiv.org/html/2605.11181#bib.bib33)\), we use the Remez algorithm\(Filipet al\.,[2018](https://arxiv.org/html/2605.11181#bib.bib6)\)to iteratively approximate the sign function using rational functions of the formx​R​\(x2​b\)=x​a\+b​x2​b1\+c​x2​bxR\(x^\{2b\}\)=x\\frac\{a\+bx^\{2b\}\}\{1\+cx^\{2b\}\}akin toAmselet al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib12)\), establishing optimality of this approach:

###### Theorem 2\.7\(Proof in[Section˜E\.1](https://arxiv.org/html/2605.11181#A5.SS1)\)\.

Leta/b∈ℚa/b\\in\\mathbb\{Q\},l\>0l\>0and supposeGGsatisfiesσmin​\(G\)≥l\\sigma\_\{\\min\}\(G\)\\geq l\. Moreover, set𝒢:=\(G​G⊤\)−ab​G\\mathcal\{G\}:=\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}G,𝒢†:=G⊤​\(G​G⊤\)ab−1\\mathcal\{G\}^\{\\dagger\}:=G^\{\\top\}\(GG^\{\\top\}\)^\{\\frac\{a\}\{b\}\-1\}andℜn,d​\(b2,b2\)\\mathfrak\{R\}\_\{n,d\}\\left\(\\frac\{b\}\{2\},\\frac\{b\}\{2\}\\right\)as in \([15](https://arxiv.org/html/2605.11181#A5.E15)\) withn,d∈ℕn,d\\in\\mathbb\{N\}such thatn\+1d\+1∈ℕ\\frac\{n\+1\}\{d\+1\}\\in\\mathbb\{N\}\. Then, there exists\{Rt\}t∈ℕ⊂ℜn,d​\(b2,b2\)\\\{R\_\{t\}\\\}\_\{t\\in\\mathbb\{N\}\}\\subset\\mathfrak\{R\}\_\{n,d\}\\left\(\\frac\{b\}\{2\},\\frac\{b\}\{2\}\\right\)such that, settingOT:=R~T∘⋯∘R~1O\_\{T\}:=\\widetilde\{R\}\_\{T\}\\circ\\cdots\\circ\\widetilde\{R\}\_\{1\}, withR~t​\(X\)=Rt​\(X​𝒢†\)​𝒢\\widetilde\{R\}\_\{t\}\(X\)=R\_\{t\}\(X\\mathcal\{G\}^\{\\dagger\}\)\\mathcal\{G\},OT​\(G\)O\_\{T\}\(G\)optimally approximates𝒢\\mathcal\{G\}for everyT∈ℕT\\in\\mathbb\{N\}\. Moreover, there existsC=C​\(a,b,l\)\>0C=C\(a,b,l\)\>0such that the following doubly exponential convergence rate holds:

‖𝒢−OT​\(G\)‖2≤C​\|1−lb\|\(n\+d\+1\)T\.\\\|\\mathcal\{G\}\-O\_\{T\}\(G\)\\\|\_\{2\}\\leq C\|1\-l^\{b\}\|^\{\(n\+d\+1\)^\{T\}\}\.

In practice,OTO\_\{T\}is computed using the coupled system presented in[Algorithm˜3](https://arxiv.org/html/2605.11181#algorithm3), with the complete algorithm summarized in[Algorithm˜4](https://arxiv.org/html/2605.11181#algorithm4)\. Unlike prior work\(Gawlik and Nakatsukasa,[2021](https://arxiv.org/html/2605.11181#bib.bib33)\), wherein the considered class of rational functions does not contain the best approximant forp≥3p\\geq 3, our recursion guarantees the existence of an optimal approximant at every step\. Moreover, the condition onn,dn,dpermits choosing higher\-degree polynomials in the numerator without increasing the denominator’s degree, and thus allowing block\-QR tricks\.

## 3Theoretical and Empirical Observations

If LMOs do not dictate optimization performance, what does? We build uponDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\), which provides a foundational theoretical lens for understanding spectral descent\. To properly diagnose these failures, we must abandon global bounds \(as inCohenet al\.\([2021](https://arxiv.org/html/2605.11181#bib.bib3)\); Islamovet al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib5)\)\) and examine the exact local expansion of the loss\.

### 3\.1The Two Quantities

Letf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}be twice differentiable, and suppose we updateXk\+1=Xk−αk​ΔkX\_\{k\+1\}=X\_\{k\}\-\\alpha\_\{k\}\\Delta\_\{k\}, withΔk=⟨G~k,Dk⟩​Dk\\Delta\_\{k\}=\\langle\\widetilde\{G\}\_\{k\},D\_\{k\}\\rangle D\_\{k\},G~k\\widetilde\{G\}\_\{k\}is a stochastic gradient offf, andDkD\_\{k\}is the direction given by the optimizer, defined for example through preconditioned spectral descent \([Section˜2\.1](https://arxiv.org/html/2605.11181#S2.SS1)\)\. By Taylor’s theorem,f​\(Xk\+1\)=f​\(Xk\)−αk​⟨Gk,Δk⟩\+αk22​⟨Δk,∇2f​\(Zk\)​\[Δk\]⟩,f\(X\_\{k\+1\}\)=f\(X\_\{k\}\)\-\\alpha\_\{k\}\\langle G\_\{k\},\\Delta\_\{k\}\\rangle\+\\frac\{\\alpha\_\{k\}^\{2\}\}\{2\}\\langle\\Delta\_\{k\},\\nabla^\{2\}f\(Z\_\{k\}\)\[\\Delta\_\{k\}\]\\rangle,withZkZ\_\{k\}on the line segment betweenXkX\_\{k\}andXk\+1X\_\{k\+1\}\. Plugging in our updateΔk=⟨G~k,Dk⟩​Dk\\Delta\_\{k\}=\\langle\\widetilde\{G\}\_\{k\},D\_\{k\}\\rangle D\_\{k\}, and factoring inherently isolates two fundamental quantities that govern the trajectory of the optimizer:

γk:=⟨Gk,Dk⟩⟨G~k,Dk⟩,⏟Batch Gradient AlignmentΦk:=⟨G~k,Dk⟩2⟨Dk,∇2f​\(Zk\)​\[Dk\]⟩⏟Local Directional Descent Potential\.\\underbrace\{\\gamma\_\{k\}:=\\frac\{\\langle G\_\{k\},D\_\{k\}\\rangle\}\{\\langle\\widetilde\{G\}\_\{k\},D\_\{k\}\\rangle\},\}\_\{\\text\{Batch Gradient Alignment\}\}\\qquad\\underbrace\{\\Phi\_\{k\}:=\\vphantom\{\\frac\{\\langle G\_\{k\},D\_\{k\}\\rangle\}\{\\langle\\widetilde\{G\}\_\{k\},D\_\{k\}\\rangle\}\}\\frac\{\\langle\\widetilde\{G\}\_\{k\},D\_\{k\}\\rangle^\{2\}\}\{\\langle D\_\{k\},\\nabla^\{2\}f\(Z\_\{k\}\)\[D\_\{k\}\]\\rangle\}\}\_\{\\text\{Local Directional Descent Potential\}\}\.\(3\)
Define in additionλk=⟨Dk,∇2f​\(Zk\)​\[Dk\]⟩\\lambda\_\{k\}=\\langle D\_\{k\},\\nabla^\{2\}f\(Z\_\{k\}\)\[D\_\{k\}\]\\rangle, a quantity which will only appear through multiplication with the learning rateαk\\alpha\_\{k\}\. Substituting these into the above yields the exact descent condition:

f​\(Xk\+1\)−f​\(Xk\)=−Φk​\(γk−\[αk​λk\]2\)​\[αk​λk\]\.f\(X\_\{k\+1\}\)\-f\(X\_\{k\}\)=\-\\Phi\_\{k\}\\left\(\\gamma\_\{k\}\-\\frac\{\[\\alpha\_\{k\}\\lambda\_\{k\}\]\}\{2\}\\right\)\[\\alpha\_\{k\}\\lambda\_\{k\}\]\.\(4\)[Equation˜4](https://arxiv.org/html/2605.11181#S3.E4)is a quadratic inαk\\alpha\_\{k\}\(assumingλk\\lambda\_\{k\}is approximately constant\), minimised atαk∗=γkλk\\alpha\_\{k\}^\{\*\}=\\frac\{\\gamma\_\{k\}\}\{\\lambda\_\{k\}\}\. At this optimum,[Equation˜4](https://arxiv.org/html/2605.11181#S3.E4)does not depend onλk\\lambda\_\{k\}, andΔ​f=−12​Φk​γk2\\Delta f=\-\\frac\{1\}\{2\}\\Phi\_\{k\}\\gamma\_\{k\}^\{2\}\. There are a number of observations we should make about this:

1. 1\.The derivation of[Equation˜4](https://arxiv.org/html/2605.11181#S3.E4)does not utilize anything apart from smoothness\. Crucially, neither quantity invokes an LMO or places any restriction on howDkD\_\{k\}is constructed\. While we cannot directly prove convergence without further assumptions \(like[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)\), it still remains useful for explaining multiple phenomena, as explained below\.
2. 2\.This directly recovers the formula ofDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\)as a special case, by boundingλk\\lambda\_\{k\}with a Lipschitz constant, settingG~k=Gk\\widetilde\{G\}\_\{k\}=G\_\{k\}andDk=lmo∥⋅∥⁡\(Gk\)D\_\{k\}=\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(G\_\{k\}\), givesγk=1\\gamma\_\{k\}=1with bound−‖Gk‖∗2/2​L∥⋅∥\-\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}/2L\_\{\\\|\\cdot\\\|\}\. But, compared with\(Davis and Drusvyatskiy,[2026](https://arxiv.org/html/2605.11181#bib.bib38)\), this includes three important generalizations: a non\-fixed step\-size, non\-LMO\-based updates, and stochasticity\.
3. 3\.While for any fixed norm, theDkD\_\{k\}maximizingΦk​λk\\Phi\_\{k\}\\lambda\_\{k\}islmo∥⋅∥⁡\(G~k\)\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(\\widetilde\{G\}\_\{k\}\), this may not be the direction that maximizesγk\\gamma\_\{k\}or minimizesλk\\lambda\_\{k\}\.
4. 4\.Because the intermediate pointZkZ\_\{k\}is unknown a priori,λk\\lambda\_\{k\}cannot be calculated during practical training\. Its exact analytical use requires quadratic models \(as explored in[Section˜3\.2](https://arxiv.org/html/2605.11181#S3.SS2)\)\.

Ultimately, optimization performance is driven by the interplay of these quantities \(γk,Φk\\gamma\_\{k\},\\Phi\_\{k\}, and the learning rateαk​λk\\alpha\_\{k\}\\lambda\_\{k\}\), rather than the strict form of the update direction itself\. In the remainder of this paper, we analyze these quantities exactly in the RF model and empirically for a GPT2 training run\.

### 3\.2Exact Asymptotics in a Random Feature Model

We take the well\-specified random feature regression model as inDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\)withWWthe weights andA∈ℝn×bA\\in\\mathbb\{R\}^\{n\\times b\}batched post\-activations of the previous layer:

minW∈ℝo×n⁡f​\(W\)=12​b​‖W​A−Y′​A‖F2\.\\min\_\{W\\in\\mathbb\{R\}^\{o\\times n\}\}f\(W\)=\\frac\{1\}\{2b\}\\\|WA\-Y^\{\\prime\}A\\\|^\{2\}\_\{F\}\.\(5\)We will consider the proportional asymptotic limit\(n,b→∞\(n,b\\rightarrow\\inftywithn/b→δ\)n/b\\rightarrow\\delta\)and assume the batched post\-activations are mean\-zero Gaussians,ai∼𝒩​\(0,C\)a\_\{i\}\\sim\\mathcal\{N\}\(0,C\), with a well\-defined limit forCC\. LetG=Y′−WG=Y^\{\\prime\}\-Wbe the oracle gradient\. Using standard random matrix theory, one can show almost sure convergence of bothΦ\\Phiandγ,\\gamma,and we give full details of this in[Appendix˜G](https://arxiv.org/html/2605.11181#A7)\. The general equations are hard to analyze, but in the case thatCCis aligned with the right singular vectors ofGG, i\.e\.,V⊤​C​V=ΛV^\{\\top\}CV=\\LambdaforΛ\\Lambdaa diagonal matrix, we can obtain certain illuminating results\.

###### Proposition 3\.1\.

If the entries ofΛ\\Lambdaare in non\-increasing order thenγSGD≤γMuon\\gamma\_\{\\operatorname\{SGD\}\}\\leq\\gamma\_\{\\operatorname\{Muon\}\}with equality if and only ifΛ=c​I\\Lambda=cIorΣ=c​I\\Sigma=cI\.

###### Proposition 3\.2\.

IfΛ=I\\Lambda=IthenΦSGD≥ΦMuon\\Phi\_\{\\operatorname\{SGD\}\}\\geq\\Phi\_\{\\operatorname\{Muon\}\}and thereforeγSGD2​ΦSGD≥γMuon2​ΦMuon\\gamma^\{2\}\_\{\\operatorname\{SGD\}\}\\Phi\_\{\\operatorname\{SGD\}\}\\geq\\gamma^\{2\}\_\{\\operatorname\{Muon\}\}\\Phi\_\{\\operatorname\{Muon\}\}

###### Theorem 3\.3\.

In the undersampled case \(δ→∞\\delta\\rightarrow\\infty\) we havelimδ→∞​γSGD2​ΦSGD≥limδ→∞​γMuon2​ΦMuon\\underset\{\\delta\\rightarrow\\infty\}\{\\lim\}\\gamma\_\{\\operatorname\{SGD\}\}^\{2\}\\Phi\_\{\\operatorname\{SGD\}\}\\geq\\underset\{\\delta\\rightarrow\\infty\}\{\\lim\}\\gamma\_\{\\operatorname\{Muon\}\}^\{2\}\\Phi\_\{\\operatorname\{Muon\}\}

![Refer to caption](https://arxiv.org/html/2605.11181v1/x5.png)Figure 5:Random\-feature regression with ReLU activations\.We consider the problem in[Equation˜5](https://arxiv.org/html/2605.11181#S3.E5), using the setup ofDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\)\.Left:training loss for GD, SpecGD, Kaon, their optimal\-step variants, and two adaptive\-exponent methods\.Centre:effective step sizes; the optimal\-GD step size highly oscillates, while optimal\-specGD step size is constant, being equal to specGD\.Right:the exponentccin the update direction\(G​G⊤\)−c​G\(GG^\{\\top\}\)^\{\-c\}Gas tracked by the two Optimal C methods during training\. Further details can be found in[Section˜H\.1](https://arxiv.org/html/2605.11181#A8.SS1)\.
### 3\.3Empirical Observations

##### The RF optimal step sizes:

The asymptotic results in[Section˜3\.2](https://arxiv.org/html/2605.11181#S3.SS2)compare update directions through local descent quantities, but they abstract away a key practical issue of step size stability\. Specifically, if we evaluate the optimal per\-step step size \(i\.e\., exact line search\), the conclusions ofDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\)completely reverse: standard GD achieves greater loss reduction than spectral descent, a phenomenon also noted byGononet al\.\([2026](https://arxiv.org/html/2605.11181#bib.bib24)\)\. However, optimal GD requires a highly oscillatory and practically unrealizable step\-size schedule \([Figure˜5](https://arxiv.org/html/2605.11181#S3.F5)\)\. Conversely, this exact line\-search model reveals the core mechanistic advantage of spectral descent: its optimal step size is constant\. Because this step size is proportional to1/‖D​A‖F21/\\\|DA\\\|^\{2\}\_\{F\}, whereDDis the polar factor ofGG, the value remains practically constant throughout training\.

##### The RF optimal norms:

Beyond step\-size stability, the RF model can be pushed further by dynamically selecting the Schatten norm exponentccthat maximizes the local directional potential\. We observe, however, that greedily optimizingccat every step yields rapid initial descent but degrades long\-term dynamics, even when paired with optimal step sizes \([Figure˜5](https://arxiv.org/html/2605.11181#S3.F5)\)\. Instead, by assuming a certain power\-law decay, we derive a simple, fixed\-exponent rule \([Section˜H\.2](https://arxiv.org/html/2605.11181#A8.SS2)\)\. This approximate choice significantly improves overall performance, overcoming the limitations of spectral descent for mean\-zero activations identified byDavis and Drusvyatskiy \([2026](https://arxiv.org/html/2605.11181#bib.bib38)\), as demonstrated in[Figure˜H\.1](https://arxiv.org/html/2605.11181#A8.F1)\.

##### The naive RF plug\-in for GPT\-2:

Attempting to map the RF model directly onto a GPT\-2 architecture exposes distinct limitations, as illustrated by the optimal exponent heatmaps in[Figure˜I\.5](https://arxiv.org/html/2605.11181#A9.F5)\(Left\) and[Figure˜I\.4](https://arxiv.org/html/2605.11181#A9.F4)\. The RF model fails to provide an accurate quadratic approximation of the layer\. Specifically, optimal exponents frequently fall into the strict non\-norm regime \(c\>1/2c\>1/2\), indicating that such updates fundamentally deviate from typical RF quadratic landscapes \([Figure˜I\.3](https://arxiv.org/html/2605.11181#A9.F3)\)\.

##### Single batch vs\. Full dataset geometry in GPT\-2:

![Refer to caption](https://arxiv.org/html/2605.11181v1/x6.png)Figure 6:Loss changeΔ​f\\Delta fover the joint\(c,η\)\(c,\\eta\)grid, evaluated on the full validation set\. Each panel shows a different training point used for the calculation of the gradient\. The dashed black curve marks the optimal learning rate percc; the green dotted line marks the actual training learning rate\.Given the inapplicability of the RF model, we evaluate both the Generalized Gauss\-Newton \(GGN\) approximation and the exact directional landscape, which qualitatively agree \([Figure˜I\.5](https://arxiv.org/html/2605.11181#A9.F5)\(ab\)\)\. For a single batch, evaluating the per\-step decrease confirms thatFreonatc=1c=1yields the optimal local decrease \([Figure˜I\.5](https://arxiv.org/html/2605.11181#A9.F5)\(c\)\), consistently across all layers \([Figure˜I\.11](https://arxiv.org/html/2605.11181#A9.F11)\)\. This shows that the local geometry is severely distorted: favoringc=1c=1is impossible in standard RF models \([Figure˜I\.3](https://arxiv.org/html/2605.11181#A9.F3)\)\. Because local geometry poorly predicts overall performance, we evaluate the global dynamics across the full dataset, revealing a highly stochastic landscape: rather than converging, the optimal exponents fluctuate dynamically over batches during training \([Figure˜6](https://arxiv.org/html/2605.11181#S3.F6),[Figure˜I\.7](https://arxiv.org/html/2605.11181#A9.F7)\)\. Thus, the actual global landscape is significantly more complex than local, single\-batch models suggest\.

##### Φ\\Phiandγ\\gammadynamics in GPT\-2:

TrackingΦ\\Phiandγ\\gammathroughout training reveals a shifting optimization hierarchy \([Figure˜I\.1](https://arxiv.org/html/2605.11181#A9.F1)\)\. Initially, methods likeTruncatedSGDexhibit strong local decrease \(highΦ\\Phi\), with mostFreonvariants similar\. At the next snapshot, this inverts:TruncatedSGDandSGDbecome heavily suboptimal with larger variance\. This mirrors the validation loss trajectories[Figure˜I\.2](https://arxiv.org/html/2605.11181#A9.F2), where the rate of decrease slows down forSGDvariants, but remains good forFreon\.

This reveals an interesting empirical observation: gradients tend to agree most in the middle of the singular value spectrum \(in the sense of signal\-to\-noise ratio\), rather than at the extremes \(c=0c=0orc=1c=1\)\. Increasingccto target this high\-SNR mid\-range \(e\.g\.,c≈3/4c\\approx 3/4\) naturally suppresses the mean and variance ofγ\\gamma, as updates align with these mid\-tier rather than maximal singular values\.

Overall, these empirical observations clarify the mechanics of[Equation˜4](https://arxiv.org/html/2605.11181#S3.E4): changing directions of descent from the gradient naturally suppressesγ\\gamma, yet this deliberate sacrifice ofγ\\gammaenables a massive amplification ofΦ\\Phi– the primary driver of loss reduction\. Tradingγ\\gammafor a boostedΦ\\Phifundamentally explains why the non\-standardc\>1/2c\>1/2regime ultimately dominates the optimization landscape\.

## 4Conclusion

This work challenges the prevailing geometric narrative of spectral optimizers likeMuon\.Freonreveals optimal updates often violate proper norm regimes, whileKaonproves randomized spectra perform just as well\. What actually drives performance is the interplay between batch gradient alignment and directional descent potential, combined with appropriate choices of step\-sizes\. We discuss limitations in[Appendix˜C](https://arxiv.org/html/2605.11181#A3), leaving a final question: if geometry is not necessary for performance, should we expect theories stemming from it to yield practically useful insights?

## References

- Handbook of mathematical functions with formulas, graphs, and mathematical tables\.National Bureau of Standards / Dover\.Note:Section 15\.1Cited by:[§E\.1\.3](https://arxiv.org/html/2605.11181#A5.SS1.SSS3.p2.2)\.
- N\. I\. Achieser \(1992\)Theory of approximation\.Dover Publications, Inc\., New York\.Note:Translated from the Russian and with a preface by Charles J\. Hyman, Reprint of the 1956 English translationExternal Links:ISBN 0\-486\-67129\-1,[MathReview Entry](https://www.ams.org/mathscinet-getitem?mr=1217081)Cited by:[§E\.1\.1](https://arxiv.org/html/2605.11181#A5.SS1.SSS1.1.p1.10),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- K\. Ahn, B\. Xu, N\. Abreu, Y\. Fan, G\. Magakyan, P\. Sharma, Z\. Zhan, and J\. Langford \(2025\)Dion: distributed orthonormalized updates\.arXiv preprint arXiv:2504\.05295\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- N\. Amsel, D\. Persson, C\. Musco, and R\. M\. Gower \(2026\)The polar express: optimal matrix sign methods and their application to the muon algorithm\.InThe Fourteenth International Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=yRtgZ1K8hO)Cited by:[Appendix B](https://arxiv.org/html/2605.11181#A2.p1.3),[§E\.1\.2](https://arxiv.org/html/2605.11181#A5.SS1.SSS2.4.p2.3),[§E\.1\.2](https://arxiv.org/html/2605.11181#A5.SS1.SSS2.p2.1),[§E\.1](https://arxiv.org/html/2605.11181#A5.SS1.p1.2),[item 1](https://arxiv.org/html/2605.11181#S1.I1.i1.p1.5),[§1](https://arxiv.org/html/2605.11181#S1.p2.1),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px1.p1.6),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- \[5\]J\. Bernstein and L\. NewhouseOld optimizer, new norm: an anthology\.InOPT 2024: Optimization for Machine Learning,Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- R\. Bhatia \(2013\)Matrix analysis\.Springer Science & Business Media\.Cited by:[§2\.4](https://arxiv.org/html/2605.11181#S2.SS4.p1.8)\.
- E\. W\. Cheney and H\. L\. Loeb \(1964\)Generalized rational approximation\.SIAM Journal on Numerical Analysis, Ser\. B1\(1\),pp\. 11–25\.Cited by:[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- J\. M\. Cohen, S\. Kaur, Y\. Li, J\. Z\. Kolter, and A\. Talwalkar \(2021\)Gradient descent on neural networks typically occurs at the edge of stability\.arXiv preprint arXiv:2103\.00065\.Cited by:[§3](https://arxiv.org/html/2605.11181#S3.p1.1)\.
- R\. Couillet and Z\. Liao \(2022\)Random matrix theory\.InRandom Matrix Methods for Machine Learning,pp\. 35–154\.Cited by:[Appendix G](https://arxiv.org/html/2605.11181#A7.2.p2.7)\.
- M\. Crawshaw, C\. Modi, M\. Liu, and R\. M\. Gower \(2025\)An exploration of non\-euclidean gradient descent: muon and its many variants\.External Links:2510\.09827,[Link](https://arxiv.org/abs/2510.09827)Cited by:[Appendix B](https://arxiv.org/html/2605.11181#A2.p1.3)\.
- D\. Davis and D\. Drusvyatskiy \(2026\)When do spectral gradient updates help in deep learning?\.External Links:2512\.04299,[Link](https://arxiv.org/abs/2512.04299)Cited by:[§H\.1](https://arxiv.org/html/2605.11181#A8.SS1.p1.10),[item 3](https://arxiv.org/html/2605.11181#S1.I1.i3.p1.3),[Figure 5](https://arxiv.org/html/2605.11181#S3.F5),[Figure 5](https://arxiv.org/html/2605.11181#S3.F5.4.2.3),[item 2](https://arxiv.org/html/2605.11181#S3.I1.i2.p1.5),[§3\.2](https://arxiv.org/html/2605.11181#S3.SS2.p1.2),[§3\.3](https://arxiv.org/html/2605.11181#S3.SS3.SSS0.Px1.p1.3),[§3\.3](https://arxiv.org/html/2605.11181#S3.SS3.SSS0.Px2.p1.2),[§3](https://arxiv.org/html/2605.11181#S3.p1.1)\.
- DeepSeek\-AI \(2026\)DeepSeek\-v4: towards highly efficient million\-token context intelligence\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- Z\. Du and W\. Su \(2026\)The newton\-muon optimizer\.External Links:2604\.01472,[Link](https://arxiv.org/abs/2604.01472)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- J\. Duchi, E\. Hazan, and Y\. Singer \(2011\)Adaptive subgradient methods for online learning and stochastic optimization\.\.Journal of machine learning research12\(7\)\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1)\.
- C\. Fan, M\. Schmidt, and C\. Thrampoulidis \(2025\)Implicit bias of spectral descent and muon on multiclass separable data\.arXiv preprint arXiv:2502\.04664\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- S\. Filip, Y\. Nakatsukasa, L\. N\. Trefethen, and B\. Beckermann \(2018\)Rational minimax approximation via adaptive barycentric representations\.SIAM Journal on Scientific Computing40\(4\),pp\. A2427–A2455\.External Links:[Document](https://dx.doi.org/10.1137/17M1132409),[Link](https://doi.org/10.1137/17M1132409),https://doi\.org/10\.1137/17M1132409Cited by:[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- E\. S\. Gawlik and Y\. Nakatsukasa \(2021\)Approximating the pth root by composite rational functions\.Journal of Approximation Theory266,pp\. 105577\.Cited by:[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p2.3)\.
- W\. Gong, J\. Zazo, Q\. Luo, P\. Wang, J\. Hensman, and C\. Ma \(2026\)ARO: a new lens on matrix optimization for large models\.External Links:2602\.09006,[Link](https://arxiv.org/abs/2602.09006)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- A\. Gonon, A\. Muşat, and N\. Boumal \(2026\)Insights on muon from simple quadratics\.External Links:2602\.11948,[Link](https://arxiv.org/abs/2602.11948)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p3.1),[§3\.3](https://arxiv.org/html/2605.11181#S3.SS3.SSS0.Px1.p1.3)\.
- W\. B\. Gragg \(1972\)The padé table and its relation to certain algorithms of numerical analysis\.SIAM Review14\(1\),pp\. 1–62\.External Links:ISSN 00361445, 10957200,[Link](http://www.jstor.org/stable/2028911)Cited by:[§E\.1\.2](https://arxiv.org/html/2605.11181#A5.SS1.SSS2.3.p1.5)\.
- V\. Gupta, T\. Koren, and Y\. Singer \(2018\)Shampoo: preconditioned stochastic tensor optimization\.InInternational Conference on Machine Learning,pp\. 1842–1850\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1)\.
- N\. J\. Higham \(2008\)Functions of matrices: theory and computation\.SIAM\.Cited by:[§E\.4\.1](https://arxiv.org/html/2605.11181#A5.SS4.SSS1.Px4.p1.3),[§E\.4\.2](https://arxiv.org/html/2605.11181#A5.SS4.SSS2.Px1.p1.1),[§E\.4\.3](https://arxiv.org/html/2605.11181#A5.SS4.SSS3.Px1.p1.1),[§E\.4\.5](https://arxiv.org/html/2605.11181#A5.SS4.SSS5.Px1.p1.1),[§E\.4](https://arxiv.org/html/2605.11181#A5.SS4.p2.1),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px1.p1.6)\.
- R\. Islamov, M\. Crawshaw, J\. Cohen, and R\. Gower \(2026\)Non\-euclidean gradient descent operates at the edge of stability\.arXiv preprint arXiv:2603\.05002\.Cited by:[§3](https://arxiv.org/html/2605.11181#S3.p1.1)\.
- X\. Jiang, A\. Semenov, and S\. U\. Stich \(2026\)Enhancing llm training via spectral clipping\.External Links:2603\.14315,[Link](https://arxiv.org/abs/2603.14315)Cited by:[§2\.2](https://arxiv.org/html/2605.11181#S2.SS2.p2.1)\.
- K\. Jordan, Y\. Jin, V\. Boza, J\. You, F\. Cesista, L\. Newhouse, and J\. Bernstein \(2024\)Muon: an optimizer for hidden layers in neural networks\.External Links:[Link](https://kellerjordan.github.io/posts/muon/)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px1.p1.6)\.
- C\. Kenney and A\. J\. Laub \(1991\)Rational iterative methods for the matrix sign function\.SIAM Journal on Matrix Analysis and Applications12\(2\),pp\. 273–291\.External Links:[Document](https://dx.doi.org/10.1137/0612020),[Link](https://doi.org/10.1137/0612020),https://doi\.org/10\.1137/0612020Cited by:[§E\.1\.3](https://arxiv.org/html/2605.11181#A5.SS1.SSS3.p2.3)\.
- J\. Kim, E\. Nichani, D\. Wu, A\. Bietti, and J\. D\. Lee \(2026\)Sharp capacity scaling of spectral optimizers in learning associative memory\.External Links:2603\.26554,[Link](https://arxiv.org/abs/2603.26554)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- D\. P\. Kingma and J\. Ba \(2017\)Adam: a method for stochastic optimization\.External Links:1412\.6980,[Link](https://arxiv.org/abs/1412.6980)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1)\.
- D\. Kovalev \(2025\)Understanding gradient orthogonalization for deep learning via non\-euclidean trust\-region optimization\.arXiv preprint arXiv:2503\.12645\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- F\. Kunstner, P\. Hennig, and L\. Balles \(2019\)Limitations of the empirical Fisher approximation for natural gradient descent\.InAdvances in Neural Information Processing Systems,Vol\.32\.External Links:[Link](https://papers.neurips.cc/paper_files/paper/2019/hash/46a558d97954d0692411c861cf78ef79-Abstract.html)Cited by:[§D\.2](https://arxiv.org/html/2605.11181#A4.SS2.SSS0.Px1.p5.6)\.
- J\. Liu, J\. Su, X\. Yao, Z\. Jiang, G\. Lai, Y\. Du, Y\. Qin, W\. Xu, E\. Lu, J\. Yan, Y\. Chen, H\. Zheng, Y\. Liu, S\. Liu, B\. Yin, W\. He, H\. Zhu, Y\. Wang, J\. Wang, M\. Dong, Z\. Zhang, Y\. Kang, H\. Zhang, X\. Xu, Y\. Zhang, Y\. Wu, X\. Zhou, and Z\. Yang \(2025\)Muon is scalable for llm training\.External Links:2502\.16982,[Link](https://arxiv.org/abs/2502.16982)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- H\. L\. Loeb \(1966\)Approximation by generalized rationals\.SIAM Journal on Numerical Analysis3\(1\),pp\. 34–55\.External Links:[Document](https://dx.doi.org/10.1137/0703003)Cited by:[§E\.1\.1](https://arxiv.org/html/2605.11181#A5.SS1.SSS1.3.p1.4)\.
- I\. Loshchilov and F\. Hutter \(2019\)Decoupled weight decay regularization\.External Links:1711\.05101,[Link](https://arxiv.org/abs/1711.05101)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1)\.
- J\. Ma, Y\. Huang, Y\. Chi, and Y\. Chen \(2026\)Preconditioning benefits of spectral orthogonalization in muon\.External Links:2601\.13474,[Link](https://arxiv.org/abs/2601.13474)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- S\. Merity, C\. Xiong, J\. Bradbury, and R\. Socher \(2016\)The WikiText long term dependency language modeling dataset\.Cited by:[Appendix B](https://arxiv.org/html/2605.11181#A2.p1.3)\.
- Y\. Nakatsukasa, Z\. Bai, and F\. Gygi \(2010\)Optimizing halley’s iteration for computing the matrix polar decomposition\.SIAM Journal on Matrix Analysis and Applications31\(5\),pp\. 2700–2720\.Cited by:[item 1](https://arxiv.org/html/2605.11181#S1.I1.i1.p1.5),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- Y\. Nakatsukasa and R\. W\. Freund \(2016\)Computing fundamental matrix decompositions accurately via the matrix sign function in two iterations: the power of zolotarev’s functions\.SIAM Review58\(3\),pp\. 461–493\.External Links:[Document](https://dx.doi.org/10.1137/140990334),[Link](https://doi.org/10.1137/140990334),https://doi\.org/10\.1137/140990334Cited by:[§E\.4\.7](https://arxiv.org/html/2605.11181#A5.SS4.SSS7.Px1.p1.1),[item 1](https://arxiv.org/html/2605.11181#S1.I1.i1.p1.5),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px2.p1.1)\.
- T\. Pethick, W\. Xie, K\. Antonakopoulos, Z\. Zhu, A\. Silveti\-Falls, and V\. Cevher \(2025\)Training deep learning models with norm\-constrained lmos\.arXiv preprint arXiv:2502\.07529\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- S\. C\. Phatak and S\. S\. Rao \(1995\)Logistic map: a possible random\-number generator\.Physical review E51\(4\),pp\. 3670\.Cited by:[§2\.3](https://arxiv.org/html/2605.11181#S2.SS3.p2.2)\.
- \[40\]A\. Riabinin, E\. Shulgin, K\. Gruntkowska, and P\. RichtárikGluon: making muon & scion great again\!\(bridging theory and practice of lmo\-based optimizers for llms\)\.InHigh\-dimensional Learning Dynamics 2025,Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- A\. Semenov, M\. Pagliardini, and M\. Jaggi \(2025\)Benchmarking optimizers for large language model pretraining\.External Links:2509\.01440,[Link](https://arxiv.org/abs/2509.01440)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p3.1)\.
- C\. Si, D\. Zhang, and W\. Shen \(2025\)Adamuon: adaptive muon optimizer\.arXiv preprint arXiv:2507\.11005\.Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- W\. Su \(2025\)Isotropic curvature model for understanding deep learning optimization: is gradient orthogonalization optimal?\.External Links:2511\.00674,[Link](https://arxiv.org/abs/2511.00674)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p3.1)\.
- 5\. Team, A\. Zeng, X\. Lv, Q\. Zheng, Z\. Hou, B\. Chen, C\. Xie, C\. Wang, D\. Yin, H\. Zeng, J\. Zhang, K\. Wang, L\. Zhong, M\. Liu, R\. Lu, S\. Cao, X\. Zhang, X\. Huang, Y\. Wei, Y\. Cheng, Y\. An, Y\. Niu, Y\. Wen, Y\. Bai, Z\. Du, Z\. Wang, Z\. Zhu, B\. Zhang, B\. Wen, B\. Wu, B\. Xu, C\. Huang, C\. Zhao, C\. Cai, C\. Yu, C\. Li, C\. Ge, C\. Huang, C\. Zhang, C\. Xu, C\. Zhu, C\. Li, C\. Yin, D\. Lin, D\. Yang, D\. Jiang, D\. Ai, E\. Zhu, F\. Wang, G\. Pan, G\. Wang, H\. Sun, H\. Li, H\. Li, H\. Hu, H\. Zhang, H\. Peng, H\. Tai, H\. Zhang, H\. Wang, H\. Yang, H\. Liu, H\. Zhao, H\. Liu, H\. Yan, H\. Liu, H\. Chen, J\. Li, J\. Zhao, J\. Ren, J\. Jiao, J\. Zhao, J\. Yan, J\. Wang, J\. Gui, J\. Zhao, J\. Liu, J\. Li, J\. Li, J\. Lu, J\. Wang, J\. Yuan, J\. Li, J\. Du, J\. Du, J\. Liu, J\. Zhi, J\. Gao, K\. Wang, L\. Yang, L\. Xu, L\. Fan, L\. Wu, L\. Ding, L\. Wang, M\. Zhang, M\. Li, M\. Xu, M\. Zhao, M\. Zhai, P\. Du, Q\. Dong, S\. Lei, S\. Tu, S\. Yang, S\. Lu, S\. Li, S\. Li, Shuang\-Li, S\. Yang, S\. Yi, T\. Yu, W\. Tian, W\. Wang, W\. Yu, W\. L\. Tam, W\. Liang, W\. Liu, X\. Wang, X\. Jia, X\. Gu, X\. Ling, X\. Wang, X\. Fan, X\. Pan, X\. Zhang, X\. Zhang, X\. Fu, X\. Zhang, Y\. Xu, Y\. Wu, Y\. Lu, Y\. Wang, Y\. Zhou, Y\. Pan, Y\. Zhang, Y\. Wang, Y\. Li, Y\. Su, Y\. Geng, Y\. Zhu, Y\. Yang, Y\. Li, Y\. Wu, Y\. Li, Y\. Liu, Y\. Wang, Y\. Li, Y\. Zhang, Z\. Liu, Z\. Yang, Z\. Zhou, Z\. Qiao, Z\. Feng, Z\. Liu, Z\. Zhang, Z\. Wang, Z\. Yao, Z\. Wang, Z\. Liu, Z\. Chai, Z\. Li, Z\. Zhao, W\. Chen, J\. Zhai, B\. Xu, M\. Huang, H\. Wang, J\. Li, Y\. Dong, and J\. Tang \(2025a\)GLM\-4\.5: agentic, reasoning, and coding \(arc\) foundation models\.External Links:2508\.06471,[Link](https://arxiv.org/abs/2508.06471)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- K\. Team, Y\. Bai, Y\. Bao, G\. Chen, J\. Chen, N\. Chen, R\. Chen, Y\. Chen, Y\. Chen, Y\. Chen, Z\. Chen, J\. Cui, H\. Ding, M\. Dong, A\. Du, C\. Du, D\. Du, Y\. Du, Y\. Fan, Y\. Feng, K\. Fu, B\. Gao, H\. Gao, P\. Gao, T\. Gao, X\. Gu, L\. Guan, H\. Guo, J\. Guo, H\. Hu, X\. Hao, T\. He, W\. He, W\. He, C\. Hong, Y\. Hu, Z\. Hu, W\. Huang, Z\. Huang, Z\. Huang, T\. Jiang, Z\. Jiang, X\. Jin, Y\. Kang, G\. Lai, C\. Li, F\. Li, H\. Li, M\. Li, W\. Li, Y\. Li, Y\. Li, Z\. Li, Z\. Li, H\. Lin, X\. Lin, Z\. Lin, C\. Liu, C\. Liu, H\. Liu, J\. Liu, J\. Liu, L\. Liu, S\. Liu, T\. Y\. Liu, T\. Liu, W\. Liu, Y\. Liu, Y\. Liu, Y\. Liu, Y\. Liu, Z\. Liu, E\. Lu, L\. Lu, S\. Ma, X\. Ma, Y\. Ma, S\. Mao, J\. Mei, X\. Men, Y\. Miao, S\. Pan, Y\. Peng, R\. Qin, B\. Qu, Z\. Shang, L\. Shi, S\. Shi, F\. Song, J\. Su, Z\. Su, X\. Sun, F\. Sung, H\. Tang, J\. Tao, Q\. Teng, C\. Wang, D\. Wang, F\. Wang, H\. Wang, J\. Wang, J\. Wang, J\. Wang, S\. Wang, S\. Wang, Y\. Wang, Y\. Wang, Y\. Wang, Y\. Wang, Y\. Wang, Z\. Wang, Z\. Wang, Z\. Wang, C\. Wei, Q\. Wei, W\. Wu, X\. Wu, Y\. Wu, C\. Xiao, X\. Xie, W\. Xiong, B\. Xu, J\. Xu, J\. Xu, L\. H\. Xu, L\. Xu, S\. Xu, W\. Xu, X\. Xu, Y\. Xu, Z\. Xu, J\. Yan, Y\. Yan, X\. Yang, Y\. Yang, Z\. Yang, Z\. Yang, Z\. Yang, H\. Yao, X\. Yao, W\. Ye, Z\. Ye, B\. Yin, L\. Yu, E\. Yuan, H\. Yuan, M\. Yuan, H\. Zhan, D\. Zhang, H\. Zhang, W\. Zhang, X\. Zhang, Y\. Zhang, Y\. Zhang, Y\. Zhang, Y\. Zhang, Y\. Zhang, Y\. Zhang, Z\. Zhang, H\. Zhao, Y\. Zhao, H\. Zheng, S\. Zheng, J\. Zhou, X\. Zhou, Z\. Zhou, Z\. Zhu, W\. Zhuang, and X\. Zu \(2025b\)Kimi k2: open agentic intelligence\.External Links:2507\.20534,[Link](https://arxiv.org/abs/2507.20534)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- \[46\]N\. Vyas, D\. Morwani, R\. Zhao, I\. Shapira, D\. Brandfonbrener, L\. Janson, and S\. M\. KakadeSOAP: improving and stabilizing shampoo using adam for language modeling\.InThe Thirteenth International Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p1.1)\.
- H\. S\. Wall \(1948\)Analytic theory of continued fractions\.D\. Van Nostrand Company, Inc\.,New York\.Note:Reprinted by Chelsea Publishing Company, Bronx, N\.Y\., 1967Cited by:[§E\.1\.2](https://arxiv.org/html/2605.11181#A5.SS1.SSS2.3.p1.5)\.
- S\. Wang, F\. Zhang, J\. Li, C\. Du, C\. Du, T\. Pang, Z\. Yang, M\. Hong, and V\. Y\. F\. Tan \(2025\)Muon outperforms adam in tail\-end associative memory learning\.External Links:2509\.26030,[Link](https://arxiv.org/abs/2509.26030)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- K\. Wen, D\. Hall, T\. Ma, and P\. Liang \(2025\)Fantastic pretraining optimizers and where to find them\.External Links:2509\.02046,[Link](https://arxiv.org/abs/2509.02046)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p3.1)\.
- R\. Xu, J\. Li, and Y\. Lu \(2026\)On the width scaling of neural optimizers under matrix operator norms i: row/column normalization and hyperparameter transfer\.External Links:2603\.09952,[Link](https://arxiv.org/abs/2603.09952)Cited by:[§2\.4\.1](https://arxiv.org/html/2605.11181#S2.SS4.SSS1.p1.18)\.
- J\. Yen, S\. Si, Z\. Meng, F\. Yu, S\. S\. Duvvuri, I\. S\. Dhillon, C\. Hsieh, and S\. Kumar \(2025\)LoRA done RITE: robust invariant transformation equilibration for loRA optimization\.InThe Thirteenth International Conference on Learning Representations,External Links:[Link](https://openreview.net/forum?id=VpWki1v2P8)Cited by:[item 1](https://arxiv.org/html/2605.11181#S1.I1.i1.p1.5)\.
- D\. Yu, H\. Tao, Y\. Wan, L\. Luo, and L\. Zhang \(2026\)Sign\-based optimizers are effective under heavy\-tailed noise\.External Links:2602\.07425,[Link](https://arxiv.org/abs/2602.07425)Cited by:[§1](https://arxiv.org/html/2605.11181#S1.p2.1)\.
- J\. Zhang, N\. Amsel, B\. Chen, and T\. Dao \(2026\)Gram newton\-schulz\.External Links:[Link](https://dao-ailab.github.io/blog/2026/gram-newton-schulz/)Cited by:[§E\.4\.2](https://arxiv.org/html/2605.11181#A5.SS4.SSS2.Px1.p1.1),[§2\.4\.3](https://arxiv.org/html/2605.11181#S2.SS4.SSS3.Px1.p1.6)\.

## Appendix ANanoGPT Tuning Details for[Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[3](https://arxiv.org/html/2605.11181#S1.F3)

The NanoGPT results in[Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[3](https://arxiv.org/html/2605.11181#S1.F3)were run through the same training pipeline on NanoGPT\. Every tuning run and learning rate sweep used 8 NVIDIA H100 80GB GPUs\.

For each optimizer, we tuned a matrix\-update learning rateηmatrix\\eta\_\{\\mathrm\{matrix\}\}and a base learning rateηbase\\eta\_\{\\mathrm\{base\}\}\. The recorded tuning grid was

\(5e\-3, 1e\-3\)\(1e\-2, 1e\-3\)\(2e\-2, 1e\-3\)\(1e\-2, 5e\-4\)\(5e\-3, 2e\-3\)\(5e\-3, 4e\-3\)\(5e\-3, 8e\-3\)\(5e\-3, 2e\-2\)\(1e\-2, 2e\-3\)\(1e\-2, 4e\-3\)\(1e\-2, 8e\-3\)\(1e\-2, 2e\-2\)\(2e\-2, 2e\-3\)\(2e\-2, 4e\-3\)\(2e\-2, 8e\-3\)\(2e\-2, 2e\-2\)for\(matrix\_lr,base\_lr\)\(\\texttt\{matrix\\\_lr\},\\texttt\{base\\\_lr\}\)\. The selected learning rates used for[Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[3](https://arxiv.org/html/2605.11181#S1.F3)are listed in[Table˜A\.1](https://arxiv.org/html/2605.11181#A1.T1)\.Kaonand bothFreonvariants used 5 Newton–Schulz/QDWH\-style matrix iterations\.

Table A\.1:Learning rates selected by the seed\-42 tuning sweep and used for[Figures˜3](https://arxiv.org/html/2605.11181#S1.F3)and[3](https://arxiv.org/html/2605.11181#S1.F3)\.[Figure˜3](https://arxiv.org/html/2605.11181#S1.F3)then jointly scales the tuned pair as\(ηmatrix,ηbase\)=ρ​\(ηmatrix⋆,ηbase⋆\)\(\\eta\_\{\\mathrm\{matrix\}\},\\eta\_\{\\mathrm\{base\}\}\)=\\rho\(\\eta\_\{\\mathrm\{matrix\}\}^\{\\star\},\\eta\_\{\\mathrm\{base\}\}^\{\\star\}\)\. The full sensitivity sweep ran 4 optimizers, 7 scale factors, and 3 seeds, for 84 runs in total\. The scale factors wereρ∈\{0\.03,0\.1,0\.3,1,3,10,30,100\}\\rho\\in\\\{0\.03,0\.1,0\.3,1,3,10,30,100\\\}and the seeds were42,43,4442,43,44\. The plotted version omitsρ=100\\rho=100, so the x\-axis in[Figure˜3](https://arxiv.org/html/2605.11181#S1.F3)shows the matrix\-update learning rates induced byρ∈\{0\.03,0\.1,0\.3,1,3,10,30\}\\rho\\in\\\{0\.03,0\.1,0\.3,1,3,10,30\\\}\. Error bars show±2\\pm 2std over the three seeds\.

[Figure˜3](https://arxiv.org/html/2605.11181#S1.F3)reruns the selected curve configurations with denser validation logging\. These dense\-validation jobs used the same tuned learning rates and the same seeds42,43,4442,43,44, with validation every 8 optimizer steps\. The dense runs generated Muon atρ=1\\rho=1,Kaonatρ=1\\rho=1andρ=3\\rho=3,Freonc=2/3c=2/3atρ=1\\rho=1, andFreonc=3/4c=3/4atρ=1\\rho=1\. The plotted dense\-validation figure uses Muon atρ=1\\rho=1,Kaonatρ=3\\rho=3, and bothFreonvariants atρ=1\\rho=1\. The plotted lines show averages across three seeds and shaded regions showing±2\\pm 2std\.

## Appendix BGPT2 Training Details

We utilize the existing codebase from\[Amselet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib12), Crawshawet al\.,[2025](https://arxiv.org/html/2605.11181#bib.bib2)\]\. All experiments train a GPT\-2 model \(12 layers, 12 attention heads, embedding dimension 768, 124M parameters\) with GELU activations and FlashAttention on the WikiText\-2 dataset\[Merityet al\.,[2016](https://arxiv.org/html/2605.11181#bib.bib1)\]\. Each run uses a context length of512512tokens, batch size88, gradient norm clipping at1\.01\.0, and a constant\-then\-linear decay schedule with a 5% warm\-up\. The combined figure[Figure˜4](https://arxiv.org/html/2605.11181#S2.F4)reports final validation loss as a function of learning rate, sweeping 8–10 learning rates per optimizer configuration\.

SGD and TruncatedSGD\.Runs use SGD with momentum0\.90\.9\. Truncation levels of0,0\.1,0\.5,1,2,5,100,0\.1,0\.5,1,2,5,10% are swept, where0% recovers vanilla SGD\. Runs with non\-zero truncation rescale the update 2\-norm to match that of the untruncated gradient\.

Muon\.Standard Muon with momentum0\.950\.95and 5 Newton–Schulz iterations \(p/q=1/2p/q=1/2\)\.

Kaon\.Same Newton–Schulz backbone as Muon \(55steps, momentum0\.950\.95\) but with the chaos iterations\.

Freon\.FreonRational with momentum0\.950\.95, 5 Newton–Schulz steps, and the exponentc=a/bc=a/bswept over0,1/4,1/3,1/2,2/3,3/4,1\{0,1/4,1/3,1/2,2/3,3/4,1\}, wherec=0c=0recovers gradient descent andc=1/2c=1/2recovers Muon\.

## Appendix CLimitations

While our analysis provides novel mechanistic insights into the behavior of spectral optimizers, there is a number of theoretical and practical limitations:

1. 1\.Scope of Empirics and Theory:Our empirical validation is currently restricted to a random feature model and a single architecture modality \(a GPT\-2 language model\)\. Furthermore, our rigorous analysis of the two core optimization quantities, batch gradient alignment and directional descent potential, is limited to the RF setting\. It remains an open question how precisely these dynamics map to other modalities \(e\.g\., vision\) or more theoretically amenable non\-quadratic loss landscapes\.
2. 2\.Theoretical Assumptions:The exact asymptotic results of[Section˜3\.2](https://arxiv.org/html/2605.11181#S3.SS2)rely strictly on the assumption of mean\-zero activations, excluding activations such as ReLU\.
3. 3\.Lack of Predictive Power:Our framework primarily serves as a diagnostic lens\. The quantities governing the optimization trajectory are highly stochastic, meaning they can only be reliably evaluated post\-factum rather than predicted a priori\.
4. 4\.Optimalcc:Because the optimal local geometry is heavily distorted and highly sensitive to individual batch dynamics, our attempts to dynamically track and utilize the optimal Schatten parameters \(p,qp,q\) step\-by\-step ultimately failed\. The optimal exponents fluctuate too wildly across batches to formulate a stable, greedily adaptive schedule\.
5. 5\.QDWH Coupled CholeskyThe current implementation of theFreonalgorithm serves primarily to demonstrate the theoretical viability of quasi\-norm updates\. While the algorithm could certainly be engineered for greater computational efficiency, it is currently unclear whether the extensive engineering effort is even necessary forFreon, given the similar behaviour across all the optimizers\.

## Appendix DMethods Proofs

### D\.1Preconditioned Spectral Descent

See[2\.1](https://arxiv.org/html/2605.11181#S2.Thmtheorem1)

###### Proof\.

Let∥⋅∥\\\|\\cdot\\\|be a unitarily invariant matrix norm and∥⋅∥∗\\\|\\cdot\\\|\_\{\*\}its dual norm\. Write the SVDGk=Uk​diag⁡\(σk\)​Vk⊤G\_\{k\}=U\_\{k\}\\operatorname\{diag\}\(\\sigma\_\{k\}\)V\_\{k\}^\{\\top\}\. Then since∥⋅∥\\\|\\cdot\\\|is unitarily invariant, there is a norm∥⋅∥σ\\\|\\cdot\\\|\_\{\\sigma\}onℝr\\mathbb\{R\}^\{r\}symmetric in the coordinates such that‖Gk‖=‖σk‖σ\\\|G\_\{k\}\\\|=\\\|\\sigma\_\{k\}\\\|\_\{\\sigma\}\. The first part of the theorem then follows from formulating steepest descent as

Xk\+1\\displaystyle X\_\{k\+1\}=Xk−η​‖Gk‖∗​lmo∥⋅∥⁡\(Gk\)\\displaystyle=\{X\}\_\{k\}\-\\eta\\left\\\|\{G\}\_\{k\}\\right\\\|\_\{\*\}\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\\left\(\{G\}\_\{k\}\\right\)=Xk−αk​⟨Gk,Dk⟩​Dk\\displaystyle=\{X\}\_\{k\}\-\\alpha\_\{k\}\\langle G\_\{k\},D\_\{k\}\\rangle D\_\{k\}withαk=η\\alpha\_\{k\}=\\etaandDk=lmo∥⋅∥⁡\(Gk\)=Uk​lmo∥⋅∥σ⁡\(σk\)​Vk⊤D\_\{k\}=\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\}\(G\_\{k\}\)=U\_\{k\}\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}\(\\sigma\_\{k\}\)V\_\{k\}^\{\\top\}\.

Now consider a preconditioned spectral descent update\. Then if it is a steepest descent update, we have as in the above argument that

pk​\(σk\)=η​‖Gk‖∗αk​⟨Gk,Dk⟩​lmo∥⋅∥σ⁡\(σk\)∝lmo∥⋅∥σ⁡\(σk\)\.p\_\{k\}\(\\sigma\_\{k\}\)=\\frac\{\\eta\\\|G\_\{k\}\\\|\_\{\*\}\}\{\\alpha\_\{k\}\\langle G\_\{k\},D\_\{k\}\\rangle\}\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}\(\\sigma\_\{k\}\)\\propto\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}\(\\sigma\_\{k\}\)\.\(6\)
WriteTk=lmo∥⋅∥σ⁡\(σk\)T\_\{k\}=\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}\(\\sigma\_\{k\}\)\. Suppose\[σk\]i≤\[σk\]j\[\\sigma\_\{k\}\]\_\{i\}\\leq\[\\sigma\_\{k\}\]\_\{j\}and\[Tk\]i\>\[Tk\]j\[T\_\{k\}\]\_\{i\}\>\[T\_\{k\}\]\_\{j\}\. DefineT~k∈ℝr\\tilde\{T\}\_\{k\}\\in\\mathbb\{R\}^\{r\}such that\[T~k\]i=\[Tk\]j\[\\tilde\{T\}\_\{k\}\]\_\{i\}=\[T\_\{k\}\]\_\{j\},\[T~k\]j=\[Tk\]i\[\\tilde\{T\}\_\{k\}\]\_\{j\}=\[T\_\{k\}\]\_\{i\}and\[T~k\]i~=\[Tk\]i~\[\\tilde\{T\}\_\{k\}\]\_\{\\tilde\{i\}\}=\[T\_\{k\}\]\_\{\\tilde\{i\}\}for alli~≠i,j\\tilde\{i\}\\neq i,j\. Then since∥⋅∥σ\\\|\\cdot\\\|\_\{\\sigma\}is symmetric in the coordinates,‖T~k‖σ=‖Tk‖σ=1\\\|\\tilde\{T\}\_\{k\}\\\|\_\{\\sigma\}=\\\|T\_\{k\}\\\|\_\{\\sigma\}=1\. But we see that⟨σk,T~k⟩\>⟨σk,Tk⟩\\langle\\sigma\_\{k\},\\tilde\{T\}\_\{k\}\\rangle\>\\langle\\sigma\_\{k\},T\_\{k\}\\rangle\. This contradictsTk=lmo∥⋅∥σ⁡\(σk\)T\_\{k\}=\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}\(\\sigma\_\{k\}\)\. Thuslmo∥⋅∥σ\\operatorname\{lmo\}\_\{\\\|\\cdot\\\|\_\{\\sigma\}\}preserves the order≤\\leqof the entries ofσk\\sigma\_\{k\}, hence so doespkp\_\{k\}\.

∎

See[2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)

###### Proof\.

AsffhasL𝒳L\_\{\\mathcal\{X\}\}\-Lipschitz continuous gradients with respect to∥⋅∥\\\|\\cdot\\\|, we can use the generalized descent lemma, which for the total updateΔk=⟨Gk,Dk⟩​Dk\\Delta\_\{k\}=\\langle G\_\{k\},D\_\{k\}\\rangle D\_\{k\}yields:

f​\(Xk\+1\)\\displaystyle f\(X\_\{k\+1\}\)≤f​\(Xk\)−αk​⟨Gk,Δk⟩\+L𝒳2​αk2​‖Δk‖2\\displaystyle\\leq f\(X\_\{k\}\)\-\\alpha\_\{k\}\\langle G\_\{k\},\\Delta\_\{k\}\\rangle\+\\frac\{L\_\{\\mathcal\{X\}\}\}\{2\}\\alpha\_\{k\}^\{2\}\\\|\\Delta\_\{k\}\\\|^\{2\}=f​\(Xk\)−αk​⟨Gk,Dk⟩2\+L𝒳2​αk2​⟨Gk,Dk⟩2​‖Dk‖2\\displaystyle=f\(X\_\{k\}\)\-\\alpha\_\{k\}\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}\+\\frac\{L\_\{\\mathcal\{X\}\}\}\{2\}\\alpha\_\{k\}^\{2\}\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}\\\|D\_\{k\}\\\|^\{2\}
Factoring out the squared inner product isolates the step parameter logic:

f​\(Xk\+1\)≤f​\(Xk\)−⟨Gk,Dk⟩2​\(αk−L𝒳2​αk2​‖Dk‖2\)f\(X\_\{k\+1\}\)\\leq f\(X\_\{k\}\)\-\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}\\left\(\\alpha\_\{k\}\-\\frac\{L\_\{\\mathcal\{X\}\}\}\{2\}\\alpha\_\{k\}^\{2\}\\\|D\_\{k\}\\\|^\{2\}\\right\)
By the absolute upper bound \(Condition 2\) and unitary invariance,‖Dk‖=‖pk​\(σk\)‖≤Mk\\\|D\_\{k\}\\\|=\\\|p\_\{k\}\(\\sigma\_\{k\}\)\\\|\\leq M\_\{k\}\. Substituting this guarantees the worst\-case penalty:

f​\(Xk\+1\)≤f​\(Xk\)−⟨Gk,Dk⟩2​\(αk−L𝒳2​αk2​Mk2\)f\(X\_\{k\+1\}\)\\leq f\(X\_\{k\}\)\-\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}\\left\(\\alpha\_\{k\}\-\\frac\{L\_\{\\mathcal\{X\}\}\}\{2\}\\alpha\_\{k\}^\{2\}M\_\{k\}^\{2\}\\right\)
Substituting the generalized step sizeαk=ηL𝒳​Mk2\\alpha\_\{k\}=\\frac\{\\eta\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}:

f​\(Xk\+1\)\\displaystyle f\(X\_\{k\+1\}\)≤f​\(Xk\)−⟨Gk,Dk⟩2​\(ηL𝒳​Mk2−L𝒳2​η2L𝒳2​Mk4​Mk2\)\\displaystyle\\leq f\(X\_\{k\}\)\-\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}\\left\(\\frac\{\\eta\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}\-\\frac\{L\_\{\\mathcal\{X\}\}\}\{2\}\\frac\{\\eta^\{2\}\}\{L\_\{\\mathcal\{X\}\}^\{2\}M\_\{k\}^\{4\}\}M\_\{k\}^\{2\}\\right\)=f​\(Xk\)−η​\(1−η2\)L𝒳​Mk2​⟨Gk,Dk⟩2\\displaystyle=f\(X\_\{k\}\)\-\\frac\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}\\langle G\_\{k\},D\_\{k\}\\rangle^\{2\}
By the scale\-invariant lower bound \(Condition 1\), the inner product is bounded below by the scaled dual norm of the gradient:⟨Gk,Dk⟩=⟨σk,pk​\(σk\)⟩≥mk​‖σk‖∗=mk​‖Gk‖∗\\langle G\_\{k\},D\_\{k\}\\rangle=\\langle\\sigma\_\{k\},p\_\{k\}\(\\sigma\_\{k\}\)\\rangle\\geq m\_\{k\}\\\|\\sigma\_\{k\}\\\|\_\{\*\}=m\_\{k\}\\\|G\_\{k\}\\\|\_\{\*\}\. Squaring this and substituting yields:

f​\(Xk\+1\)≤f​\(Xk\)−η​\(1−η2\)​mk2L𝒳​Mk2​‖Gk‖∗2f\(X\_\{k\+1\}\)\\leq f\(X\_\{k\}\)\-\\frac\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)m\_\{k\}^\{2\}\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}
Rearranging and summing overKKiterations gives a telescoping sum:

η​\(1−η2\)L𝒳​∑k=0K−1\(mkMk\)2​‖Gk‖∗2≤f​\(X0\)−f​\(XK\)≤f​\(X0\)−fmin\\frac\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)\}\{L\_\{\\mathcal\{X\}\}\}\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}\\leq f\(X\_\{0\}\)\-f\(X\_\{K\}\)\\leq f\(X\_\{0\}\)\-f\_\{\\min\}
We can lower\-bound the sum on the left by replacing each gradient norm with the minimum gradient norm encountered over theKKiterations:

\[min0≤k<K⁡‖Gk‖∗2\]​∑k=0K−1\(mkMk\)2≤∑k=0K−1\(mkMk\)2​‖Gk‖∗2\\left\[\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}\\right\]\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\\leq\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}
Substituting this into the inequality and isolating the minimum gradient norm perfectly recovers the finite\-time rate:

min0≤k<K⁡‖Gk‖∗2≤L𝒳​\(f​\(X0\)−fmin\)η​\(1−η2\)​∑k=0K−1\(mkMk\)2\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}\\leq\\frac\{L\_\{\\mathcal\{X\}\}\(f\(X\_\{0\}\)\-f\_\{\\min\}\)\}\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\}
Because the series∑\(mk/Mk\)2\\sum\(m\_\{k\}/M\_\{k\}\)^\{2\}diverges to infinity asK→∞K\\to\\infty, the right\-hand side converges to exactly zero\. Therefore,lim infk→∞‖Gk‖∗=0\\liminf\_\{k\\to\\infty\}\\\|G\_\{k\}\\\|\_\{\*\}=0\. ∎

###### Theorem D\.1\(Convergence of Stochastic Spectral Descent\)\.

Let the setting be identical to[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)\. Supposepk​\(σ\)p\_\{k\}\(\\sigma\)is generated by a stochastic process over\(k,σ\)∈ℕ×ℝ\>0r\(k,\\sigma\)\\in\\mathbb\{N\}\\times\\mathbb\{R\}^\{r\}\_\{\>0\}, and denote by𝔼k−1​\[⋅\]\\mathbb\{E\}\_\{k\-1\}\[\\cdot\]the conditional expectation given the history up to stepk−1k\-1\. Assume that the random variablesmk,Mkm\_\{k\},M\_\{k\}from[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)satisfy the following bounds almost surely,mk≥0m\_\{k\}\\geq 0andMk∈\(0,∞\)M\_\{k\}\\in\(0,\\infty\)\.

Letμk=𝔼k−1​\[\(mk/Mk\)2\]\\mu\_\{k\}=\\mathbb\{E\}\_\{k\-1\}\[\(m\_\{k\}/M\_\{k\}\)^\{2\}\], and assume∑k=0∞μk=∞\\sum\_\{k=0\}^\{\\infty\}\\mu\_\{k\}=\\inftyalmost surely\. Then, for step size chosen asαk=ηL𝒳​Mk2\\alpha\_\{k\}=\\frac\{\\eta\}\{L\_\{\\mathcal\{X\}\}M\_\{k\}^\{2\}\}forη∈\(0,2\)\\eta\\in\(0,2\), almost surely, the sequence achieves the finite\-time convergence rate:

min0≤k<K⁡‖Gk‖∗2=𝒪​\(1∑k=0K−1μk\)andlim infk→∞‖Gk‖∗=0\.\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}=\\mathcal\{O\}\\left\(\\frac\{1\}\{\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}\}\\right\)\\qquad\\text\{and\}\\qquad\\liminf\_\{k\\to\\infty\}\\\|G\_\{k\}\\\|\_\{\*\}=0\.

###### Proof\.

First, using Cauchy\-Schwarz, we have⟨σ,pk​\(σ\)⟩≤‖σ‖∗​‖pk​\(σ\)‖\\langle\\sigma,p\_\{k\}\(\\sigma\)\\rangle\\leq\\\|\\sigma\\\|\_\{\*\}\\\|p\_\{k\}\(\\sigma\)\\\|\. Substituting our assumed almost\-sure bounds into this inequality yields:

mk​‖σ‖∗≤⟨σ,pk​\(σ\)⟩≤‖σ‖∗​‖pk​\(σ\)‖≤Mk​‖σ‖∗m\_\{k\}\\\|\\sigma\\\|\_\{\*\}\\leq\\langle\\sigma,p\_\{k\}\(\\sigma\)\\rangle\\leq\\\|\\sigma\\\|\_\{\*\}\\\|p\_\{k\}\(\\sigma\)\\\|\\leq M\_\{k\}\\\|\\sigma\\\|\_\{\*\}Dividing by‖σ‖∗\\\|\\sigma\\\|\_\{\*\}reveals thatmk≤Mkm\_\{k\}\\leq M\_\{k\}almost surely\. Therefore, the sequence of random variablesZk=\(mk/Mk\)2Z\_\{k\}=\(m\_\{k\}/M\_\{k\}\)^\{2\}is non\-negative and uniformly bounded:Zk∈\[0,1\]Z\_\{k\}\\in\[0,1\]\.

We now invoke Lévy’s Extended Borel\-Cantelli Lemma for conditional expectations\. For a sequence of uniformly bounded, non\-negative random variablesZkZ\_\{k\}, the lemma states that the sum diverges almost surely and that their asymptotic ratio converges to exactly 1 almost surely if and only if the sum of their conditional expectations diverges:

∑k=0∞Zk=∞​and​limK→∞∑k=0K−1Zk∑k=0K−1𝔼k−1​\[Zk\]=1​a\.s\.⟺∑k=0∞𝔼k−1​\[Zk\]=∞\\sum\_\{k=0\}^\{\\infty\}Z\_\{k\}=\\infty\\;\\;\\text\{and\}\\;\\;\\lim\_\{K\\to\\infty\}\\frac\{\\sum\_\{k=0\}^\{K\-1\}Z\_\{k\}\}\{\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}\_\{k\-1\}\[Z\_\{k\}\]\}=1\\;\\;\\text\{a\.s\.\}\\quad\\Longleftrightarrow\\quad\\sum\_\{k=0\}^\{\\infty\}\\mathbb\{E\}\_\{k\-1\}\[Z\_\{k\}\]=\\inftyBy assumption, the sum of the conditional expectations diverges:∑k=0∞μk=∞\\sum\_\{k=0\}^\{\\infty\}\\mu\_\{k\}=\\infty\. By the ratio limit above, the actual observed sum of the random variables grows asymptotically at the exact same rate:

∑k=0K−1\(mkMk\)2=Ω​\(∑k=0K−1μk\)a\.s\.\\sum\_\{k=0\}^\{K\-1\}\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}=\\Omega\\left\(\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}\\right\)\\quad\\text\{a\.s\.\}
Since the random mappingspk​\(σ\)p\_\{k\}\(\\sigma\)almost surely satisfy the deterministic bounds of[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3)at every step, we can directly substitute this equivalent asymptotic growth into the denominator of the finite\-time rate derived in[Theorem˜2\.3](https://arxiv.org/html/2605.11181#S2.Thmtheorem3):

min0≤k<K⁡‖Gk‖∗2≤L𝒳​\(f​\(X0\)−fmin\)η​\(1−η2\)​Ω​\(∑k=0K−1μk\)=𝒪​\(1∑k=0K−1μk\)\\min\_\{0\\leq k<K\}\\\|G\_\{k\}\\\|\_\{\*\}^\{2\}\\leq\\frac\{L\_\{\\mathcal\{X\}\}\(f\(X\_\{0\}\)\-f\_\{\\min\}\)\}\{\\eta\(1\-\\frac\{\\eta\}\{2\}\)\\Omega\\left\(\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}\\right\)\}=\\mathcal\{O\}\\left\(\\frac\{1\}\{\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}\}\\right\)AsK→∞K\\to\\infty, the diverging sum in the denominator guarantees that this bound converges to zero, ensuring asymptotic convergence with probability11:lim infk→∞‖Gk‖∗=0\\liminf\_\{k\\to\\infty\}\\\|G\_\{k\}\\\|\_\{\*\}=0\. ∎

See[2\.5](https://arxiv.org/html/2605.11181#S2.Thmtheorem5)

###### Proof\.

Because all matrix norms are equivalent in finite dimensions, there exists a universal constantc\>0c\>0such that for any vectorv∈ℝrv\\in\\mathbb\{R\}^\{r\},‖v‖1≥c​‖v‖∗\\\|v\\\|\_\{1\}\\geq c\\\|v\\\|\_\{\*\}\. The upper bound is trivially constant from normalization:Mk=1M\_\{k\}=1\. Next, we verify the scale\-invariant lower bound\. The inner product is:

⟨σk,pk​\(σk\)⟩=∑i=1r\[σk\]i​\[Ek\]i‖Ek‖≥\[Ek\]min‖Ek‖​‖σk‖1\\langle\\sigma\_\{k\},p\_\{k\}\(\\sigma\_\{k\}\)\\rangle=\\sum\_\{i=1\}^\{r\}\[\\sigma\_\{k\}\]\_\{i\}\\frac\{\[E\_\{k\}\]\_\{i\}\}\{\\\|E\_\{k\}\\\|\}\\geq\\frac\{\[E\_\{k\}\]\_\{\\min\}\}\{\\\|E\_\{k\}\\\|\}\\\|\\sigma\_\{k\}\\\|\_\{1\}where\[Ek\]min=mini\[Ek\]i\[E\_\{k\}\]\_\{\\min\}=\\min\_\{i\}\[E\_\{k\}\]\_\{i\}\. Applying the norm equivalence‖σk‖1≥c​‖σk‖∗\\\|\\sigma\_\{k\}\\\|\_\{1\}\\geq c\\\|\\sigma\_\{k\}\\\|\_\{\*\}gives:

⟨σk,pk​\(σk\)⟩≥\(c​\[Ek\]min‖Ek‖\)​‖σk‖∗\\langle\\sigma\_\{k\},p\_\{k\}\(\\sigma\_\{k\}\)\\rangle\\geq\\left\(\\frac\{c\[E\_\{k\}\]\_\{\\min\}\}\{\\\|E\_\{k\}\\\|\}\\right\)\\\|\\sigma\_\{k\}\\\|\_\{\*\}We setmk=c​\[Ek\]min‖Ek‖m\_\{k\}=\\frac\{c\[E\_\{k\}\]\_\{\\min\}\}\{\\\|E\_\{k\}\\\|\}\. Because the vectorsEkE\_\{k\}are drawn entirely independently of the algorithmic history up to stepk−1k\-1, the conditional expectation equals the unconditional expectation\.

Because\[Ek\]min\[E\_\{k\}\]\_\{\\min\}is the minimum of a finite numberrrof positive random variables,mk\>0m\_\{k\}\>0\. Furthermore, sincemk/Mkm\_\{k\}/M\_\{k\}is bounded, its expectation exists and evaluates to a strictly positive constantμ\\mu:

μk=𝔼k−1​\[\(mkMk\)2\]=𝔼​\[\(c​\[Ek\]min‖Ek‖\)2\]=μ\>0\\mu\_\{k\}=\\mathbb\{E\}\_\{k\-1\}\\left\[\\left\(\\frac\{m\_\{k\}\}\{M\_\{k\}\}\\right\)^\{2\}\\right\]=\\mathbb\{E\}\\left\[\\left\(\\frac\{c\[E\_\{k\}\]\_\{\\min\}\}\{\\\|E\_\{k\}\\\|\}\\right\)^\{2\}\\right\]=\\mu\>0
Becauseμk=μ\\mu\_\{k\}=\\muis constant, the sum of expectations evaluates to∑k=0K−1μk=K​μ\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}=K\\mu, which diverges asK→∞K\\to\\infty\. The stochastic mapping thus satisfies all conditions of[Theorem˜D\.1](https://arxiv.org/html/2605.11181#A4.Thmtheorem1)\. Substituting∑μk=K​μ\\sum\\mu\_\{k\}=K\\muinto the generalized rate yields exactly:

𝒪​\(1∑k=0K−1μk\)=𝒪​\(1K​μ\)=𝒪​\(1K\)\\mathcal\{O\}\\left\(\\frac\{1\}\{\\sum\_\{k=0\}^\{K\-1\}\\mu\_\{k\}\}\\right\)=\\mathcal\{O\}\\left\(\\frac\{1\}\{K\\mu\}\\right\)=\\mathcal\{O\}\\left\(\\frac\{1\}\{K\}\\right\)Therefore, setting the step sizeαk=ηL𝒳\\alpha\_\{k\}=\\frac\{\\eta\}\{L\_\{\\mathcal\{X\}\}\}almost surely guarantees the𝒪​\(1/K\)\\mathcal\{O\}\(1/K\)finite\-time rate and asymptotic convergence\. ∎

### D\.2Equivariance ofFreon\(c=1\)\(c=1\)

###### Theorem D\.2\(Equivariance ofFreon\(c=1\)\(c=1\)\)\.

ForW∈ℝdW\\in\\mathbb\{R\}^\{d\}, letϕW\\phi\_\{W\}be an arbitrary function \(the neural network\) with a layered structure:W=\(W1,…,WL\)W=\(W\_\{1\},\\dots,W\_\{L\}\)withWl∈ℝnl×mlW\_\{l\}\\in\\mathbb\{R\}^\{n\_\{l\}\\times m\_\{l\}\}\. Letf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}be the differentiable loss function which depends onWWonly through the neural networkϕW\\phi\_\{W\}, and letGl∈ℝnl×mlG\_\{l\}\\in\\mathbb\{R\}^\{n\_\{l\}\\times m\_\{l\}\}be the gradient offffor layerll\. Then updates of the formWl−η​\(Gl​Gl⊤\)−1​GlW\_\{l\}\-\\eta\(G\_\{l\}G\_\{l\}^\{\\top\}\)^\{\-1\}G\_\{l\}are equivariant under symmetries ofϕW\\phi\_\{W\}which act on layer matrices by

Wl↦Ll​Wl​RlforLl∈ℝ∗​O​\(nl\),Rl∈ℝ∗​O​\(ml\)W\_\{l\}\\mapsto L\_\{l\}W\_\{l\}R\_\{l\}\\quad\\text\{for\}\\quad L\_\{l\}\\in\\mathbb\{R\}^\{\*\}O\(n\_\{l\}\),\\;R\_\{l\}\\in\\mathbb\{R\}^\{\*\}O\(m\_\{l\}\)whereℝ∗​O​\(ml\)=\{λ​Q:λ≠0,Q∈O​\(ml\)\}\\mathbb\{R\}^\{\*\}O\(m\_\{l\}\)=\\\{\\lambda Q:\\lambda\\neq 0,Q\\in O\(m\_\{l\}\)\\\}\. Therefore, the optimization trajectory ofϕW\\phi\_\{W\}is invariant under such symmetries\.

IfGlG\_\{l\}has full row \(column\) rank, the above holds more generally for symmetries withLl∈G​Lnl​\(ℝ\)L\_\{l\}\\in GL\_\{n\_\{l\}\}\(\\mathbb\{R\}\)\(Rl∈G​Lml​\(ℝ\)R\_\{l\}\\in GL\_\{m\_\{l\}\}\(\\mathbb\{R\}\)respectively\)\.

###### Proof\.

Write such a symmetrygg\. We have

f​\(g⋅W\)=f​\(W\)andg⋅\(W1,…,WL\)=\(L1​W1​R1,…,LL​WL​RL\)f\(g\\cdot W\)=f\(W\)\\quad\\text\{and\}\\quad g\\cdot\(W\_\{1\},\\dots,W\_\{L\}\)=\(\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{1\}W\_\{1\}\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{1\},\\dots,\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{L\}W\_\{L\}\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{L\}\)whereW=\(W1,…,WL\)W=\(W\_\{1\},\\dots,W\_\{L\}\)and someLl∈G​Lnl​\(ℝ\)\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\\in GL\_\{n\_\{l\}\}\(\\mathbb\{R\}\)andRl∈G​Lml​\(ℝ\)\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}\\in GL\_\{m\_\{l\}\}\(\\mathbb\{R\}\), and we use colours to denote the matrix symmetries\.

For1≤l≤L1\\leq l\\leq L, writeGl​\[W\]:=∇lf​\(W\)G\_\{l\}\[W\]:=\\nabla\_\{l\}f\(W\)where∇l\\nabla\_\{l\}denotes the gradient in thellthcoordinate\. By the chain rule we have

Gl​\[W\]=∇lf​\(W\)=∇Wlf​\(W\)=∇Wlf​\(g⋅W\)=Ll⊤​∇lf​\(g⋅W\)​Rl⊤=Ll⊤​Gl​\[g⋅W\]​Rl⊤\.G\_\{l\}\[W\]=\\nabla\_\{l\}f\(W\)=\\nabla\_\{W\_\{l\}\}f\(W\)=\\nabla\_\{W\_\{l\}\}f\(g\\cdot W\)=\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}^\{\\top\}\\nabla\_\{l\}f\(g\\cdot W\)\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}^\{\\top\}=\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}^\{\\top\}G\_\{l\}\[g\\cdot W\]\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}^\{\\top\}\.
IfLl∈O​\(nl\)\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\\in O\(n\_\{l\}\)andRl∈O​\(ml\)\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}\\in O\(m\_\{l\}\), note that

\(Gl​\[g⋅W\]​Gl​\[g⋅W\]⊤\)−1​Gl​\[g⋅W\]\\displaystyle\(G\_\{l\}\[g\\cdot W\]G\_\{l\}\[g\\cdot W\]^\{\\top\}\)^\{\-1\}G\_\{l\}\[g\\cdot W\]=\(Ll−⊤​Gl​\[W\]​Rl−⊤​Rl−1​Gl​\[W\]⊤​Ll−1\)−1​Ll−⊤​Gl​\[W\]​Rl−⊤\\displaystyle=\\left\(\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}^\{\-\\top\}G\_\{l\}\[W\]\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}^\{\-\\top\}R\_\{l\}^\{\-1\}G\_\{l\}\[W\]^\{\\top\}\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}^\{\-1\}\\right\)^\{\-1\}\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}^\{\-\\top\}G\_\{l\}\[W\]\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}^\{\-\\top\}\(7\)=Ll​\(Gl​\[W\]​Gl​\[W\]⊤\)−1​Gl​\[W\]​Rl\\displaystyle=\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\(G\_\{l\}\[W\]G\_\{l\}\[W\]^\{\\top\}\)^\{\-1\}G\_\{l\}\[W\]\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}where pseudoinverses are understood wherever necessary\. The same holds more generally forLl∈ℝ∗​O​\(nl\)\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\\in\\mathbb\{R\}^\{\*\}O\(n\_\{l\}\)andRl∈ℝ∗​O​\(ml\)\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}\\in\\mathbb\{R\}^\{\*\}O\(m\_\{l\}\), by cancelling out the scalar multiples which commute around all matrices\.

IfGl​\[W\]G\_\{l\}\[W\], or equivalentlyGl​\[g⋅W\]G\_\{l\}\[g\\cdot W\], has full row rank, then we have a genuine inverse in \([7](https://arxiv.org/html/2605.11181#A4.E7)\), so we can takeLl\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}out of the inverse for anyLl∈G​Lnl​\(ℝ\)\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\\in GL\_\{n\_\{l\}\}\(\\mathbb\{R\}\), and obtain the same result for suchLl\\color\[rgb\]\{1,\.5,0\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{1,\.5,0\}L\_\{l\}\.

IfGl​\[W\]G\_\{l\}\[W\]has full column rank, we can apply the same trick on the other side forRl∈G​Lml​\(ℝ\)\\color\[rgb\]\{0,1,1\}\\definecolor\[named\]\{pgfstrokecolor\}\{rgb\}\{0,1,1\}\\pgfsys@color@cmyk@stroke\{1\}\{0\}\{0\}\{0\}\\pgfsys@color@cmyk@fill\{1\}\{0\}\{0\}\{0\}R\_\{l\}\\in GL\_\{m\_\{l\}\}\(\\mathbb\{R\}\)by first using the identity

\(Gl​\[g⋅W\]​Gl​\[g⋅W\]⊤\)−1​Gl​\[g⋅W\]=Gl​\[g⋅W\]​\(Gl​\[g⋅W\]⊤​Gl​\[g⋅W\]\)−1\.\(G\_\{l\}\[g\\cdot W\]G\_\{l\}\[g\\cdot W\]^\{\\top\}\)^\{\-1\}G\_\{l\}\[g\\cdot W\]=G\_\{l\}\[g\\cdot W\]\(G\_\{l\}\[g\\cdot W\]^\{\\top\}G\_\{l\}\[g\\cdot W\]\)^\{\-1\}\.
Thus, under such symmetries, the updates are equivariant in weight space, and hence they are invariant in output space\. ∎

###### Note D\.3\.

These encapsulate ReLU scaling symmetries, node permutation symmetries, attention linear symmetries, etc, and most particularly transformations that are symmetries only on sampled data, even if the parameterization itself does not admit it globally\.

##### Comparison with Newton’s method\.

The preconditionerG​G⊤GG^\{\\top\}to the update step inFreon\(1/1\)\(1/1\)can seem reminiscent to the Fisher information matrix used as preconditioner in Newton’s method\. It is well\-known that Newton’s steps are equivariant under*all*linear symmetries \(not just layerwise ones\)\. Here, we briefly clarify the connection and distinction between the two optimizers\.

Consider a general deep networkϕW:ℝnin→ℝnout\\phi\_\{W\}:\\mathbb\{R\}^\{n\_\{\\text\{in\}\}\}\\to\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\}, with parametersWW\. As before,f:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}is the loss function\.

In the following we abstract away the multivalued data and considerϕW\\phi\_\{W\}implicitly evaluated solely at one inputx∈ℝninx\\in\\mathbb\{R\}^\{n\_\{\\text\{in\}\}\}, the general case follows by summing or averaging over such inputs\. The lossffis viewed as a function ofϕW\\phi\_\{W\}\.

Define the quantitiesgϕ:=∇WϕW∈ℝnoutg\_\{\\phi\}:=\\nabla\_\{W\}\\phi\_\{W\}\\in\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\},hf:=∇ϕW∇ϕW⊤⁡f∈ℝnout×nouth\_\{f\}:=\\nabla\_\{\\phi\_\{W\}\}\\nabla\_\{\\phi\_\{W\}\}^\{\\top\}f\\in\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\\times n\_\{\\text\{out\}\}\},gf:=∇Wf∈ℝdg\_\{f\}:=\\nabla\_\{W\}f\\in\\mathbb\{R\}^\{d\},Jϕ:=∇WϕW∈ℝnout×dJ\_\{\\phi\}:=\\nabla\_\{W\}\\phi\_\{W\}\\in\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\\times d\},Gl:=∇Wlf∈ℝnl×mlG\_\{l\}:=\\nabla\_\{W\_\{l\}\}f\\in\\mathbb\{R\}^\{n\_\{l\}\\times m\_\{l\}\},Jl:=∇WlϕW∈ℝnout×\(nl×ml\)J\_\{l\}:=\\nabla\_\{W\_\{l\}\}\\phi\_\{W\}\\in\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\\times\(n\_\{l\}\\times m\_\{l\}\)\}\.

By the chain rulegϕ=Jϕ⊤​gfg\_\{\\phi\}=J\_\{\\phi\}^\{\\top\}g\_\{f\}\. For the usuall2l\_\{2\}or softmax\-cross entropy losses, we have the probabilistic interpretationf​\(W\)=−log⁡p​\(y∗∣x,W\)f\(W\)=\-\\log p\(y^\{\*\}\\mid x,W\)wherey∗∈ℝnouty^\{\*\}\\in\\mathbb\{R\}^\{n\_\{\\text\{out\}\}\}is the label associated toxx\. The Fisher information matrix \(or GGN\) is given by

𝔼p​\(y∣x,W\)​\[−∇W∇W⊤log⁡p​\(y∣x,W\)\]=Jϕ⊤​hϕ​Jϕ\.\\mathbb\{E\}\_\{p\(y\\mid x,W\)\}\[\-\\nabla\_\{W\}\\nabla\_\{W\}^\{\\top\}\\log p\(y\\mid x,W\)\]=J\_\{\\phi\}^\{\\top\}h\_\{\\phi\}J\_\{\\phi\}\.A common approximation is the empirical Fisher\[Kunstneret al\.,[2019](https://arxiv.org/html/2605.11181#bib.bib67)\]:

𝔼p​\(y∣x,W\)​\[−∇W∇W⊤log⁡p​\(y∣x,W\)\]\\displaystyle\\mathbb\{E\}\_\{p\(y\\mid x,W\)\}\[\-\\nabla\_\{W\}\\nabla\_\{W\}^\{\\top\}\\log p\(y\\mid x,W\)\]=𝔼p​\(y∣x,W\)​\[∇Wlog⁡p​\(y∣x,W\)​∇Wlog⁡p​\(y∣x,W\)⊤\]\\displaystyle=\\mathbb\{E\}\_\{p\(y\\mid x,W\)\}\[\\nabla\_\{W\}\\log p\(y\\mid x,W\)\\nabla\_\{W\}\\log p\(y\\mid x,W\)^\{\\top\}\]\(8\)≈∇Wlog⁡p​\(y∗∣x,W\)​∇Wlog⁡p​\(y∗∣x,W\)⊤\\displaystyle\\approx\\nabla\_\{W\}\\log p\(y^\{\*\}\\mid x,W\)\\nabla\_\{W\}\\log p\(y^\{\*\}\\mid x,W\)^\{\\top\}=gf​gf⊤\\displaystyle=g\_\{f\}g\_\{f\}^\{\\top\}=Jϕ⊤​gf​gf⊤​Jϕ\\displaystyle=J\_\{\\phi\}^\{\\top\}g\_\{f\}g\_\{f\}^\{\\top\}J\_\{\\phi\}This is the preconditioner used in Newton’s method\.

Now by the chain ruleGl=Mat\(Jl⊤​gf\)G\_\{l\}=\\operatorname\*\{Mat\}\(J\_\{l\}^\{\\top\}g\_\{f\}\), where theMat\\operatorname\*\{Mat\}operator makes the output a matrix of the appropriate shape\. The preconditioner inFreon\(1/1\)\(1/1\)is

Gl​Gl⊤=Mat\(Jl⊤​gf\)​Mat\(gf⊤​Jl\)\.G\_\{l\}G\_\{l\}^\{\\top\}=\\operatorname\*\{Mat\}\(J\_\{l\}^\{\\top\}g\_\{f\}\)\\operatorname\*\{Mat\}\(g\_\{f\}^\{\\top\}J\_\{l\}\)\.\(9\)
Compare \([9](https://arxiv.org/html/2605.11181#A4.E9)\) with \([8](https://arxiv.org/html/2605.11181#A4.E8)\)\.

## Appendix EFrom Polar to Hölder Express

### E\.1Fundamental Properties of thea/ba/b\-Method

We aim in this section to generalise the results inAmselet al\.\[[2026](https://arxiv.org/html/2605.11181#bib.bib12)\]to our setup\. To this end, letr,p∈ℕr,p\\in\\mathbb\{N\}and consider the following subsets of the polynomial ringℙ\\mathbb\{P\}

ℙodd​\(r\):=span​\{x2​i​r\+1:i∈ℕ\},\\displaystyle\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\):=\\mathrm\{span\}\\\{x^\{2ir\+1\}:i\\in\\mathbb\{N\}\\\},\(10\)ℙeven​\(r\):=span​\{x2​i​r:i∈ℕ\},\\displaystyle\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\):=\\mathrm\{span\}\\\{x^\{2ir\}:i\\in\\mathbb\{N\}\\\},\(11\)ℙpodd​\(r\):=\{P∈ℙodd​\(r\):deg⁡\(P\)≤2​r​p\+1\},\\displaystyle\\mathbb\{P\}^\{\\mathrm\{odd\}\}\_\{p\}\(r\):=\\\{P\\in\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\):\\deg\(P\)\\leq 2rp\+1\\\},\(12\)ℙpeven​\(r\):=\{P∈ℙeven​\(r\):deg⁡\(P\)≤2​r​p\}\.\\displaystyle\\mathbb\{P\}^\{\\mathrm\{even\}\}\_\{p\}\(r\):=\\\{P\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\):\\deg\(P\)\\leq 2rp\\\}\.\(13\)Moreover, for anyn,d∈ℕn,d\\in\\mathbb\{N\}, we also define the following sets of rational functions

ℜ​\(r,s\):=\{ND:N∈ℙodd​\(r\),D∈ℙeven​\(s\),D\>0\},\\displaystyle\\mathfrak\{R\}\(r,s\):=\\left\\\{\\frac\{N\}\{D\}:N\\in\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\),D\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(s\),\\,D\>0\\right\\\},\(14\)ℜn,d​\(r,s\):=\{ND:N∈ℙno​d​d​\(r\),D∈ℙde​v​e​n​\(s\),D\>0\}\.\\displaystyle\\mathfrak\{R\}\_\{n,d\}\(r,s\):=\\left\\\{\\frac\{N\}\{D\}:N\\in\\mathbb\{P\}\_\{n\}^\{odd\}\(r\),D\\in\\mathbb\{P\}\_\{d\}^\{even\}\(s\),\\,D\>0\\right\\\}\.\(15\)
#### E\.1\.1Existence and Characterization of the Optimal Approximants\.

The setℜ​\(r,s\)\\mathfrak\{R\}\(r,s\)considered above can be understood as a particular subset of generalized rational functions of the form

R​\(x\)=∑i=0nai​gi​\(x\)∑j=0dbj​hj​\(x\),∑j=0dbj​hj​\(x\)\>0,R\(x\)=\\frac\{\\sum\_\{i=0\}^\{n\}a\_\{i\}g\_\{i\}\(x\)\}\{\\sum\_\{j=0\}^\{d\}b\_\{j\}h\_\{j\}\(x\)\},\\qquad\\sum\_\{j=0\}^\{d\}b\_\{j\}h\_\{j\}\(x\)\>0,wheregi​\(x\)=x2​r​i\+1g\_\{i\}\(x\)=x^\{2ri\+1\}andhj​\(x\)=x2​s​jh\_\{j\}\(x\)=x^\{2sj\}\. For these objects, unlike the classical polynomial case, the existence of an optimal approximant is not guaranteed in general\. Furthermore, additional issues arise from the fact that the numerator and denominator of the approximantsRRcan have a common factor, which leads to simplifications that might break the constraintR∈ℜn,d​\(r\)R\\in\\mathfrak\{R\}\_\{n,d\}\(r\)\. This phenomenon is studied in the literature through the so called deficiency index𝔡\\mathfrak\{d\}which strictly depends on the class of functions we want to approximate\. In this regards, let us define the following class of functions

###### Definition E\.1\.

Leta,b∈ℝa,b\\in\\mathbb\{R\}\. We define the set of0\-deficiency continuous functionsC𝔡=00​\(\[a,b\]\)⊂C0​\(\[a,b\]\)C\_\{\\mathfrak\{d\}=0\}^\{0\}\(\[a,b\]\)\\subset C^\{0\}\(\[a,b\]\)as

C𝔡=00​\(\[a,b\]\):=\{f∈C0​\(\[a,b\]\):𝔡f=0\}\.C\_\{\\mathfrak\{d\}=0\}^\{0\}\(\[a,b\]\):=\\\{f\\in C^\{0\}\(\[a,b\]\):\\mathfrak\{d\}\_\{f\}=0\\\}\.

We can prove the following result\.

###### Lemma E\.2\.

Letr∈ℕr\\in\\mathbb\{N\}and considerℜ​\(r,r\)\\mathfrak\{R\}\(r,r\)as in \([15](https://arxiv.org/html/2605.11181#A5.E15)\)\. Moreover, let0<a≤b∈ℝ0<a\\leq b\\in\\mathbb\{R\}and define

𝒞0​\(\[a,b\]\):=\{f∈C0​\(\[a,b\]\):f​\(x12​r\)x12​r∈C𝔡=00​\(\[a,b\]\)\}\.\\mathscr\{C\}^\{0\}\(\[a,b\]\):=\\left\\\{f\\in C^\{0\}\(\[a,b\]\):\\frac\{f\(x^\{\\frac\{1\}\{2r\}\}\)\}\{x^\{\\frac\{1\}\{2r\}\}\}\\in C\_\{\\mathfrak\{d\}=0\}^\{0\}\(\[a,b\]\)\\right\\\}\.Then, for everyf∈𝒞0​\(\[a,b\]\)f\\in\\mathscr\{C\}^\{0\}\(\[a,b\]\), there existsR∗∈ℜn,d​\(r,r\)R^\{\*\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)such that

‖f−R∗‖C0​\(\[a,b\]\)≤‖f−R‖C0​\(\[a,b\]\),∀R∈ℜn,d​\(r\)\.\\\|f\-R^\{\*\}\\\|\_\{C^\{0\}\(\[a,b\]\)\}\\leq\\\|f\-R\\\|\_\{C^\{0\}\(\[a,b\]\)\},\\qquad\\forall R\\in\\mathfrak\{R\}\_\{n,d\}\(r\)\.

###### Proof\.

Let us observe that anyR∈ℜn,d​\(r,r\)R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)can be rewritten as

R\(x\)=∑i=0nai​x2​r​i\+1∑j=0dbj​x2​r​j=x∑i=0nai​x2​r​i∑j=0dbj​x2​r​j=:xP​\(x2​r\)Q​\(x2​r\)R\(x\)=\\frac\{\\sum\_\{i=0\}^\{n\}a\_\{i\}x^\{2ri\+1\}\}\{\\sum\_\{j=0\}^\{d\}b\_\{j\}x^\{2rj\}\}=x\\frac\{\\sum\_\{i=0\}^\{n\}a\_\{i\}x^\{2ri\}\}\{\\sum\_\{j=0\}^\{d\}b\_\{j\}x^\{2rj\}\}=:x\\frac\{P\(x^\{2r\}\)\}\{Q\(x^\{2r\}\)\}\(16\)whereP∈ℙnP\\in\\mathbb\{P\}\_\{n\}andQ∈ℙdQ\\in\\mathbb\{P\}\_\{d\}\. ByAchieser \[[1992](https://arxiv.org/html/2605.11181#bib.bib51)\], we know that for everyg∈C0​\(\[a,b\]\)g\\in C^\{0\}\(\[a,b\]\), there always existP∗∈ℙnP^\{\*\}\\in\\mathbb\{P\}\_\{n\},Q∗∈ℙdQ^\{\*\}\\in\\mathbb\{P\}\_\{d\}such thatP∗Q∗\\frac\{P^\{\*\}\}\{Q^\{\*\}\}best approximatesgg\. Then, since0<a≤b0<a\\leq b, consideringf∈C0​\(\[a,b\]\)f\\in C^\{0\}\(\[a,b\]\)and setting

g~​\(x\)=f​\(x12​r\)x12​r,\\widetilde\{g\}\(x\)=\\frac\{f\\left\(x^\{\\frac\{1\}\{2r\}\}\\right\)\}\{x^\{\\frac\{1\}\{2r\}\}\},\(17\)we get a best approximation candidateR∗​\(x\)=x​P∗​\(x\)Q∗​\(x\)R^\{\*\}\(x\)=x\\frac\{P^\{\*\}\(x\)\}\{Q^\{\*\}\(x\)\}forff\.

To conclude, by Definition[E\.1](https://arxiv.org/html/2605.11181#A5.Thmtheorem1), whenf∈𝒞0​\(\[a,b\]\)f\\in\\mathscr\{C\}^\{0\}\(\[a,b\]\), we have𝔡g~=0\\mathfrak\{d\}\_\{\\widetilde\{g\}\}=0and we can ensureR∗∈ℜn,d​\(r,r\)R^\{\*\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\. ∎

Finally, mimicking the polynomial case, we are interested in providing a characterization of the optimal approximants\. The following holds\.

###### Lemma E\.3\.

Let us assume the same setup as in Lemma[E\.2](https://arxiv.org/html/2605.11181#A5.Thmtheorem2)\. Then, for everyf∈𝒞0​\(\[a,b\]\)f\\in\\mathscr\{C\}^\{0\}\(\[a,b\]\), anyR∗∈ℜn,d​\(r,r\)R^\{\*\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)such that

\#​\{x∈\[a,b\]:\|f​\(x\)−R∗​\(x\)\|=‖f−R∗‖C0\}≥n\+d\+2\\\#\\\{x\\in\[a,b\]:\|f\(x\)\-R^\{\*\}\(x\)\|=\\\|f\-R^\{\*\}\\\|\_\{C^\{0\}\}\\\}\\geq n\+d\+2\(18\)is an optimal approximant offf\.

###### Proof\.

Givenf∈𝒞0​\(\[a,b\]\)f\\in\\mathscr\{C\}^\{0\}\(\[a,b\]\), exploiting \([16](https://arxiv.org/html/2605.11181#A5.E16)\), any optimalR∗∈ℜn,d​\(r,r\)R^\{\*\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)fulfils the same properties thatP∗/Q∗P^\{\*\}/Q^\{\*\}satisfies approximatingg~\\widetilde\{g\}\. Then, by\[Loeb,[1966](https://arxiv.org/html/2605.11181#bib.bib50), Theorem 1\], we have the thesis\. ∎

#### E\.1\.2Optimality

We will now address the problem of the optimality of the \. To start with, the following closure property holds\.

###### Lemma E\.4\.

Letr∈ℕr\\in\\mathbb\{N\}\. Then,ℙodd​\(r\)\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\),ℙeven​\(r\)\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)andℜ​\(r,r\)\\mathfrak\{R\}\(r,r\)are closed with respect to composition\.

###### Proof\.

Let us considerℙodd​\(r\)\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\)\. By induction on the monomials, it suffices to show that

p=p1∘p2∈ℙodd​\(r\)p=p\_\{1\}\\circ p\_\{2\}\\in\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\)forp1=x2​l​r\+1p\_\{1\}=x^\{2lr\+1\},p2=x2​n​r\+1\+x2​m​r\+1p\_\{2\}=x^\{2nr\+1\}\+x^\{2mr\+1\}\. Recalling now Newton’s binomial formula, it holds

p​\(x\)\\displaystyle p\(x\)=∑i=02​l​r\+1\(x2​n​r\+1\)i​\(x2​m​b\+1\)2​l​b\+1−i\\displaystyle=\\sum\_\{i=0\}^\{2lr\+1\}\\left\(x^\{2nr\+1\}\\right\)^\{i\}\\left\(x^\{2mb\+1\}\\right\)^\{2lb\+1\-i\}=∑i=02​l​r\+1x2​n​r​i\+i\+2​m​\(2​l​r\+1−i\)​r\+k−i\\displaystyle=\\sum\_\{i=0\}^\{2lr\+1\}x^\{2nri\+i\+2m\(2lr\+1\-i\)r\+k\-i\}=∑i=02​l​r\+1x2​\[n​r​i\+m​\(2​l​r\+1−i\)\+2​l\]​r\+1\.\\displaystyle=\\sum\_\{i=0\}^\{2lr\+1\}x^\{2\\left\[nri\+m\(2lr\+1\-i\)\+2l\\right\]r\+1\}\.Then, we conclude by observing that, sincen​r​i\+m​\(2​l​r\+1−i\)\+2​l∈ℕnri\+m\(2lr\+1\-i\)\+2l\\in\\mathbb\{N\}, every term in the above expansion is inℙodd​\(r\)\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\)\. Exploiting an analogous argument as above, we can also prove the following

1. 1\.ifp1p\_\{1\},p2∈ℙeven​\(r\)p\_\{2\}\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\), thenp∈ℙeven​\(r\)p\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)\(which implies the closure ofℙeven​\(r\)\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)as claimed\);
2. 2\.ifp1∈ℙodd​\(r\)p\_\{1\}\\in\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\),p2∈ℙeven​\(r\)p\_\{2\}\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)\(or vice versa\), thenp∈ℙeven​\(r\)p\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)\.

Let us considerq1q\_\{1\},q2∈ℜn,d​\(r,r\)q\_\{2\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\. We have

D12​n​r\+1​\(x\)​\(N2∘q1\)​\(x\)\\displaystyle D\_\{1\}^\{2nr\+1\}\(x\)\(N\_\{2\}\\circ q\_\{1\}\)\(x\)=∑i=1nci​N12​i​r\+1​\(x\)​D12​r​\(n−i\)​\(x\),\\displaystyle=\\sum\_\{i=1\}^\{n\}c\_\{i\}N\_\{1\}^\{2ir\+1\}\(x\)D\_\{1\}^\{2r\(n\-i\)\}\(x\),D12​d​r​\(x\)​\(D2∘q1\)​\(x\)\\displaystyle D\_\{1\}^\{2dr\}\(x\)\(D\_\{2\}\\circ q\_\{1\}\)\(x\)=∑i=1dei​N12​i​r​\(x\)​D12​r​\(d−i\)​\(x\)\.\\displaystyle=\\sum\_\{i=1\}^\{d\}e\_\{i\}N\_\{1\}^\{2ir\}\(x\)D\_\{1\}^\{2r\(d\-i\)\}\(x\)\.The arguments above implies then thatN12​i​r\+1∈ℙodd​\(r\)N\_\{1\}^\{2ir\+1\}\\in\\mathbb\{P\}^\{\\mathrm\{odd\}\}\(r\)whileD12​r​\(n−i\)D\_\{1\}^\{2r\(n\-i\)\},N12​i​rN\_\{1\}^\{2ir\},D12​r​\(d−i\)∈ℙeven​\(r\)D\_\{1\}^\{2r\(d\-i\)\}\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\)\. Finally, recalling that the product of even polynomials is even and the product of an even polynomial and an odd one is odd, since onlyrr\-powers monomial appears, we immediately getN2∘q1∈ℜ​\(r,r\)N\_\{2\}\\circ q\_\{1\}\\in\\mathfrak\{R\}\(r,r\)andD12​d​r​\(D2∘q1\)∈ℙeven​\(r\)D\_\{1\}^\{2dr\}\(D\_\{2\}\\circ q\_\{1\}\)\\in\\mathbb\{P\}^\{\\mathrm\{even\}\}\(r\), so that

\(q1∘q2\)​\(x\)=\(N2∘q1\)​\(x\)D12​\(d−n\)​r−1​\(x\)​\(D2∘q1\)​\(x\)∈ℜn,d​\(r,r\),\(q\_\{1\}\\circ q\_\{2\}\)\(x\)=\\frac\{\(N\_\{2\}\\circ q\_\{1\}\)\(x\)\}\{D\_\{1\}^\{2\(d\-n\)r\-1\}\(x\)\(D\_\{2\}\\circ q\_\{1\}\)\(x\)\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\),for allnn,d∈ℕd\\in\\mathbb\{N\}\. ∎

We can then generalise\[Amselet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib12), Theorem 3\.1\]with the following optimality result\.

###### Theorem E\.5\.

Letl1=l,u1=u∈\(0,1\)l\_\{1\}=l,u\_\{1\}=u\\in\(0,1\),r,T∈ℕr,T\\in\\mathbb\{N\},n,d∈ℕn,d\\in\\mathbb\{N\}and consider

Rt\+1=argminR∈ℜn,d​\(r,r\)​maxx∈\[lt,ut\]⁡\|1−R​\(x\)\|,∀t∈\[1,T\],R\_\{t\+1\}=\\underset\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\{\\mathrm\{argmin\}\}\\max\_\{x\\in\[l\_\{t\},u\_\{t\}\]\}\|1\-R\(x\)\|,\\qquad\\forall t\\in\[1,T\],withlt\+1=minx∈\[lt,ut\]⁡Rt​\(x\)l\_\{t\+1\}=\\min\_\{x\\in\[l\_\{t\},u\_\{t\}\]\}R\_\{t\}\(x\),ut\+1=maxx∈\[lt,ut\]⁡Rt​\(x\)u\_\{t\+1\}=\\max\_\{x\\in\[l\_\{t\},u\_\{t\}\]\}R\_\{t\}\(x\)\. Then,R∗:=RT∘RT−1∘⋯∘R1R^\{\*\}:=R\_\{T\}\\circ R\_\{T\-1\}\\circ\\cdots\\circ R\_\{1\}is optimal and

maxx∈\[lT,uT\]⁡\|1−R∗​\(x\)\|=1−lT\.\\max\_\{x\\in\[l\_\{T\},u\_\{T\}\]\}\|1\-R^\{\*\}\(x\)\|=1\-l\_\{T\}\.Moreover,

\{lt\+1=Rt​\(lt\),ut\+1=2−Rt​\(lt\),maxx∈\[lt,ut\]⁡\|1−Rt​\(x\)\|=1−lt,∀t∈\[1,T\]\.\\left\\\{\\begin\{aligned\} &l\_\{t\+1\}=R\_\{t\}\(l\_\{t\}\),\\\\ &u\_\{t\+1\}=2\-R\_\{t\}\(l\_\{t\}\),\\\\ &\\max\_\{x\\in\[l\_\{t\},u\_\{t\}\]\}\|1\-R\_\{t\}\(x\)\|=1\-l\_\{t\},\\end\{aligned\}\\right\.\\qquad\\forall t\\in\[1,T\]\.\(19\)

###### Proof\.

Let us first observe that, byWall \[[1948](https://arxiv.org/html/2605.11181#bib.bib34)\], we have that the Hankel determinants ofx↦xαx\\mapsto x^\{\\alpha\},α∉ℕ\\alpha\\not\\in\\mathbb\{N\}, are all non\-zero\. Thus, by\[Gragg,[1972](https://arxiv.org/html/2605.11181#bib.bib10), Corollary 2\],f​\(x\)≡1∈𝒞0​\(\[a,b\]\)f\(x\)\\equiv 1\\in\\mathscr\{C\}^\{0\}\(\[a,b\]\), so that, thanks to Lemma[E\.2](https://arxiv.org/html/2605.11181#A5.Thmtheorem2),RtR\_\{t\}\(and consequentlyR∗R^\{\*\}\) is well defined\.

We will now adapt the argument used in\[Amselet al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib12), Theorem 3\.1\], and split our discussion into two parts: first, we prove \([19](https://arxiv.org/html/2605.11181#A5.E19)\) by showing thatltl\_\{t\}is a minimizer ofRtR\_\{t\}; second, we argue on the optimality ofR∗R^\{\*\}\.

##### ltl\_\{t\}is a minimizer ofRtR\_\{t\}\.

By direct computation, we have

Rt′​\(x\)\\displaystyle R\_\{t\}^\{\\prime\}\(x\)=N′​\(x\)​D​\(x\)−D′​\(x\)​P​\(x\)D2​\(x\)\\displaystyle=\\frac\{N^\{\\prime\}\(x\)D\(x\)\-D^\{\\prime\}\(x\)P\(x\)\}\{D^\{2\}\(x\)\}=\(∑i=0nai​\(2​r​i\+1\)​x2​r​i\)​\(∑j=0dbj​x2​r​j\)−\(∑j=0d2​r​j​bj​x2​r​j−1\)​\(∑i=0nai​x2​r​i\+1\)D2​\(x\)\\displaystyle=\\frac\{\\left\(\\sum\_\{i=0\}^\{n\}a\_\{i\}\(2ri\+1\)x^\{2ri\}\\right\)\\left\(\\sum\_\{j=0\}^\{d\}b\_\{j\}x^\{2rj\}\\right\)\-\\left\(\\sum\_\{j=0\}^\{d\}2rjb\_\{j\}x^\{2rj\-1\}\\right\)\\left\(\\sum\_\{i=0\}^\{n\}a\_\{i\}x^\{2ri\+1\}\\right\)\}\{D^\{2\}\(x\)\}=∑i=0n∑j=0d\(2​r​\(i−j\)\+1\)​ai​bj​x2​r​\(i\+j\)D2​\(x\)=N~​\(x2​r\)D~​\(x2​r\)\\displaystyle=\\frac\{\\sum\_\{i=0\}^\{n\}\\sum\_\{j=0\}^\{d\}\(2r\(i\-j\)\+1\)a\_\{i\}b\_\{j\}x^\{2r\(i\+j\)\}\}\{D^\{2\}\(x\)\}=\\frac\{\\widetilde\{N\}\(x^\{2r\}\)\}\{\\widetilde\{D\}\(x^\{2r\}\)\}whereN~∈ℙn\\widetilde\{N\}\\in\\mathbb\{P\}\_\{n\},D~∈ℙd\\widetilde\{D\}\\in\\mathbb\{P\}\_\{d\}\. Thus,RtR\_\{t\}has at most2​\(n\+d\)2\(n\+d\)extremal points, of which, due to the symmetry ofN~​\(x2​r\)\\widetilde\{N\}\(x^\{2r\}\), at mostn\+dn\+dcan be contained in\[lt,ut\]⊂\[0,1\]\[l\_\{t\},u\_\{t\}\]\\subset\[0,1\]\. Exploiting now Lemma[E\.3](https://arxiv.org/html/2605.11181#A5.Thmtheorem3), we know that, beingRtR\_\{t\}an optimal approximant, there exist at leastn\+d\+2n\+d\+2extremal points for\|1−Rt\|\|1\-R\_\{t\}\|on\[lt,ut\]\[l\_\{t\},u\_\{t\}\]\. In particular, we can conclude that the limiting pointsltl\_\{t\}andutu\_\{t\}need to be extremal points\.

SinceRt​\(0\)=0R\_\{t\}\(0\)=0, beingRtR\_\{t\}optimal, we haveRt\>0R\_\{t\}\>0on\[lt,ut\]\[l\_\{t\},u\_\{t\}\]otherwise,Rt=0R\_\{t\}=0would be a better approximant\. In particular, there existsx^∈\[0,lt\]\\hat\{x\}\\in\[0,l\_\{t\}\]such thatRt′​\(x^\)\>0R^\{\\prime\}\_\{t\}\(\\hat\{x\}\)\>0\. If we suppose now thatltl\_\{t\}is a maximum ofRtR\_\{t\}, thenRt′​\(lt\)≤0R^\{\\prime\}\_\{t\}\(l\_\{t\}\)\\leq 0and, by the intermediate value theorem, there must existx¯∈\[x^,lt\]\\overline\{x\}\\in\[\\hat\{x\},l\_\{t\}\]such thatRt​\(x¯\)=0R\_\{t\}\(\\overline\{x\}\)=0\. This is, however, a contradiction since all the nonnegative extremal points ofRtR\_\{t\}are in\[lt,ut\]\[l\_\{t\},u\_\{t\}\]so thatltl\_\{t\}must be a minimum ofRtR\_\{t\}\. Finally, by direct computation, \([19](https://arxiv.org/html/2605.11181#A5.E19)\) immediately follows\.

##### Optimality ofR∗R^\{\*\}\.

We will argue by induction\. Let us observe that Lemma[E\.2](https://arxiv.org/html/2605.11181#A5.Thmtheorem2)and the previous paragraph ensure the base step\. Then, let us assume that the thesis holds forT−1T\-1and considerf:\[l,u\]→ℝf:\[l,u\]\\to\\mathbb\{R\}given byx↦R¯T−1∘⋯∘R¯1​\(x\)x\\mapsto\\overline\{R\}\_\{T\-1\}\\circ\\cdots\\circ\\overline\{R\}\_\{1\}\(x\), for someR¯1,…​R¯T−1∈ℜn,d​\(r,r\)\\overline\{R\}\_\{1\},\.\.\.\\overline\{R\}\_\{T\-1\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\. Moreover, letting\[a,b\]\[a,b\]be the image of\[l,u\]\[l,u\]throughff, for any constantccwe immediately have

maxx∈\[l,u\]⁡\|1−c​f​\(x\)\|=max⁡\{1−c​a,c​b−1\}\.\\max\_\{x\\in\[l,u\]\}\|1\-cf\(x\)\|=\\max\\\{1\-ca,cb\-1\\\}\.\(20\)Let us suppose now thatab\>ltut\\frac\{a\}\{b\}\>\\frac\{l\_\{t\}\}\{u\_\{t\}\}then the following holds

ab\>lt2−lt⇔1−lt\>b−aa\+b\.\\displaystyle\\frac\{a\}\{b\}\>\\frac\{l\_\{t\}\}\{2\-l\_\{t\}\}\\quad\\Leftrightarrow\\quad 1\-l\_\{t\}\>\\frac\{b\-a\}\{a\+b\}\.Then, since by the previous discussionlt=minx∈\[l,u\]⁡RT−1∘⋯∘R1​\(x\)l\_\{t\}=\\min\_\{x\\in\[l,u\]\}R\_\{T\-1\}\\circ\\cdots\\circ R\_\{1\}\(x\), by choosingc¯=2a\+b\\overline\{c\}=\\frac\{2\}\{a\+b\}in \([20](https://arxiv.org/html/2605.11181#A5.E20)\), the above implies

maxx∈\[l,u\]⁡\|1−RT−1∘⋯∘R1​\(x\)\|\>maxx∈\[l,u\]⁡\|1−c¯​f​\(x\)\|\\max\_\{x\\in\[l,u\]\}\|1\-R\_\{T\-1\}\\circ\\cdots\\circ R\_\{1\}\(x\)\|\>\\max\_\{x\\in\[l,u\]\}\|1\-\\overline\{c\}f\(x\)\|contradicting the inductive hypothesis\. Finally, sinceab<ltut\\frac\{a\}\{b\}<\\frac\{l\_\{t\}\}\{u\_\{t\}\}, for anyR¯T∈ℜn,d​\(r,r\)\\overline\{R\}\_\{T\}\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\), we have

maxx∈\[l,u\]⁡\|1−R¯T​\(f​\(x\)\)\|\\displaystyle\\max\_\{x\\in\[l,u\]\}\|1\-\\overline\{R\}\_\{T\}\(f\(x\)\)\|≥minR∈ℜn,d​\(r,r\)⁡maxx∈\[a,b\]⁡\|1−R​\(x\)\|\\displaystyle\\geq\\min\_\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\\max\_\{x\\in\[a,b\]\}\|1\-R\(x\)\|=\(∗\)​minR∈ℜn,d​\(r,r\)⁡maxx∈\[a/b,1\]⁡\|1−R​\(x\)\|\\displaystyle\\overset\{\(\*\)\}\{=\}\\min\_\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\\max\_\{x\\in\[a/b,1\]\}\|1\-R\(x\)\|≥minR∈ℜn,d​\(r,r\)⁡maxx∈\[lt/ut,1\]⁡\|1−R​\(x\)\|\\displaystyle\\geq\\min\_\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\\max\_\{x\\in\[l\_\{t\}/u\_\{t\},1\]\}\|1\-R\(x\)\|=minR∈ℜn,d​\(r,r\)⁡maxx∈\[lt,ut\]⁡\|1−R​\(x\)\|\\displaystyle=\\min\_\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\\max\_\{x\\in\[l\_\{t\},u\_\{t\}\]\}\|1\-R\(x\)\|=minR∈ℜn,d​\(r,r\)⁡maxx∈\[l,u\]⁡\|1−R​\(RT−1∘⋯∘R1​\(x\)\)\|\\displaystyle=\\min\_\{R\\in\\mathfrak\{R\}\_\{n,d\}\(r,r\)\}\\max\_\{x\\in\[l,u\]\}\|1\-R\(R\_\{T\-1\}\\circ\\cdots\\circ R\_\{1\}\(x\)\)\|=maxx∈\[l,u\]⁡\|1−R∗​\(x\)\|\\displaystyle=\\max\_\{x\\in\[l,u\]\}\|1\-R^\{\*\}\(x\)\|where\(∗\)\(\*\)is justified by the invariance ofℜn,d​\(r,r\)\\mathfrak\{R\}\_\{n,d\}\(r,r\)under rescaling\. ∎

Finally, we are ready for the following

###### Proof of Theorem[2\.7](https://arxiv.org/html/2605.11181#S2.Thmtheorem7)\- Optimality\.

Let us considerRt∈ℜ​\(r,r\)R\_\{t\}\\in\\mathfrak\{R\}\(r,r\)and defineR~t​\(X\):=Rt​\(X​𝒢†\)​𝒢\\widetilde\{R\}\_\{t\}\(X\):=R\_\{t\}\(X\\mathcal\{G\}^\{\\dagger\}\)\\mathcal\{G\}\. It is immediate to verify that, ifRtR\_\{t\}is an optimal approximant off≡Idf\\equiv\\mathrm\{Id\}according to Theorem[E\.5](https://arxiv.org/html/2605.11181#A5.Thmtheorem5), thenR~t\\widetilde\{R\}\_\{t\}is a best approximant of𝒢\\mathcal\{G\}\. More generally, consideringR∗R^\{\*\}as in Theorem[E\.5](https://arxiv.org/html/2605.11181#A5.Thmtheorem5), let us observe that

R∗​\(x\)=R~T​\(R~T−1​\(⋯​R~1​\(x\)\)\)​𝒢†\.\\displaystyle R^\{\*\}\(x\)=\\widetilde\{R\}\_\{T\}\(\\widetilde\{R\}\_\{T\-1\}\(\\cdots\\widetilde\{R\}\_\{1\}\(x\)\)\)\\mathcal\{G\}^\{\\dagger\}\.Then, it holdsOT​\(X\)=R∗​\(X​𝒢†\)​𝒢O\_\{T\}\(X\)=R^\{\*\}\(X\\mathcal\{G\}^\{\\dagger\}\)\\mathcal\{G\}, which implies thatOT​\(G\)O\_\{T\}\(G\)is the best approximant of𝒢\\mathcal\{G\}for everyT∈ℕT\\in\\mathbb\{N\}\. ∎

#### E\.1\.3Approximation Errors

We will now provide an error estimate for the proposed method\.

To start with, let us point out that, callingF12\{\{\}\_\{2\}F\}\_\{1\}the generalized hypergeometric function given by

F12​\(α,β,γ,ξ\)=∑n=0∞\(α\)n​\(β\)nn\!​\(γ\)n​ξn,\{\{\}\_\{2\}F\}\_\{1\}\(\\alpha,\\beta,\\gamma,\\xi\)=\\sum\_\{n=0\}^\{\\infty\}\\frac\{\(\\alpha\)\_\{n\}\(\\beta\)\_\{n\}\}\{n\!\(\\gamma\)\_\{n\}\}\\xi^\{n\},with\(α\)n:=α​\(α\+1\)​⋯​\(α\+n−1\)\(\\alpha\)\_\{n\}:=\\alpha\(\\alpha\+1\)\\cdots\(\\alpha\+n\-1\), by\[Abramowitz and Stegun,[1964](https://arxiv.org/html/2605.11181#bib.bib11), Chapter 15, 15\.1\.8\], it holds

\(1−x\)−12​r=F12​\(12​r,1,1,x\)\.\(1\-x\)^\{\-\\frac\{1\}\{2r\}\}=\{\{\}\_\{2\}F\}\_\{1\}\\left\(\\frac\{1\}\{2r\},1,1,x\\right\)\.Then, we immediately get the following generalisation of\[Kenney and Laub,[1991](https://arxiv.org/html/2605.11181#bib.bib64), Theorem 3\]\.

###### Lemma E\.7\.

Letr∈ℕr\\in\\mathbb\{N\}and considerP​\(x\)∈ℙne​v​e​n​\(r\)P\(x\)\\in\\mathbb\{P\}\_\{n\}^\{even\}\(r\),Q​\(x\)∈ℙde​v​e​n​\(r\)Q\(x\)\\in\\mathbb\{P\}\_\{d\}^\{even\}\(r\)such thath​\(x\):=P​\(x\)Q​\(x\)h\(x\):=\\frac\{P\(x\)\}\{Q\(x\)\}is the\[n/d\]\[n/d\]Padé\-approximant of\(1−x\)−12​r\(1\-x\)^\{\-\\frac\{1\}\{2r\}\}\. Then, for everyt∈ℕt\\in\\mathbb\{N\}, setting

xt\+1=xt​h​\(1−xt2​r\),x\_\{t\+1\}=x\_\{t\}h\(1\-x\_\{t\}^\{2r\}\),we have

\|1−xt2​r\|≤\|1−x02​r\|\(n\+d\+1\)⊤\.\|1\-x\_\{t\}^\{2r\}\|\\leq\|1\-x\_\{0\}^\{2r\}\|^\{\(n\+d\+1\)^\{\\top\}\}\.

Finally, the following approximation result holds \(concluding the proof of Theorem[2\.7](https://arxiv.org/html/2605.11181#S2.Thmtheorem7)\)\.

###### Theorem E\.8\.

LetT≥0T\\geq 0,l∈\(0,1\)l\\in\(0,1\),σmin​\(G\)≥l\\sigma\_\{\\min\}\(G\)\\geq land considerR∗R^\{\*\}as in Theorem[2\.7](https://arxiv.org/html/2605.11181#S2.Thmtheorem7)\. Then,XTX\_\{T\}given by \([23](https://arxiv.org/html/2605.11181#A5.E23)\) satisfies the following estimate bound

‖\(G​G⊤\)−ab​G−OT‖2≤\|1−lb\|\(n\+d\+1\)Tl2​a−bb\.\\\|\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}G\-O\_\{T\}\\\|\_\{2\}\\leq\\frac\{\|1\-l^\{b\}\|^\{\(n\+d\+1\)^\{T\}\}\}\{l^\{\\frac\{2a\-b\}\{b\}\}\}\.

###### Proof\.

First, let us observe that, settingZt:=\(G​G⊤\)ab−12​OtZ\_\{t\}:=\(GG^\{\\top\}\)^\{\\frac\{a\}\{b\}\-\\frac\{1\}\{2\}\}O\_\{t\},t∈\[1,T\]t\\in\[1,T\], it holds

l2​a−bb​‖\(G​G⊤\)−ab​G−OT‖2\\displaystyle l^\{\\frac\{2a\-b\}\{b\}\}\\left\\\|\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}G\-O\_\{T\}\\right\\\|\_\{2\}≤σm​i​n​\(G\)2​a−bb​‖\(G​G⊤\)−ab​G−OT‖2\\displaystyle\\leq\\sigma\_\{min\}\(G\)^\{\\frac\{2a\-b\}\{b\}\}\\left\\\|\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}G\-O\_\{T\}\\right\\\|\_\{2\}≤‖\(G​G⊤\)ab−12​\(\(G​G⊤\)−ab​G−OT\)‖2\\displaystyle\\leq\\left\\\|\(GG^\{\\top\}\)^\{\\frac\{a\}\{b\}\-\\frac\{1\}\{2\}\}\\left\(\(GG^\{\\top\}\)^\{\-\\frac\{a\}\{b\}\}G\-O\_\{T\}\\right\)\\right\\\|\_\{2\}=‖sign​\(G\)−Zt‖2\\displaystyle=\\\|\\mathrm\{sign\}\(G\)\-Z\_\{t\}\\\|\_\{2\}It is fundamental to observe that, by direct computation, we haveXt∈ℜ​\(b2,b2\)X\_\{t\}\\in\\mathfrak\{R\}\\left\(\\frac\{b\}\{2\},\\frac\{b\}\{2\}\\right\)\.

Let nowhhbe the\[n/d\]\[n/d\]Padé\-approximant of\(1−x\)−1b\(1\-x\)^\{\-\\frac\{1\}\{b\}\}andR~​\(x\)=x​h​\(1−xb\)\\widetilde\{R\}\(x\)=xh\(1\-x^\{b\}\)\. Moreover, by direct computation, we can immediately check thath​\(1−x2​r\)=P​\(xb\)Q​\(xb\)h\(1\-x^\{2r\}\)=\\frac\{P\(x^\{b\}\)\}\{Q\(x^\{b\}\)\}for someP∈ℙnP\\in\\mathbb\{P\}\_\{n\},Q∈ℙdQ\\in\\mathbb\{P\}\_\{d\}\. Thus,R~∈ℜn,d​\(b2,b2\)\\widetilde\{R\}\\in\\mathfrak\{R\}\_\{n,d\}\(\\frac\{b\}\{2\},\\frac\{b\}\{2\}\), and, by Lemma[E\.4](https://arxiv.org/html/2605.11181#A5.Thmtheorem4), it is well defined

f:=R~​∘⋯∘⏟Ttimes​R~∈ℜ​\(b2,b2\)\.f:=\\widetilde\{R\}\\;\\underbrace\{\\circ\\cdots\\circ\}\_\{\\text\{$T$ times\}\}\\;\\widetilde\{R\}\\in\\mathfrak\{R\}\(\\frac\{b\}\{2\},\\frac\{b\}\{2\}\)\.Finally, consideringR∗R^\{\*\}as in Theorem[E\.5](https://arxiv.org/html/2605.11181#A5.Thmtheorem5), by Lemma[E\.7](https://arxiv.org/html/2605.11181#A5.Thmtheorem7), we have

‖sign​\(G\)−ZT‖2\\displaystyle\\\|\\mathrm\{sign\}\(G\)\-Z\_\{T\}\\\|\_\{2\}≤maxx∈\[l,1\]⁡\|1−R⋆​\(x\)\|≤maxx∈\[l,1\]⁡\|1−f​\(x\)\|\\displaystyle\\leq\\max\_\{x\\in\[l,1\]\}\|1\-R^\{\\star\}\(x\)\|\\leq\\max\_\{x\\in\[l,1\]\}\|1\-f\(x\)\|≤maxx∈\[l,1\]⁡\|1−f​\(x\)b\|\(n\+1\)T∑i=0b−1f​\(x\)i≤maxx∈\[l,1\]⁡\|1−xb\|\(n\+d\+1\)T∑i=0b−1f​\(x\)i\\displaystyle\\leq\\max\_\{x\\in\[l,1\]\}\\frac\{\|1\-f\(x\)^\{b\}\|^\{\(n\+1\)^\{T\}\}\}\{\\sum\_\{i=0\}^\{b\-1\}f\(x\)^\{i\}\}\\leq\\max\_\{x\\in\[l,1\]\}\\frac\{\|1\-x^\{b\}\|^\{\(n\+d\+1\)^\{T\}\}\}\{\\sum\_\{i=0\}^\{b\-1\}f\(x\)^\{i\}\}≤\|1−lb\|\(n\+d\+1\)T,\\displaystyle\\leq\|1\-l^\{b\}\|^\{\(n\+d\+1\)^\{T\}\},and we can conclude\. ∎

### E\.2Accuracy of approximations

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_2_3_cond2.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_3_4_cond2.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_1_1_cond2.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_2_3_cond4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_3_4_cond4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_1_1_cond4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_2_3_cond6.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_3_4_cond6.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/sv_comparison_1_1_cond6.png)Figure E\.1:Output singular values of three methods for computing\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}Gon a synthetic256×128256\\times 128matrixG=U​Σ​V⊤G=U\\Sigma V^\{\\top\}, where the singular values ofΣ\\Sigmaare spaced logarithmically between11andκ−1\\kappa^\{\-1\}to give condition numberκ\\kappa\. Columns correspond to exponentsa/b∈\{2/3,3/4,1\}a/b\\in\\\{2/3,\\,3/4,\\,1\\\}; rows to condition numbersκ∈\{102,104,106\}\\kappa\\in\\\{10^\{2\},\\,10^\{4\},\\,10^\{6\}\\\}\. Each panel shows three sub\-panels for float16, float32, and float64 \(left to right\)\. Within each sub\-panel the solid black line is the exact targetσi​\(G\)1−2​a/b\\sigma\_\{i\}\(G\)^\{1\-2a/b\}\(sorted descending by index\), the dashed black line is the regularized targetσi​\(G\)​\(σi​\(G\)2\+ε\)−a/b\\sigma\_\{i\}\(G\)\\,\(\\sigma\_\{i\}\(G\)^\{2\}\+\\varepsilon\)^\{\-a/b\}whereε\\varepsilonis the dtype\-dependent regularization threshold scaled by the spectral norm ofGG, and the coloured curves are the output singular values ofdirect\(red, polynomial Newton–Schulz applied directly toGG\),coupled\_pq\(blue, coupled polynomial iteration trackingGGand its Gram matrix\), andcoupled\_chol\(green, rational Newton–Schulz via Cholesky tracking\)\. All methods use 10 iterations\. A method that diverges \(NaN/Inf output\) is shown as a flat phantom line in the legend only\.coupled\_cholremains stable across all conditions and precisions, whiledirectandcoupled\_pqdegrade or diverge at high condition numbers in low precision\.
### E\.3Instability of Polynomial Iterations

[Figure˜E\.2](https://arxiv.org/html/2605.11181#A5.F2)compares three iterative methods for computing\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}Gacross matrix sizes, condition numbers, and floating\-point precisions\. Grey cells indicate divergence \(NaN or Inf output\)\. Bothdirectandcoupled\_pqexhibit systematic divergence at high condition numbers, particularly infloat16andbfloat16, with divergence spreading tofloat32at extreme conditioning\. Regularization delays but does not eliminate this behaviour\. In contrast,coupled\_chol, remains numerically stable across all tested conditions and dtypes, achieving near\-machine\-precision accuracy throughout\. This stability comes without any regularization: the Cholesky structure enforces positive definiteness at every step, preventing the eigenvalue collapse that causes the other methods to diverge\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/heatmaps_direct_sv_ratio_err_noreg.png)

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/heatmaps_coupled_pq_sv_ratio_err_noreg.png)

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/heatmaps_coupled_chol_sv_ratio_err_noreg.png)

Figure E\.2:Singular\-value ratio errorεsv=maxi⁡\|σi​\(Z^\)⋅σi​\(G\)\(2​a−b\)/b−1\|\\varepsilon\_\{\\mathrm\{sv\}\}=\\max\_\{i\}\|\\sigma\_\{i\}\(\\hat\{Z\}\)\\cdot\\sigma\_\{i\}\(G\)^\{\(2a\-b\)/b\}\-1\|without regularization for three iterative methods computing\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}G:direct\(left\),coupled\_pq\(centre\), andcoupled\_chol\(right\)\. Each panel is a4×34\\times 3grid of heatmaps indexed by exponenta/b∈1/1,1/2,3/4,2/3a/b\\in\{1/1,1/2,3/4,2/3\}\(rows\) and matrix size\(m,n\)∈64×32,256×128,512×256\(m,n\)\\in\{64\\times 32,256\\times 128,512\\times 256\}\(columns\); within each cell theyy\-axis encodes condition numberκ∈101,102,…,1012,1016\\kappa\\in\{10^\{1\},10^\{2\},\\dots,10^\{12\},10^\{16\}\}\(increasing downward\) and thexx\-axis encodes floating\-point type \(f16, bf16, f32, f64, increasing precision left to right\)\. Test matrices are constructed asG=U​Σ​V⊤G=U\\Sigma V^\{\\top\}, whereU∈ℝm×mU\\in\\mathbb\{R\}^\{m\\times m\}andV∈ℝn×nV\\in\\mathbb\{R\}^\{n\\times n\}are independent Haar\-distributed orthogonal matrices obtained by QR\-decomposing Gaussian random matrices infloat64, andΣ=diag⁡\(σ1,…,σk\)\\Sigma=\\operatorname\{diag\}\(\\sigma\_\{1\},\\dots,\\sigma\_\{k\}\)withk=min⁡\(m,n\)k=\\min\(m,n\)singular values spaced logarithmically from11toκ−1\\kappa^\{\-1\}; each matrix is then cast to the target dtype before the iterative method is applied\. The error metric equals zero for a perfect computation and is invariant to the overall scale ofGG; no regularization is applied, so grey cells indicate that the method produced non\-finite values for that configuration\.directandcoupled\_pqeach perform 40 steps of a degree\-\(p,q\)\(p,q\)polynomial iteration whose coefficients are precomputed by Zolotarev\-optimal minimax approximation on\[ℓ,1\]\[\\ell,1\]withℓ=10−16\\ell=10^\{\-16\}, run entirely in the target dtype\.coupled\_cholperforms 25 steps of a rational Remez iteration with Cholesky\-based stabilization\. Colour encodes error on a logarithmic scale \(viridis\); darker shades indicate smaller error\. Error grows with condition number and worsens at reduced precision, whilecoupled\_pqachieves consistently lower error thandirectacross all settings\. Most notablycoupled\_cholremains stable at all condition numbers and precisions\.
### E\.4Iterative Methods for\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}\\,G

Throughout this section, letG∈ℝm×nG\\in\\mathbb\{R\}^\{m\\times n\}\(m≤nm\\leq n\) denote the \(Frobenius\-normalised\) gradient matrix\. LetPtP\_\{t\}be a step\-varying degree\-2 polynomial with Remez\-approximated coefficients, and letk∈ℤ\+k\\in\\mathbb\{Z\}^\{\+\}satisfyP​\(x\)k​xk​a/b≈1P\(x\)^\{k\}x^\{ka/b\}\\approx 1on the fitting interval\. All methods target the fixed pointX∞=\(G​G⊤\)−a/b​GX\_\{\\infty\}=\(GG^\{\\top\}\)^\{\-a/b\}\\,G\.

A critical note on numerical stability:It is important to state upfront that, as identified by HighamHigham \[[2008](https://arxiv.org/html/2605.11181#bib.bib36)\], and numerically illustrated in[Section˜E\.3](https://arxiv.org/html/2605.11181#A5.SS3), most of the polynomial and coupled iterations detailed below fail to converge in finite\-precision arithmetic due to severe numerical stability issues \(e\.g\., condition\-number squaring and diverging coupled sequences\)\. In practice,*only the final method presented*\(theCoupled\-Chol Iteration\) exhibits the forward stability required to reliably converge\. The preceding methods are documented primarily for their theoretical connections and historical context\.

For cost accounting, we distinguish*G\-MMs*\(fullm×nm\\times ngradient products,O​\(m2​n\)O\(m^\{2\}n\)\) from*S\-MMs*\(m×mm\\times msquare products,O​\(m3\)O\(m^\{3\}\)\) and*QRs*\(thin QR of a2​m×m2m\\times mblock,O​\(m3\)O\(m^\{3\}\)\)\. Integer powersAkA^\{k\}require⌈log2⁡k⌉\\lceil\\log\_\{2\}k\\rceilS\-MMs via binary exponentiation\.

#### E\.4\.1Direct Iteration

The direct method iteratesOtO\_\{t\}without maintaining separate state for the polynomial argumentAtA\_\{t\}\. Instead,AtA\_\{t\}is recomputed at each step from the current iterate\. Three cases arise depending ona/ba/b:

##### General case \(a/b∉\{1/2,1\}a/b\\notin\\\{1/2,1\\\}\)\.

A frozen Gram power\(G​G⊤\)k​\(a/b−1/2\)\(GG^\{\\top\}\)^\{k\(a/b\-1/2\)\}is precomputed once; each step forms:

\{Ot\+1=Rt​\(At\)​Ot,O0=G,At=\(Ot​Ot⊤\)k/2​\(G​G⊤\)k​\(a/b−1/2\)\.\\begin\{cases\}O\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] A\_\{t\}=\(O\_\{t\}O\_\{t\}^\{\\top\}\)^\{k/2\}\\,\(GG^\{\\top\}\)^\{k\(a/b\\,\-\\,1/2\)\}\.\\end\{cases\}\(21\)

##### Polar case \(a/b=1/2a/b=1/2\)\.

The frozen factor vanishes \(k​\(a/b−1/2\)=0k\(a/b\-1/2\)=0\) andk/2=1k/2=1, simplifying the argument to the plain Gram matrix of the current iterate:

\{Ot\+1=Rt​\(At\)​Ot,O0=G,At=Ot​Ot⊤\.\\begin\{cases\}O\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] A\_\{t\}=O\_\{t\}O\_\{t\}^\{\\top\}\.\\end\{cases\}\(22\)

##### Freon case \(a/b=1a/b=1,k=1k=1\)\.

The general formula \([21](https://arxiv.org/html/2605.11181#A5.E21)\) would require an expensive matrix square root\(Xt​Xt⊤\)1/2\(X\_\{t\}X\_\{t\}^\{\\top\}\)^\{1/2\}at every step\. Instead,AtA\_\{t\}is defined via a mixed product with the*frozen*initial transposeG⊤G^\{\\top\}:

\{Ot\+1=Rt​\(At\)​Ot,O0=G,At=Ot​G0⊤\.\\begin\{cases\}O\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] A\_\{t\}=O\_\{t\}\\,G\_\{0\}^\{\\top\}\.\\end\{cases\}\(23\)UsingG0⊤G\_\{0\}^\{\\top\}avoids both the matrix square root and condition\-number squaring, sinceκ​\(Ot​G⊤\)≤κ​\(Ot\)​κ​\(G\)\\kappa\(O\_\{t\}G^\{\\top\}\)\\leq\\kappa\(O\_\{t\}\)\\kappa\(G\)compared toκ​\(Ot​Ot⊤\)=κ​\(Ot\)2\\kappa\(O\_\{t\}O\_\{t\}^\{\\top\}\)=\\kappa\(O\_\{t\}\)^\{2\}\.

In all three cases, convergence toO∞​O∞⊤=\(G​G⊤\)1−2​a/bO\_\{\\infty\}O\_\{\\infty\}^\{\\top\}=\(GG^\{\\top\}\)^\{1\-2a/b\}ensuresA∞=IA\_\{\\infty\}=I\(the fixed point ofPP\)\.

##### Classical connection\.

Fora/b=1/2a/b=1/2and degree\-1PtP\_\{t\}, the update reduces to the Newton–Schulz iteration for the polar decompositionHigham \[[2008](https://arxiv.org/html/2605.11181#bib.bib36)\]\. The higher power inside the polynomial allows for generala/b≠1/2a/b\\neq 1/2\.

##### Complexity\.

- •Init:1 G\-MM;⌈log2⁡\(k​\(a/b−1/2\)\)⌉\\lceil\\log\_\{2\}\(k\(a/b\-1/2\)\)\\rceilS\-MMs \(zero ifa/b=1/2a/b=1/2\)\.
- •Per step:2 G\-MMs;\(1\+⌈log2⁡\(k/2\)⌉\+𝟏\[a/b≠1/2\]\)\(1\+\\lceil\\log\_\{2\}\(k/2\)\\rceil\+\\mathbf\{1\}\_\{\[a/b\\neq 1/2\]\}\)S\-MMs\.

#### E\.4\.2M\-Accumulator Iteration

Instead of iteratingOtO\_\{t\}directly, this method maintains a preconditionerMt∈ℝm×mM\_\{t\}\\in\\mathbb\{R\}^\{m\\times m\}, applied to the frozen gradientGGonly upon convergence:

\{Mt\+1=Rt​\(At\)​Mt,M0=I,At=Mtk​\(G​G⊤\)k​a/b\.\\begin\{cases\}M\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)\\,M\_\{t\},&M\_\{0\}=I,\\\\\[4\.0pt\] A\_\{t\}=M\_\{t\}^\{k\}\\,\(GG^\{\\top\}\)^\{ka/b\}\.\\end\{cases\}\(24\)The final result isOT=MT​GO\_\{T\}=M\_\{T\}G\. At convergence,M∞=\(G​G⊤\)−a/bM\_\{\\infty\}=\(GG^\{\\top\}\)^\{\-a/b\}andA∞=IA\_\{\\infty\}=I\.

##### Classical connection\.

Fora/b=1/2a/b=1/2, this is the uncoupled polynomial analogue of the Denman–Beavers iteration for the matrix square rootHigham \[[2008](https://arxiv.org/html/2605.11181#bib.bib36)\]\. This is also the method of\[Zhanget al\.,[2026](https://arxiv.org/html/2605.11181#bib.bib35)\]\.

##### Complexity\.

- •Init:1 G\-MM;⌈log2⁡\(k​a/b\)⌉\\lceil\\log\_\{2\}\(ka/b\)\\rceilS\-MMs\.
- •Per step:0 G\-MMs;\(⌈log2⁡k⌉\+3\)\(\\lceil\\log\_\{2\}k\\rceil\+3\)S\-MMs\.
- •Final:1 G\-MM \(MT​G0M\_\{T\}G\_\{0\}\)\.

M\-accumulator is the cheapest per\-step method whenn≫mn\\gg m, deferring the large\-matrix product to the end\.

#### E\.4\.3Coupled Iteration

The coupled method eliminates the per\-step recomputation ofAtA\_\{t\}by tracking it as its own state\. BecausePt​\(At\)P\_\{t\}\(A\_\{t\}\)commutes withAtA\_\{t\}, the system closes exactly:

\{Ot\+1=Rt​\(At\)​Ot,O0=G,At\+1=Rt​\(At\)k​At,A0=\(G​G⊤\)k​a/b\.\\begin\{cases\}O\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] A\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)^\{k\}\\,A\_\{t\},&A\_\{0\}=\(GG^\{\\top\}\)^\{ka/b\}\.\\end\{cases\}\(25\)
##### Classical connection\.

Tracking the argument matrix alongside the iterand echoes Schulz\-type coupled iterations for matrix inversionHigham \[[2008](https://arxiv.org/html/2605.11181#bib.bib36)\]\.

##### Complexity\.

- •Init:1 G\-MM;⌈log2⁡\(k​a/b\)⌉\\lceil\\log\_\{2\}\(ka/b\)\\rceilS\-MMs\.
- •Per step:1 G\-MM;\(⌈log2⁡k⌉\+2\)\(\\lceil\\log\_\{2\}k\\rceil\+2\)S\-MMs\.

Compared to M\-accumulator, this adds one G\-MM per step but saves one S\-MM, making it favorable only whenm≈nm\\approx n\.

#### E\.4\.4Coupled\-a​babIteration

A variant of coupled iteration whereAtA\_\{t\}initializes at the raw Gram matrixG​G⊤GG^\{\\top\}, avoiding the fractional matrix power during initialization:

\{Ot\+1=Rt​\(At\)a​Ot,O0=G,At\+1=Rt​\(At\)b​At,A0=G​G⊤\.\\begin\{cases\}O\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)^\{a\}\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] A\_\{t\+1\}=R\_\{t\}\(A\_\{t\}\)^\{b\}\\,A\_\{t\},&A\_\{0\}=GG^\{\\top\}\.\\end\{cases\}\(26\)The invariantR​\(x\)b​x≈1R\(x\)^\{b\}x\\approx 1ensuresAt→IA\_\{t\}\\to IwhileOt→\(G​G⊤\)−a/b​GO\_\{t\}\\to\(GG^\{\\top\}\)^\{\-a/b\}G\.

##### Complexity\.

- •Init:1 G\-MM; 0 S\-MMs\.
- •Per step:1 G\-MM;\(2\+⌈log2⁡a⌉\+⌈log2⁡b⌉\)\(2\+\\lceil\\log\_\{2\}a\\rceil\+\\lceil\\log\_\{2\}b\\rceil\)S\-MMs\.

#### E\.4\.5Dual\-a​babIteration

This method maintains a factored representationAt=Xt​YtA\_\{t\}=X\_\{t\}Y\_\{t\}to avoid formingG​G⊤GG^\{\\top\}explicitly at each step:

\{At=Ot​Yt,Rt=Rt​\(At\),Ot\+1=Rta​Ot,O0=G,Yt\+1=Yt​Rtb−a,Y0=G⊤\.\\begin\{cases\}A\_\{t\}=O\_\{t\}Y\_\{t\},\\quad R\_\{t\}=R\_\{t\}\(A\_\{t\}\),\\\\\[4\.0pt\] O\_\{t\+1\}=R\_\{t\}^\{a\}\\,O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] Y\_\{t\+1\}=Y\_\{t\}\\,R\_\{t\}^\{b\-a\},&Y\_\{0\}=G^\{\\top\}\.\\end\{cases\}\(27\)A balancing stepOt←γt​OtO\_\{t\}\\leftarrow\\gamma\_\{t\}O\_\{t\},Yt←Yt/γtY\_\{t\}\\leftarrow Y\_\{t\}/\\gamma\_\{t\}withγt=‖Yt‖F/‖Xt‖F\\gamma\_\{t\}=\\sqrt\{\\\|Y\_\{t\}\\\|\_\{F\}/\\\|X\_\{t\}\\\|\_\{F\}\}prevents mismatched sequence growth\.

##### Classical connection\.

The factored form averts condition\-number squaring, mirroring product formulas for evaluating matrix functions of productsHigham \[[2008](https://arxiv.org/html/2605.11181#bib.bib36)\]\.

##### Complexity\.

- •Init:0 G\-MMs \(transpose only\)\.
- •Per step:3 G\-MMs;\(1\+⌈log2⁡a⌉\+⌈log2⁡\(b−a\)⌉\)\(1\+\\lceil\\log\_\{2\}a\\rceil\+\\lceil\\log\_\{2\}\(b\-a\)\\rceil\)S\-MMs\.

#### E\.4\.6Coupled\-Dual Iteration

This method runs two coupled masters in parallel for exponentsa/ba/band1−a/b1\-a/b\. A cross\-stabilizerEt=I−Xt​Yt→0E\_\{t\}=I\-X\_\{t\}Y\_\{t\}\\to 0corrects the slave iterates directly, avoiding fractional matrix powers inside the loop:

\{RA=Rt​\(At\),RB=Rt​\(Bt\),Et=I−Ot​Yt,Ot\+1=RA​Ot\+γ​Et​Ot,O0=G,Yt\+1=Yt​RB\+γ​Yt​Et,Y0=G⊤,At\+1=RAk​At,A0=\(G​G⊤\)k​a/b,Bt\+1=RBk​Bt,B0=\(G​G⊤\)k​\(b−a\)/b,\\begin\{cases\}R\_\{A\}=R\_\{t\}\(A\_\{t\}\),\\quad R\_\{B\}=R\_\{t\}\(B\_\{t\}\),\\quad E\_\{t\}=I\-O\_\{t\}Y\_\{t\},\\\\\[4\.0pt\] O\_\{t\+1\}=R\_\{A\}\\,O\_\{t\}\+\\gamma\\,E\_\{t\}O\_\{t\},&O\_\{0\}=G,\\\\\[4\.0pt\] Y\_\{t\+1\}=Y\_\{t\}\\,R\_\{B\}\+\\gamma\\,Y\_\{t\}E\_\{t\},&Y\_\{0\}=G^\{\\top\},\\\\\[4\.0pt\] A\_\{t\+1\}=R\_\{A\}^\{k\}\\,A\_\{t\},&A\_\{0\}=\(GG^\{\\top\}\)^\{ka/b\},\\\\\[4\.0pt\] B\_\{t\+1\}=R\_\{B\}^\{k\}\\,B\_\{t\},&B\_\{0\}=\(GG^\{\\top\}\)^\{k\(b\-a\)/b\},\\end\{cases\}\(28\)whereγ\>0\\gamma\>0is a scalar hyperparameter\.

##### Complexity\.

- •Init:1 G\-MM;⌈log2⁡\(k​a/b\)⌉\+⌈log2⁡\(k​\(b−a\)/b\)⌉\\lceil\\log\_\{2\}\(ka/b\)\\rceil\+\\lceil\\log\_\{2\}\(k\(b\-a\)/b\)\\rceilS\-MMs\.
- •Per step:3 G\-MMs;\(4\+2​⌈log2⁡k⌉\)\(4\+2\\lceil\\log\_\{2\}k\\rceil\)S\-MMs\.

#### E\.4\.7Rational\-Chol \(QDWH\-style\) Iteration

Instead of polynomials, this method applies a rational functionRt​\(x\)=\(αt\+βt​x\)/\(1\+δt​x\)R\_\{t\}\(x\)=\(\\alpha\_\{t\}\+\\beta\_\{t\}x\)/\(1\+\\delta\_\{t\}x\)targetingx−1/bx^\{\-1/b\}\. It maintains a preconditionerMt→\(G​G⊤\)−1/bM\_\{t\}\\to\(GG^\{\\top\}\)^\{\-1/b\}:

Mt\+1=Rt​\(Mtb​G​G⊤\)​Mt,M0=I\.M\_\{t\+1\}=R\_\{t\}\(M\_\{t\}^\{b\}\\,GG^\{\\top\}\)\\,M\_\{t\},\\quad M\_\{0\}=I\.Writingb=2​rb=2rand taking a thin QRL​L⊤=G​G⊤LL^\{\\top\}=GG^\{\\top\}, we defineYt=Mtr​LY\_\{t\}=M\_\{t\}^\{r\}L\. A block QR avoids forming\(I\+δt​Yt​Yt⊤\)−1\(I\+\\delta\_\{t\}Y\_\{t\}Y\_\{t\}^\{\\top\}\)^\{\-1\}:

\{Kt=\[δt​Yt⊤Im\],\[Q1;Q2\]​R=Kt​\(thin QR\),Mt\+1=γt​Mt\+\(αt−γt\)​Q2​Q2⊤​Mt,\\begin\{cases\}K\_\{t\}=\\begin\{bmatrix\}\\sqrt\{\\delta\_\{t\}\}\\,Y\_\{t\}^\{\\top\}\\\\ I\_\{m\}\\end\{bmatrix\},\\quad\[Q\_\{1\};\\,Q\_\{2\}\]\\,R=K\_\{t\}\\;\\text\{\(thin QR\)\},\\\\\[6\.0pt\] M\_\{t\+1\}=\\gamma\_\{t\}\\,M\_\{t\}\+\(\\alpha\_\{t\}\-\\gamma\_\{t\}\)\\,Q\_\{2\}\\,Q\_\{2\}^\{\\top\}M\_\{t\},\\end\{cases\}\(29\)whereγt=βt/δt\\gamma\_\{t\}=\\beta\_\{t\}/\\delta\_\{t\}\. The final result isMTa​GM\_\{T\}^\{a\}\\,G\.

##### Classical connection\.

Fora/b=1/2a/b=1/2, this is the QDWH algorithm for polar decomposition, using Zolotarev approximants via block QRs to avoid condition\-number squaringNakatsukasa and Freund \[[2016](https://arxiv.org/html/2605.11181#bib.bib8)\]\.

##### Complexity\.

- •Init:1 thin QR ofG⊤G^\{\\top\}\(equivalent to 1 G\-MM\)\.
- •Per step:1 QR;\(3\+⌈log2⁡r⌉\)\(3\+\\lceil\\log\_\{2\}r\\rceil\)S\-MMs\.
- •Final:⌈log2⁡a⌉\\lceil\\log\_\{2\}a\\rceilS\-MMs; 1 G\-MM \(MTa​GM\_\{T\}^\{a\}G\)\.

#### E\.4\.8Coupled\-Chol Iteration

This method merges Cholesky\-factor tracking \(from Coupled\-a​bab\) with the block QR technique \(from Rational\-Chol\)\. It tracks a Cholesky\-like factorLt∈ℝm×mL\_\{t\}\\in\\mathbb\{R\}^\{m\\times m\}\(Lt​Lt⊤→IL\_\{t\}L\_\{t\}^\{\\top\}\\to I\) directly alongside an accumulatorCt∈ℝm×mC\_\{t\}\\in\\mathbb\{R\}^\{m\\times m\}:

\{Wt=Rt​\(Lt​Lt⊤\)\(via block QR on​Lt​\),Lt\+1=Wtr​Lt,L0=L​\(thin QR factor of​G⊤​\),Ct\+1=Wt​Ct,C0=I\.\\begin\{cases\}W\_\{t\}=R\_\{t\}\(L\_\{t\}L\_\{t\}^\{\\top\}\)&\\text\{\(via block QR on \}L\_\{t\}\\text\{\)\},\\\\\[4\.0pt\] L\_\{t\+1\}=W\_\{t\}^\{r\}\\,L\_\{t\},&L\_\{0\}=L\\;\\text\{\(thin QR factor of \}G^\{\\top\}\\text\{\)\},\\\\\[4\.0pt\] C\_\{t\+1\}=W\_\{t\}\\,C\_\{t\},&C\_\{0\}=I\.\\end\{cases\}\(30\)The final result isCTa​GC\_\{T\}^\{a\}\\,G\.

##### Classical connection\.

This directly generalizes QDWH\. Instead of exponentiating the accumulated preconditionerMtrM\_\{t\}^\{r\}, it exponents only the well\-conditioned, single\-step updateWtrW\_\{t\}^\{r\}, improving forward\-stability\.

##### Complexity\.

- •Init:1 thin QR ofG⊤G^\{\\top\}\.
- •Per step:1 QR;\(3\+⌈log2⁡r⌉\)\(3\+\\lceil\\log\_\{2\}r\\rceil\)S\-MMs\.
- •Final:⌈log2⁡a⌉\\lceil\\log\_\{2\}a\\rceilS\-MMs; 1 G\-MM \(CTa​GC\_\{T\}^\{a\}G\)\.

Costs match Rational\-Chol, but the numerical stability is superior becauseWt≈IW\_\{t\}\\approx I\.

#### E\.4\.9Remez Visualisations

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_poly_composition_1_2.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_rational_composition_1_2.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_poly_composition_3_4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_rational_composition_3_4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_poly_composition_2_3.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_rational_composition_2_3.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_poly_composition_1_4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_rational_composition_1_4.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_poly_composition_1_1.png)![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/paper_rational_composition_1_1.png)Figure E\.3:Convergence of polynomial \(left columns\) and rational \(right columns\) Remez iterations for computing\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}G, shown fora/b∈1/2,3/4,2/3,1/4,1/1a/b\\in\{1/2,3/4,2/3,1/4,1/1\}\. Each pair of panels shares anxx\-axis representing a scalarx∈\[ℓ,1\]x\\in\[\\ell,1\], which stands in for the singular values of the \(normalised\) input matrix after mapping to the unit interval\.*Left panel of each pair*: the raw iterated compositionPT​\(x\)P\_\{T\}\(x\)\(polynomial\) orRT​\(x\)R\_\{T\}\(x\)\(rational\), which is initialised at the identity and should converge to11on\[ℓ,1\]\[\\ell,1\]as the number of stepsTTincreases; curves for successive values ofTTare overlaid, progressing from light to dark\.*Right panel of each pair*: the rescaled quantityx1−2​a/b​PT​\(x\)x^\{1\-2a/b\}P\_\{T\}\(x\)orx1−2​a/b​RT​\(x\)x^\{1\-2a/b\}R\_\{T\}\(x\)plotted against the targetx1−2​a/bx^\{1\-2a/b\}; since the singular values of the exact output\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}Gareσi1−2​a/b\\sigma\_\{i\}^\{1\-2a/b\}, convergence of the composition to11is equivalent to convergence of the rescaled curve to the target\. Polynomial coefficients are obtained obtained by the Remez algorithm on\[ℓ,1\]\[\\ell,1\]withℓ=10−16\\ell=10^\{\-16\}; rational coefficients are obtained by the Remez algorithm withℓ=10−5\\ell=10^\{\-5\}\. For rational approximations, exponentsa/ba/bwith odd denominator \(e\.g\.2/32/3\), the rational iteration internally doubles to4/64/6to ensure the companion\-form construction remains valid\.

### E\.5Picking Parameters for Freon Iterations: Cushion and l

The rational Newton–Schulz iteration requires two hyperparameters: a cushion factorδ\\deltathat stabilises the Remez denominator coefficients, and a lower spectral boundllthat determines how wide a range the iteration must cover\. Both must be set as functions ofbbto avoid either numerical blow\-up or wasted iterations\. Given the iterations wildly oscillate untilllgets mapped close to11, which are then stabilizedllis chosen such that the 5th iteration stabilizes\. Cushion primarily controls the size of the resulting coefficients out of the fitting, with smaller cushions resulting in ever larger values ofγ\\gamma\. For this reason we introduce an artificial numerics based bound of10510^\{5\}\. From empirical sweeps \([Figure˜E\.4](https://arxiv.org/html/2605.11181#A5.F4)\) we derive closed\-form schedulesδ​\(b\)=\(1\.84×10−8\)1/b\\delta\(b\)=\(1\.84\\times 10^\{\-8\}\)^\{1/b\}andl​\(b\)=\(10−11\)2/bl\(b\)=\(10^\{\-11\}\)^\{2/b\}, which we use throughout without further tuning\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/x7.png)

![Refer to caption](https://arxiv.org/html/2605.11181v1/x8.png)

Figure E\.4:Hyperparameter selection for the rational Newton–Schulz iteration\.\(Left\)Minimum cushion factor required to keep the denominator coefficients bounded \(maxt⁡γt≤105\\max\_\{t\}\\gamma\_\{t\}\\leq 10^\{5\}\) as a function of effectiveqq\(where oddqqis replaced by2​q2q\)\. The fitted formula\(1\.84×10−8\)1/q\(1\.84\\times 10^\{\-8\}\)^\{1/q\}safely exceeds the critical threshold for allqqand replaces the old fixed value of0\.0240\.024\.\(Right\)Smallest lower boundllsuch that the iteration fully converges within five Newton–Schulz steps, i\.e\. the final interval width reaches its asymptotic fixed point\. The formula\(10−11\)2/q\(10^\{\-11\}\)^\{2/q\}matches the empirical optimum across all tested values ofqqand ensures all steps contribute to convergence\.Input:Gradient

G∈ℝk×nG\\in\\mathbb\{R\}^\{k\\times n\}\(

k≤nk\\leq n\), powers

a,b∈ℕa,b\\in\\mathbb\{N\}, rational coefficients

\{\(αt,βt,γt\)\}t=1N\\\{\(\\alpha\_\{t\},\\beta\_\{t\},\\gamma\_\{t\}\)\\\}\_\{t=1\}^\{N\}, number of steps

NN, regularization

ϵ\\epsilon\.

Output:Approximation to

\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}G\.

//Ensurebbis even to yield an integerrrfor the block QR step

if*b\(mod2\)≠0b\\pmod\{2\}\\neq 0*then

a←2​a,b←2​ba\\leftarrow 2a,\\quad b\\leftarrow 2b
end if

r←b/2r\\leftarrow b/2

//Normalize gradient to improve numerical stability

ν←‖G‖F\+ϵ\\nu\\leftarrow\\\|G\\\|\_\{F\}\+\\epsilon

G←G/νG\\leftarrow G/\\nu

//InitializeLLvia thin QR of augmented gradient

Q0,R0←QR​\(\[G⊤ϵ​Ik\]\)Q\_\{0\},R\_\{0\}\\leftarrow\\text\{QR\}\\left\(\\begin\{bmatrix\}G^\{\\top\}\\\\ \\sqrt\{\\epsilon\}I\_\{k\}\\end\{bmatrix\}\\right\)

L←R0⊤L\\leftarrow R\_\{0\}^\{\\top\}

//YieldsL​L⊤=G​G⊤\+ϵ​ILL^\{\\top\}=GG^\{\\top\}\+\\epsilon I

L0←LL\_\{0\}\\leftarrow L

//Saved for log\-determinant;logdet\(GG⊤\)=2∑ilog\[L0\]i​i\+2klogν\\log\\det\(GG^\{\\top\}\)=2\\sum\_\{i\}\\log\[L\_\{0\}\]\_\{ii\}\+2k\\log\\nu

C←IkC\\leftarrow I\_\{k\}

//Accumulator

for*t=1t=1toNN*do

//ComputeW=αt​I\+βt​At​\(I\+γt​At\)−1W=\\alpha\_\{t\}I\+\\beta\_\{t\}A\_\{t\}\(I\+\\gamma\_\{t\}A\_\{t\}\)^\{\-1\}whereAt=L​L⊤A\_\{t\}=LL^\{\\top\}

ρt←βt/γt\\rho\_\{t\}\\leftarrow\\beta\_\{t\}/\\gamma\_\{t\}
//Construct Block MatrixKKoptimally to avoid numeric blowup

if*γ≤1\.0\\gamma\\leq 1\.0*then

K←\[γt​L⊤Ik\]K\\leftarrow\\begin\{bmatrix\}\\sqrt\{\\gamma\_\{t\}\}L^\{\\top\}\\\\ I\_\{k\}\\end\{bmatrix\}
else

K←\[L⊤1γt​Ik\]K\\leftarrow\\begin\{bmatrix\}L^\{\\top\}\\\\ \\frac\{1\}\{\\sqrt\{\\gamma\_\{t\}\}\}I\_\{k\}\\end\{bmatrix\}
end if

QK,RK←QR​\(K\)Q\_\{K\},R\_\{K\}\\leftarrow\\text\{QR\}\(K\)
Q2←bottom​k×k​block of​QKQ\_\{2\}\\leftarrow\\text\{bottom \}k\\times k\\text\{ block of \}Q\_\{K\}
V←Q2​Q2⊤V\\leftarrow Q\_\{2\}Q\_\{2\}^\{\\top\}
//Identity:V=\(I\+γt​L​L⊤\)−1V=\(I\+\\gamma\_\{t\}LL^\{\\top\}\)^\{\-1\}

W←ρt​I\+\(αt−ρt\)​VW\\leftarrow\\rho\_\{t\}I\+\(\\alpha\_\{t\}\-\\rho\_\{t\}\)V
W←12​\(W\+W⊤\)W\\leftarrow\\frac\{1\}\{2\}\(W\+W^\{\\top\}\)
//Ensure numerical symmetry

//Coupled state updates

L←Wr​LL\\leftarrow W^\{r\}L
//Equivalent to:Lt\+1​Lt\+1⊤=Wb​AtL\_\{t\+1\}L\_\{t\+1\}^\{\\top\}=W^\{b\}A\_\{t\}

C←W​CC\\leftarrow WC
//Accumulate product∏Wt\\prod W\_\{t\}

end for

//Rescale and apply to gradient

return*ν1−2​a/b​Ca​G\\nu^\{1\-2a/b\}C^\{a\}G,L0L\_\{0\}*

Algorithm 4Coupled QDWH Cholesky Iteration for\(G​G⊤\)−a/b​G\(GG^\{\\top\}\)^\{\-a/b\}G

## Appendix FFurther Optimizer Details

### F\.1Kaon

![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/pdf_grid.png)\(a\)Empirical stationary distributions of the Kaon mapxn\+1=λ​xn​\(1−xn2\)2x\_\{n\+1\}=\\lambda\\,x\_\{n\}\(1\-x\_\{n\}^\{2\}\)^\{2\}for ten values ofλ\\lambdauniformly spaced in\[3\.5,4\.25\]\[3\.5,4\.25\]\. Each histogram is obtained by evolving5,0005\{,\}000particles for 500 burn\-in iterations followed by 200 collection steps\.
![Refer to caption](https://arxiv.org/html/2605.11181v1/Freon/figures/cobweb_c4_1.png)\(b\)Cobweb diagram of the Kaon map atλ=4\.1\\lambda=4\.1, showing three trajectories of 40 iterations each from independent random initial conditions in\(0\.1,0\.9\)\(0\.1,\\,0\.9\)\. The map curvey=f​\(x\)y=f\(x\)\(blue\) and the diagonaly=xy=x\(dashed\) are overlaid; the chaotic wandering of the iterates is evident\.

Figure F\.1:Dynamics of the Kaon polynomial mapf​\(x\)=λ​x​\(1−x2\)2f\(x\)=\\lambda\\,x\(1\-x^\{2\}\)^\{2\}used in the Kaon optimiser as a Newton\-Schulz\-style spectral normalisation step\.##### Stationary Distribution

The Kaon optimiser applies the polynomial mapf​\(x\)=λ​x​\(1−x2\)2f\(x\)=\\lambda x\(1\-x^\{2\}\)^\{2\}iteratively to the singular values of the gradient matrix as a cheap alternative to the exact polar factorisation used by Muon\. Unlike Newton–Schulz iterations, which converge the singular values to±1\\pm 1, the Kaon map atλ=4\.1\\lambda=4\.1operates in a chaotic regime: rather than converging, the iterates explore a broad stationary distribution supported on\(0,1\.175\)\(0,1\.175\)\.[Figure˜1\(a\)](https://arxiv.org/html/2605.11181#A6.F1.sf1)displays the empirical stationary distributions for ten values ofλ\\lambdaspanning\[3\.5,4\.25\]\[3\.5,4\.25\], estimated by running5,0005\{,\}000particles through 500 burn\-in followed by 200 collection steps of the map\. The distributions shift and spread asλ\\lambdaincreases, with the support and shape varying considerably across the range\. Figure[1\(b\)](https://arxiv.org/html/2605.11181#A6.F1.sf2)shows a cobweb diagram atλ=4\.1\\lambda=4\.1, confirming the chaotic character of the dynamics: trajectories starting from different initial conditions in\(0\.1,0\.9\)\(0\.1,0\.9\)neither converge to a fixed point nor settle into a periodic orbit, but instead wander erratically through the attractor\. This stochasticity in the effective preconditioner is an intrinsic feature of Kaon rather than a deficiency\.

##### Cost of Kaon

Kaonuses the same number of stepsTTasMuon’s Newton–Schulz iteration and matches its matrix\-multiplication cost\. Each step of both methods formsX​X⊤XX^\{\\top\}, squares oner×rr\\times rmatrix, and multiplies the resulting degree\-two matrix polynomial byXX, costingT​\(4​r2​s\+2​r3\)T\(4r^\{2\}s\+2r^\{3\}\)FLOPs to leading order\.

## Appendix GRandom Feature Asymptotics

In this section we show how to compute our two quantities in this asymptotic framework\.

###### Definition G\.1\.

Theiith normalized trace is defined by

τi=1n​tr⁡\(Ci\)\\tau\_\{i\}=\\frac\{1\}\{n\}\\operatorname\{tr\}\(C^\{i\}\)

###### Theorem G\.2\.

Assume thatA=C1/2​ZA=C^\{1/2\}ZwithCCsymmetric nonnegative definite withlimsup‖C‖o​p<∞\\lim\\sup\\\|C\\\|\_\{op\}<\\infty, andZZa matrix with entries which are iid Gaussians\.

LetH=1b​A​A⊤H=\\frac\{1\}\{b\}AA^\{\\top\}\. Then there exists universal polynomialspkp\_\{k\}such that we have almost surely

‖V⊤​Hk​V−V⊤​\(pk​\(C\)\)​V‖o​p→0\\\|V^\{\\top\}H^\{k\}V\-V^\{\\top\}\(p\_\{k\}\(C\)\)V\\\|\_\{op\}\\rightarrow 0
The first three such polynomials are

p1​\(C\)\\displaystyle p\_\{1\}\(C\)=C\\displaystyle=Cp2​\(C\)\\displaystyle p\_\{2\}\(C\)=C2\+δ​τ1​C\\displaystyle=C^\{2\}\+\\delta\\tau\_\{1\}Cp3​\(C\)\\displaystyle p\_\{3\}\(C\)=C3\+δ​τ1​C2\+\(δ​τ2\+δ2​τ12\)​C\\displaystyle=C^\{3\}\+\\delta\\tau\_\{1\}C^\{2\}\+\(\\delta\\tau\_\{2\}\+\\delta^\{2\}\\tau\_\{1\}^\{2\}\)C

###### Proof\.

Choose an orthonormal basiseie\_\{i\}ofℝr\\mathbb\{R\}^\{r\}, asVVis ann×rn\\times rmatrix\. It is sufficient to bound this operator norm by bounding the inner product over all vectors forming the orthonormal basis\. We thus fixx=V​eix=Ve\_\{i\}andy=V​ejy=Ve\_\{j\}and compute the inner product⟨x,H​y⟩\\langle x,Hy\\rangle\. By Cauchy’s Integral Formula, we have

⟨x,Hk​y⟩=12​π​i​∮Γzk​⟨x,\(W−z​I\)−1​y⟩​𝑑z\\langle x,H^\{k\}y\\rangle=\\frac\{1\}\{2\\pi i\}\\oint\_\{\\Gamma\}z^\{k\}\\langle x,\(W\-zI\)^\{\-1\}y\\rangle dzwhereΓ\\Gammais a contour which encloses the spectrum ofHH\. First we explicitly choose the contour for almost everyHH\. Note that

H=1b​C1/2​Z​Z⊤​C1/2H=\\frac\{1\}\{b\}C^\{1/2\}ZZ^\{\\top\}C^\{1/2\}whereZZhas iid standard Gaussian columns\. Therefore

‖H‖o​p≤‖C‖o​p​‖1b​Z​Z⊤‖o​p=K​‖1b​Z​Z⊤‖o​p\\\|H\\\|\_\{op\}\\leq\\\|C\\\|\_\{op\}\\\|\\frac\{1\}\{b\}ZZ^\{\\top\}\\\|\_\{op\}=K\\\|\\frac\{1\}\{b\}ZZ^\{\\top\}\\\|\_\{op\}By standard Marchenko\-Pasteur theory, the spectrum of1b​Z​Z⊤\\frac\{1\}\{b\}ZZ^\{\\top\}is almost surely eventually contained in the Marchenko\-Pasteur bulk\[\(1−δ\)2,\(1\+δ\)2\]\[\(1\-\\sqrt\{\\delta\}\)^\{2\},\(1\+\\sqrt\{\\delta\}\)^\{2\}\]\. Therefore almost surely

‖H‖o​p≤K​\(\(1\+δ\)2\+o​\(1\)\)\\\|H\\\|\_\{op\}\\leq K\(\(1\+\\sqrt\{\\delta\}\)^\{2\}\+o\(1\)\)ChooseR\>‖H‖o​pR\>\\\|H\\\|\_\{op\}\. Take the contourΓ\\Gammato be the circle of radiusRR\.

Now we use\[Couillet and Liao,[2022](https://arxiv.org/html/2605.11181#bib.bib69), Theorem 2\.6\]\. Then we have, similar to Equation 2\.46 of loc\. cit,

⟨x,Hk​y⟩=12​π​i​∮Γzk−1​⟨x,\(I\+m~p​\(z\)​C\)−1​y⟩​𝑑z\+o​\(1\)\\langle x,H^\{k\}y\\rangle=\\frac\{1\}\{2\\pi i\}\\oint\_\{\\Gamma\}z^\{k\-1\}\\langle x,\(I\+\\tilde\{m\}\_\{p\}\(z\)C\)^\{\-1\}y\\rangle dz\+o\(1\)wherem~p​\(z\)\\tilde\{m\}\_\{p\}\(z\)solves the equation

z=−1m~p​\(z\)\+δn​tr⁡\(C​\(I\+m~p​\(z\)​C\)−1\)z=\-\\frac\{1\}\{\\tilde\{m\}\_\{p\}\(z\)\}\+\\frac\{\\delta\}\{n\}\\operatorname\{tr\}\(C\(I\+\\tilde\{m\}\_\{p\}\(z\)C\)^\{\-1\}\)\(31\)By expandingm~p​\(z\)\\tilde\{m\}\_\{p\}\(z\)in a power series inz−1z^\{\-1\}we obtain a Laurent series\. Combining this with the expansion for\(I\+m~p​\(z\)​C\)−1\(I\+\\tilde\{m\}\_\{p\}\(z\)C\)^\{\-1\}gives a Laurent series inz−1z^\{\-1\}with coefficients which are polynomials inCC\. By standard complex analysis, the integral simply selects one of these coefficients\. This gives us the result\.

Now we explicitly compute the first three polynomials\. Expanding the second term of[Equation˜31](https://arxiv.org/html/2605.11181#A7.E31)we have

1n​tr⁡\(C​\(I\+m~p​\(z\)​C\)−1\)=τ1−m~p​\(z\)​τ2\+m~p​\(z\)2​τ3\+O​\(m~p​\(z\)3\)\\frac\{1\}\{n\}\\operatorname\{tr\}\(C\(I\+\\tilde\{m\}\_\{p\}\(z\)C\)^\{\-1\}\)=\\tau\_\{1\}\-\\tilde\{m\}\_\{p\}\(z\)\\tau\_\{2\}\+\\tilde\{m\}\_\{p\}\(z\)^\{2\}\\tau\_\{3\}\+O\(\\tilde\{m\}\_\{p\}\(z\)^\{3\}\)
This allows us to solve

m~p​\(z\)=−1z−δ​τ1z2−δ​τ2\+δ2​τ12z3\+O​\(z−4\)\\tilde\{m\}\_\{p\}\(z\)=\-\\frac\{1\}\{z\}\-\\frac\{\\delta\\tau\_\{1\}\}\{z^\{2\}\}\-\\frac\{\\delta\\tau\_\{2\}\+\\delta^\{2\}\\tau\_\{1\}^\{2\}\}\{z^\{3\}\}\+O\(z^\{\-4\}\)
For largezz,m~p​\(z\)\\tilde\{m\}\_\{p\}\(z\)is small hence we may expand

\(I\+m~p​\(z\)​C\)−1=I−m~p​\(z\)​C\+\(m~p​\(z\)​C\)2−\(m~p​\(z\)​C\)3\+O​\(C4\)\(I\+\\tilde\{m\}\_\{p\}\(z\)C\)^\{\-1\}=I\-\\tilde\{m\}\_\{p\}\(z\)C\+\(\\tilde\{m\}\_\{p\}\(z\)C\)^\{2\}\-\(\\tilde\{m\}\_\{p\}\(z\)C\)^\{3\}\+O\(C^\{4\}\)Then by plugging in the expansion form~p​\(z\)\\tilde\{m\}\_\{p\}\(z\)we obtain exactly the polynomials written above\. ∎

In order to actually evaluate the limit in this theorem we need to assume that fori=1,2,3i=1,2,3the following limit exists

V⊤​Ci​V→ΓiV^\{\\top\}C^\{i\}V\\rightarrow\\Gamma\_\{i\}
Then by the above theorem we have that fori=1,2,3i=1,2,3, and the polynomials as defined above

V⊤​Hi​V→pi​\(Γi,…,Γ1\)V^\{\\top\}H^\{i\}V\\rightarrow p\_\{i\}\(\\Gamma\_\{i\},\\dots,\\Gamma\_\{1\}\)We define fori=2,3i=2,3,Γ^i=pi​\(Γi,…,Γ1\)\\widehat\{\\Gamma\}\_\{i\}=p\_\{i\}\(\\Gamma\_\{i\},\\dots,\\Gamma\_\{1\}\)\.

Then almost surely

γSGD→tr⁡\(Σ2​Γ1\)tr⁡\(Σ2​Γ^2\)\\gamma\_\{\\texttt\{SGD\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma^\{2\}\\Gamma\_\{1\}\)\}\{\\operatorname\{tr\}\(\\Sigma^\{2\}\\widehat\{\\Gamma\}\_\{2\}\)\}\(32\)ΦSGD→tr\(Σ2Γ^\)2tr⁡\(Σ2​Γ^3\)\\Phi\_\{\\texttt\{SGD\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma^\{2\}\\widehat\{\\Gamma\}\)^\{2\}\}\{\\operatorname\{tr\}\(\\Sigma^\{2\}\\widehat\{\\Gamma\}\_\{3\}\)\}\(33\)
For Muon we have

D=polar⁡\(G​H\)=\(G​H2​G⊤\)−1/2​G​HD=\\operatorname\{polar\}\(GH\)=\(GH^\{2\}G^\{\\top\}\)^\{\-1/2\}GHDefineTi=V⊤​Hi​VT^\{i\}=V^\{\\top\}H^\{i\}V, so that almost surelyT1→Γ1T^\{1\}\\rightarrow\\Gamma\_\{1\}and fori=2,3i=2,3almost surelyTi→G^iT^\{i\}\\rightarrow\\widehat\{G\}\_\{i\}\. Then

γMuon\\displaystyle\\gamma\_\{\\texttt\{Muon\}\}=tr⁡\(\(G​H2​G⊤\)−1/2​G​H​G⊤\)tr⁡\(\(G​H2​G⊤\)1/2\)=tr⁡\(\(Σ​T2​Σ\)−1/2​Σ​T1​Σ\)tr⁡\(\(Σ​T2​Σ\)1/2\)\\displaystyle=\\frac\{\\operatorname\{tr\}\(\(GH^\{2\}G^\{\\top\}\)^\{\-1/2\}GHG^\{\\top\}\)\}\{\\operatorname\{tr\}\(\(GH^\{2\}G^\{\\top\}\)^\{1/2\}\)\}=\\frac\{\\operatorname\{tr\}\(\(\\Sigma T\_\{2\}\\Sigma\)^\{\-1/2\}\\Sigma T\_\{1\}\\Sigma\)\}\{\\operatorname\{tr\}\(\(\\Sigma T\_\{2\}\\Sigma\)^\{1/2\}\)\}\(34\)→tr⁡\(\(Σ​Γ^2​Σ\)−1/2​Σ​Γ1​Σ\)tr⁡\(\(Σ​Γ^2​Σ\)1/2\)\\displaystyle\\rightarrow\\frac\{\\operatorname\{tr\}\(\(\\Sigma\\widehat\{\\Gamma\}\_\{2\}\\Sigma\)^\{\-1/2\}\\Sigma\\Gamma\_\{1\}\\Sigma\)\}\{\\operatorname\{tr\}\(\(\\Sigma\\widehat\{\\Gamma\}\_\{2\}\\Sigma\)^\{1/2\}\)\}\(35\)ΦMuon\\displaystyle\\Phi\_\{\\texttt\{Muon\}\}=tr\(\(GH2G⊤\)1/2\)2tr⁡\(\(G​H2​G⊤\)−1​G​H3​G⊤\)=tr\(\(ΣT2Σ\)1/2\)2tr⁡\(\(Σ​T2​Σ\)−1​Σ​T3​Σ\)\\displaystyle=\\frac\{\\operatorname\{tr\}\(\(GH^\{2\}G^\{\\top\}\)^\{1/2\}\)^\{2\}\}\{\\operatorname\{tr\}\(\(GH^\{2\}G^\{\\top\}\)^\{\-1\}GH^\{3\}G^\{\\top\}\)\}=\\frac\{\\operatorname\{tr\}\(\(\\Sigma T\_\{2\}\\Sigma\)^\{1/2\}\)^\{2\}\}\{\\operatorname\{tr\}\(\(\\Sigma T\_\{2\}\\Sigma\)^\{\-1\}\\Sigma T\_\{3\}\\Sigma\)\}\(36\)→tr\(\(ΣΓ^2Σ\)1/2\)2tr\(\(ΣΓ^2Σ\)−1\(ΣΓ^3Σ\)\\displaystyle\\rightarrow\\frac\{\\operatorname\{tr\}\(\(\\Sigma\\widehat\{\\Gamma\}\_\{2\}\\Sigma\)^\{1/2\}\)^\{2\}\}\{\\operatorname\{tr\}\(\(\\Sigma\\widehat\{\\Gamma\}\_\{2\}\\Sigma\)^\{\-1\}\(\\Sigma\\widehat\{\\Gamma\}\_\{3\}\\Sigma\)\}\(37\)
### G\.1The Diagonal Case

If we assume thatCCis always aligned with the right singular vector basis ofRR, i\.e\. thatV⊤​C​V=ΛV^\{\\top\}CV=\\Lambdafor some diagonal matrixΛ\\Lambdathen we can say a lot more about these various formulas\.

Indeed we obtain,

Γi=limV⊤​Ci​V=Λi\\Gamma\_\{i\}=\\lim V^\{\\top\}C^\{i\}V=\\Lambda^\{i\}
We obtain thusly

γSGD→tr⁡\(Σ2​Λ\)tr⁡\(Σ2​\(Λ2\+δ​Λ\)\)\\gamma\_\{\\texttt\{SGD\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma^\{2\}\\Lambda\)\}\{\\operatorname\{tr\}\(\\Sigma^\{2\}\(\\Lambda^\{2\}\+\\delta\\Lambda\)\)\}\(38\)ΦSGD→tr\(Σ2\(Λ2\+δΛ\)\)2tr⁡\(Σ2​\(Λ3\+2​δ​Λ2\+\(δ\+δ2\)​Λ\)\)\\Phi\_\{\\texttt\{SGD\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma^\{2\}\(\\Lambda^\{2\}\+\\delta\\Lambda\)\)^\{2\}\}\{\\operatorname\{tr\}\(\\Sigma^\{2\}\(\\Lambda^\{3\}\+2\\delta\\Lambda^\{2\}\+\(\\delta\+\\delta^\{2\}\)\\Lambda\)\)\}\(39\)γMuon→tr⁡\(Σ​\(Λ2\+δ​Λ\)−1/2​Λ\)tr⁡\(Σ​\(Λ2\+δ​Λ\)1/2\)\\gamma\_\{\\texttt\{Muon\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma\(\\Lambda^\{2\}\+\\delta\\Lambda\)^\{\-1/2\}\\Lambda\)\}\{\\operatorname\{tr\}\(\\Sigma\(\\Lambda^\{2\}\+\\delta\\Lambda\)^\{1/2\}\)\}\(40\)ΦMuon→tr\(Σ\(Λ2\+δΛ\)1/2\)2tr⁡\(\(Λ2\+δ​Λ\)−1​\(Λ3\+2​δ​Λ2\+\(δ\+δ2\)​Λ\)\)\\Phi\_\{\\texttt\{Muon\}\}\\rightarrow\\frac\{\\operatorname\{tr\}\(\\Sigma\(\\Lambda^\{2\}\+\\delta\\Lambda\)^\{1/2\}\)^\{2\}\}\{\\operatorname\{tr\}\(\(\\Lambda^\{2\}\+\\delta\\Lambda\)^\{\-1\}\(\\Lambda^\{3\}\+2\\delta\\Lambda^\{2\}\+\(\\delta\+\\delta^\{2\}\)\\Lambda\)\)\}\(41\)
Definexi=σi​λi​\(λi\+δ\)x\_\{i\}=\\sigma\_\{i\}\\sqrt\{\\lambda\_\{i\}\(\\lambda\_\{i\}\+\\delta\)\}andhi=1λi\+δh\_\{i\}=\\frac\{1\}\{\\lambda\_\{i\}\+\\delta\}\. Then we rewrite the alignments as

γMuon=∑xi​hi∑xi\\gamma\_\{\\texttt\{Muon\}\}=\\frac\{\\sum x\_\{i\}h\_\{i\}\}\{\\sum x\_\{i\}\}\(42\)γSGD=∑xi2​hi∑xi2\\gamma\_\{\\texttt\{SGD\}\}=\\frac\{\\sum x\_\{i\}^\{2\}h\_\{i\}\}\{\\sum x\_\{i\}^\{2\}\}\(43\)
ForΦ\\Phidefine

Di=λi3\+2​δ​λi2\+\(δ\+δ2\)​λiλi2\+δ​λiD\_\{i\}=\\frac\{\\lambda\_\{i\}^\{3\}\+2\\delta\\lambda\_\{i\}^\{2\}\+\(\\delta\+\\delta^\{2\}\)\\lambda\_\{i\}\}\{\\lambda\_\{i\}^\{2\}\+\\delta\\lambda\_\{i\}\}Then we get

ΦMuon=\(∑xi\)2∑Di\\Phi\_\{\\texttt\{Muon\}\}=\\frac\{\(\\sum x\_\{i\}\)^\{2\}\}\{\\sum D\_\{i\}\}\(44\)ΦSGD=\(∑xi2\)2∑xi2​Di\\Phi\_\{\\texttt\{SGD\}\}=\\frac\{\(\\sum x\_\{i\}^\{2\}\)^\{2\}\}\{\\sum x\_\{i\}^\{2\}D\_\{i\}\}\(45\)
See[3\.1](https://arxiv.org/html/2605.11181#S3.Thmtheorem1)

###### Proof\.

Because we have that theλi\\lambda\_\{i\}are in non\-increasing order, and theσi\\sigma\_\{i\}are singular values hence are in non\-increasing order, we obtain that thexix\_\{i\}are also in non\-increasing order\. Finally, thehih\_\{i\}are in non\-decreasing order ashih\_\{i\}is decreasing inλi\\lambda\_\{i\}\.

Then

γMuon−γSGD=∑i,jxi​xj2​hi−∑i,jxi​xj2​hj=∑i<jxi​xj​\(xj−xi\)​\(hi−hj\)≥0\\gamma\_\{\\texttt\{Muon\}\}\-\\gamma\_\{\\texttt\{SGD\}\}=\\sum\_\{i,j\}x\_\{i\}x\_\{j\}^\{2\}h\_\{i\}\-\\sum\_\{i,j\}x\_\{i\}x\_\{j\}^\{2\}h\_\{j\}=\\sum\_\{i<j\}x\_\{i\}x\_\{j\}\(x\_\{j\}\-x\_\{i\}\)\(h\_\{i\}\-h\_\{j\}\)\\geq 0Equality holds if and only if one of the two sequencexix\_\{i\}orhih\_\{i\}are not strictly increasing or decreasing\. This is equivalent to one of the sequenceσi\\sigma\_\{i\}andλi\\lambda\_\{i\}being not strictly increasing\. ∎

See[3\.2](https://arxiv.org/html/2605.11181#S3.Thmtheorem2)

###### Proof\.

We simply compute that

Di=1\+3​δ\+δ21\+δD\_\{i\}=\\frac\{1\+3\\delta\+\\delta^\{2\}\}\{1\+\\delta\}Therefore

ΦSGD=\(1\+δ\)21\+3​δ\+δ2​∑σi2\\Phi\_\{\\texttt\{SGD\}\}=\\frac\{\(1\+\\delta\)^\{2\}\}\{1\+3\\delta\+\\delta^\{2\}\}\\sum\\sigma\_\{i\}^\{2\}ΦMuon=\(1\+δ\)21\+3​δ\+δ2​\(∑σi\)2d\\Phi\_\{\\texttt\{Muon\}\}=\\frac\{\(1\+\\delta\)^\{2\}\}\{1\+3\\delta\+\\delta^\{2\}\}\\frac\{\(\\sum\\sigma\_\{i\}\)^\{2\}\}\{d\}Then Cauchy\-Schwartz givesΦMuon≤ΦSGD\\Phi\_\{\\texttt\{Muon\}\}\\leq\\Phi\_\{\\texttt\{SGD\}\}\. ∎

See[3\.3](https://arxiv.org/html/2605.11181#S3.Thmtheorem3)

###### Proof\.

Forδ→∞\\delta\\rightarrow\\inftywe have

limδ→∞γMuon2​ΦMuon=1δ2​\(∑σi​λi1/2\)2d\\lim\_\{\\delta\\rightarrow\\infty\}\\gamma\_\{\\texttt\{Muon\}\}^\{2\}\\Phi\_\{\\texttt\{Muon\}\}=\\frac\{1\}\{\\delta^\{2\}\}\\frac\{\(\\sum\\sigma\_\{i\}\\lambda\_\{i\}^\{1/2\}\)^\{2\}\}\{d\}limδ→∞γSGD2​ΦSGD=∑σi2​λi\\lim\_\{\\delta\\rightarrow\\infty\}\\gamma\_\{\\texttt\{SGD\}\}^\{2\}\\Phi\_\{\\texttt\{SGD\}\}=\\sum\\sigma\_\{i\}^\{2\}\\lambda\_\{i\}Then by Cauchy\-Schwartz, we have

limδ→∞γMuon2​ΦMuon≤limδ→∞γSGD2​ΦSGD\\lim\_\{\\delta\\rightarrow\\infty\}\\gamma\_\{\\texttt\{Muon\}\}^\{2\}\\Phi\_\{\\texttt\{Muon\}\}\\leq\\lim\_\{\\delta\\rightarrow\\infty\}\\gamma\_\{\\texttt\{SGD\}\}^\{2\}\\Phi\_\{\\texttt\{SGD\}\}∎

## Appendix HRandom Feature Experiments

### H\.1Detailed Setup

For the Random Feature experiments considered in this paper, we utilize the setup and code fromDavis and Drusvyatskiy \[[2026](https://arxiv.org/html/2605.11181#bib.bib38)\], provided in[https://github\.com/damek/specgd/](https://github.com/damek/specgd/)\. The specific setup is as follows: We consider a random feature regression problem with objectivef​\(W\)=1N​d​‖W​A−Y‖F2f\(W\)=\\frac\{1\}\{N\\sqrt\{d\}\}\\\|WA\-Y\\\|\_\{F\}^\{2\}, whereW∈ℝo×dW\\in\\mathbb\{R\}^\{o\\times d\}is the weight matrix,A∈ℝd×NA\\in\\mathbb\{R\}^\{d\\times N\}is a fixed random feature matrix, andY=W⋆​AY=W^\{\\star\}Afor a randomly initialised ground\-truthW⋆W^\{\\star\}\. Since the loss is exactly quadratic inWW, the gradient isG=∇Wf=2N​d​\(W−W⋆\)​A​A⊤G=\\nabla\_\{W\}f=\\frac\{2\}\{N\\sqrt\{d\}\}\(W\-W^\{\\star\}\)AA^\{\\top\}and the Hessian action along any directionDDis⟨D,∇2f​\[D\]⟩=2N​d​‖D​A‖F2\\langle D,\\nabla^\{2\}f\[D\]\\rangle=\\frac\{2\}\{N\\sqrt\{d\}\}\\\|DA\\\|\_\{F\}^\{2\}\. Consequently, the exact optimal step size for any update directionDDis

η⋆​\(D\)=⟨G,D⟩1d​‖D​A‖F2\.\\eta^\{\\star\}\(D\)=\\frac\{\\langle G,D\\rangle\}\{\\frac\{1\}\{\\sqrt\{d\}\}\\\|DA\\\|\_\{F\}^\{2\}\}\.\(46\)
We run two experiments with this setup\. In theReLUexperiment,AAis obtained by applying a ReLU activation to a random linear map ofN=400N=400inputs, yieldingW∈ℝ120×100W\\in\\mathbb\{R\}^\{120\\times 100\},A∈ℝ100×400A\\in\\mathbb\{R\}^\{100\\times 400\}, stable rankst​\(A\)≈3\.06\\mathrm\{st\}\(A\)\\approx 3\.06,‖A‖op≈80\.87\\\|A\\\|\_\{\\mathrm\{op\}\}\\approx 80\.87, and‖A‖F≈141\.40\\\|A\\\|\_\{F\}\\approx 141\.40; training runs for10001000iterations\. In theSwiGLUexperiment the same dimensions are used but the activation is SwiGLU, givingst​\(A\)≈29\.96\\mathrm\{st\}\(A\)\\approx 29\.96,‖A‖op≈22\.15\\\|A\\\|\_\{\\mathrm\{op\}\}\\approx 22\.15, and‖A‖F≈121\.22\\\|A\\\|\_\{F\}\\approx 121\.22; training runs for400400iterations\. Both experiments use base learning rateη=10−2\\eta=10^\{\-2\}, double precision \(float64\), and fixed random seed\. The Optimal C methods selectccandη\\etajointly at each step by maximising the quadratic lower boundΔ​f​\(c,η\)=−η​a​\(c\)\+η2​b​\(c\)\\Delta f\(c,\\eta\)=\-\\eta\\,a\(c\)\+\\eta^\{2\}b\(c\)over a grid, wherea​\(c\)=⟨G,Dc⟩a\(c\)=\\langle G,D\_\{c\}\\rangleandb​\(c\)=1d​‖Dc​A‖F2b\(c\)=\\frac\{1\}\{\\sqrt\{d\}\}\\\|D\_\{c\}A\\\|\_\{F\}^\{2\}withDc=\(G​G⊤\)−c​GD\_\{c\}=\(GG^\{\\top\}\)^\{\-c\}G\.Freonis implemented using SVD, whileKaonis implemented with SVD with random uniform noise\.

### H\.2Exponent Approximation for Optimal Schatten Norm

We are in the random feature model setting and assume to haveG∈ℝm×nG\\in\\mathbb\{R\}^\{m\\times n\}be gradient with a compact Singular Value DecompositionG=U​S​V⊤G=USV^\{\\top\}, whereS=diag​\(σ1,…,σr\)S=\\text\{diag\}\(\\sigma\_\{1\},\\dots,\\sigma\_\{r\}\)\. We seek a spectrally scaled matrixD=U​Sa​V⊤D=US^\{a\}V^\{\\top\}that maximizes optimal decreaseΦk\\Phi\_\{k\}\.

The objective is to maximize decrease at optimal step size, for‖D​A‖F2=∑ki​σi2​a\\\|DA\\\|\_\{F\}^\{2\}=\\sum k\_\{i\}\\sigma\_\{i\}^\{2a\}, whereki=‖A⊤​vi‖2k\_\{i\}=\\\|A^\{\\top\}v\_\{i\}\\\|^\{2\}represents the projection energy of an activation matrixAA:

maxa⁡g​\(a\)=\(∑i=1rσia\+1\)2∑i=1rki​σi2​a\\max\_\{a\}g\(a\)=\\frac\{\\left\(\\sum\_\{i=1\}^\{r\}\\sigma\_\{i\}^\{a\+1\}\\right\)^\{2\}\}\{\\sum\_\{i=1\}^\{r\}k\_\{i\}\\sigma\_\{i\}^\{2a\}\}\(47\)Setting the derivative oflog⁡g​\(a\)\\log g\(a\)to zero yields the stationarity condition for optimality:

∑i=1rσia\+1​log⁡σi∑i=1rσia\+1=∑i=1rki​σi2​a​log⁡σi∑i=1rki​σi2​a\\frac\{\\sum\_\{i=1\}^\{r\}\\sigma\_\{i\}^\{a\+1\}\\log\\sigma\_\{i\}\}\{\\sum\_\{i=1\}^\{r\}\\sigma\_\{i\}^\{a\+1\}\}=\\frac\{\\sum\_\{i=1\}^\{r\}k\_\{i\}\\sigma\_\{i\}^\{2a\}\\log\\sigma\_\{i\}\}\{\\sum\_\{i=1\}^\{r\}k\_\{i\}\\sigma\_\{i\}^\{2a\}\}\(48\)
To approximate an optimalaavalue, we will assume the singular values ofGGand the penalty weights follow a power\-law decay:

σi∝i−α\(α\>0\),ki∝i−β\\sigma\_\{i\}\\propto i^\{\-\\alpha\}\\quad\(\\alpha\>0\),\\qquad k\_\{i\}\\propto i^\{\-\\beta\}\(49\)In this case, equality is achieved by equating exponents of LHS and RHS, i\.e\.α​\(a\+1\)=2​a​α\+β\\alpha\(a\+1\)=2a\\alpha\+\\beta, meaning that

a=1−βαa=1\-\\frac\{\\beta\}\{\\alpha\}

### H\.3Optimal C methods

Both adaptive\-exponent methods take the update directionDc=\(G​G⊤\)−c​G=U​diag​\(σ~1−2​c\)​V⊤D\_\{c\}=\(GG^\{\\top\}\)^\{\-c\}G=U\\,\\mathrm\{diag\}\(\\tilde\{\\sigma\}^\{1\-2c\}\)\\,V^\{\\top\}, whereG=U​Σ​V⊤G=U\\Sigma V^\{\\top\}is the SVD of the gradient andσ~i=σi/μc\\tilde\{\\sigma\}\_\{i\}=\\sigma\_\{i\}/\\mu\_\{c\}are the singular values normalised by the power meanμc=\(1r​∑iσi2​\(1−c\)\)1/\(2​\(1−c\)\)\\mu\_\{c\}=\\bigl\(\\tfrac\{1\}\{r\}\\sum\_\{i\}\\sigma\_\{i\}^\{2\(1\-c\)\}\\bigr\)^\{1/\(2\(1\-c\)\)\}\. The optimal step sizeη⋆​\(c\)=n​a​\(c\)/b​\(c\)\\eta^\{\\star\}\(c\)=n\\,a\(c\)/b\(c\)is used exactly at every step, wherea​\(c\)=∑iσi​σ~i1−2​ca\(c\)=\\sum\_\{i\}\\sigma\_\{i\}\\,\\tilde\{\\sigma\}\_\{i\}^\{1\-2c\}andb​\(c\)=‖Dc​A‖F2/db\(c\)=\\\|D\_\{c\}A\\\|\_\{F\}^\{2\}/\\sqrt\{d\}\.

##### Optimal C \[Greedy\]

determinesccby exhaustive search: at each iteration a grid of 41 valuesc∈\[−0\.5,1\.5\]c\\in\[\-0\.5,1\.5\]is evaluated and the value minimising the quadratic predicted loss decreaseΔ​f⋆​\(c\)=−n​a​\(c\)2/\(2​b​\(c\)\)\\Delta f^\{\\star\}\(c\)=\-n\\,a\(c\)^\{2\}/\(2\\,b\(c\)\)is selected\.

##### Optimal C \[Scaling\]

instead infersccanalytically from the spectral structure of the gradient and the feature matrix based on[Section˜H\.2](https://arxiv.org/html/2605.11181#A8.SS2)\. Fittinglog⁡σi≈−α​log⁡i\\log\\sigma\_\{i\}\\approx\-\\alpha\\log iandlog⁡ki≈−β​log⁡i\\log k\_\{i\}\\approx\-\\beta\\log iby least squares, whereki=‖vi⊤​A‖2k\_\{i\}=\\\|v\_\{i\}^\{\\top\}A\\\|^\{2\}is the squared projection of theii\-th right singular vector onto the feature matrix, the theoretically optimal exponent is

c⋆=β2​α,c^\{\\star\}=\\frac\{\\beta\}\{2\\alpha\},\(50\)which is clamped to\[−0\.5,1\.5\]\[\-0\.5,1\.5\]and used with step sizeη⋆​\(c⋆\)\\eta^\{\\star\}\(c^\{\\star\}\)\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/x9.png)Figure H\.1:Random\-feature regression with SwiGLU activations\.Same setup as Figure[5](https://arxiv.org/html/2605.11181#S3.F5)but with a SwiGLU random\-feature matrixA∈ℝ100×400A\\in\\mathbb\{R\}^\{100\\times 400\}\(st​\(A\)≈29\.96\\mathrm\{st\}\(A\)\\approx 29\.96, substantially higher than the ReLU case\), trained for 400 iterations\.*Left:*training loss\.*Centre:*effective step sizes; Optimal GD exhibits strongly oscillatory step sizes throughout training\.*Right:*exponentccin\(G​G⊤\)−c​G\(GG^\{\\top\}\)^\{\-c\}Gfor the two Optimal C methods\. The higher stable rank ofAAunder SwiGLU changes the curvature landscape relative to the ReLU setting, and both Optimal C methods converge to a smaller value ofcc\(closer to GD\) than in the ReLU experiment\.

## Appendix IGPT2 Additional Experiments

![Refer to caption](https://arxiv.org/html/2605.11181v1/x10.png)Figure I\.1:Per\-layer curvature alignment quantitiesΦk\\Phi\_\{k\}\(top row, linear scale\) andγk\\gamma\_\{k\}\(bottom row, linear scale\) for the weight matrixh\.3\.attn\.c\_attnat four training checkpoints \(steps 1, 1797, 3594, 5391\)\. Each bar shows the mean±\\pmone standard deviation over validation batches, separately for each candidate update directionDkD\_\{k\}\. The directions includeFreon’s power\-mean normalised descent direction at interpolation levelsc∈0,0\.25,0\.5,0\.75,1c\\in\{0,0\.25,0\.5,0\.75,1\}\(wherec=0c=0matchesSGD, andc=1c=1is the fullFreondirection\),TruncatedSGD\(batch gradient with the top 1% of singular values zeroed\), andKaon\. Formally, for a given directionDkD\_\{k\}and batch gradientgbg\_\{b\}, we calculate⟨gb,Dk⟩2\\langle g\_\{b\},D\_\{k\}\\rangle^\{2\}, the directional curvatureλk=⟨Dk,ℱval​\[Dk\]⟩\\lambda\_\{k\}=\\langle D\_\{k\},\\mathcal\{F\}\_\{\\mathrm\{val\}\}\[D\_\{k\}\]\\ranglewhereℱval\\mathcal\{F\}\_\{\\mathrm\{val\}\}is the Gauss–Newton \(GGN\) matrix estimated over all validation batches via Jacobian–vector products, and the batch gradient alignmentγk=⟨Gk,Dk⟩/⟨gb,Dk⟩\\gamma\_\{k\}=\\langle G\_\{k\},D\_\{k\}\\rangle/\\langle g\_\{b\},D\_\{k\}\\ranglewhereGkG\_\{k\}is the gradient averaged over all validation batches\. A largeΦk\\Phi\_\{k\}indicates that the direction is strongly aligned with the batch gradient and lies in a low\-curvature region, supporting a larger step, whileγk≈1\\gamma\_\{k\}\\approx 1indicates that the batch gradient faithfully represents the global gradient in that direction\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x11.png)Figure I\.2:Validation loss as a function of training progress for each optimiser at its best learning rate on WikiText\-2\. Curves shown for Freon with exponentc∈\{0,14,12,34,1\}c\\in\\\{0,\\tfrac\{1\}\{4\},\\tfrac\{1\}\{2\},\\tfrac\{3\}\{4\},1\\\}, 5% TruncatedSGD, and Kaon\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x12.png)Figure I\.3:Loss changeΔ​f\\Delta fover the joint\(c,η\)\(c,\\eta\)grid for the quadratic model, evaluated at three training iterations of the random feature experiment\. Each panel uses the gradient computed at that iteration to evaluateΔ​f​\(c,η\)=−η​a​\(c\)\+η2​b​\(c\)\\Delta f\(c,\\eta\)=\-\\eta\\,a\(c\)\+\\eta^\{2\}b\(c\), wherea​\(c\)=⟨G,Dc⟩a\(c\)=\\langle G,D\_\{c\}\\rangleandb​\(c\)=12​n​‖Dc​A‖F2b\(c\)=\\frac\{1\}\{2n\}\\\|D\_\{c\}A\\\|\_\{F\}^\{2\}\. The dashed black curve marks the optimal learning rate percc; the green dotted line marks the actual training learning rateη=0\.01\\eta=0\.01\.Δp=‖G‖dual1⁡\(p\)2‖A‖dual2⁡\(p\)2,\\Delta\_\{p\}=\\tfrac\{\\\|G\\\|^\{2\}\_\{\\operatorname\{dual\}\_\{1\}\(p\)\}\}\{\\\|A\\\|^\{2\}\_\{\{\\operatorname\{dual\}\_\{2\}\(p\)\}\}\},\(51\)for1p\+1dual1⁡p=1\\tfrac\{1\}\{p\}\+\\tfrac\{1\}\{\\operatorname\{dual\}\_\{1\}p\}=1and1p\+1dual2⁡p=12\\tfrac\{1\}\{p\}\+\\tfrac\{1\}\{\\operatorname\{dual\}\_\{2\}p\}=\\frac\{1\}\{2\}\.

![Refer to caption](https://arxiv.org/html/2605.11181v1/x13.png)Figure I\.4:Layer\-wise optimalccat three training snapshots \(RandFeat quadratic model\)\.Each cell shows the value ofc∈\[0,1\]c\\in\[0,1\]that minimises the predicted one\-step loss decrease for that layer under a random\-feature approximation to the curvature\. Panels correspond to the beginning, middle, and end of training\. Rows correspond to the four weight matrices per transformer block \(WQ​K​VW\_\{QKV\},WOW\_\{O\},WMLP,1W\_\{\\mathrm\{MLP\},1\},WMLP,2W\_\{\\mathrm\{MLP\},2\}\); columns index the 12 transformer blocks\.c=0c=0recovers gradient descent,c=12c=\\tfrac\{1\}\{2\}recovers Muon, andc=1c=1recovers the pseudoinverse\-like endpoint\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x14.png)Figure I\.5:Per\-step loss decrease vs\. learning rate andccunder three local quadratic models\(layerh\.2\.attn\.c\_attn, step 800\)\.*Left \(Frobenius\):*The naive random feature model predictsc≈2/3c\\approx 2/3as optimal\.*Centre \(JVP\):*Under the GGN quadratic approximation on a single batch,c=1c=1achieves the largest per\-step decrease – similar to the Newton step being optimal for a fixed\-batch quadratic\.*Right \(Actual\):*Evaluating on the actual single\-batch loss agrees with the GGN picture, motivating the inter\-batch SNR analysis of[Section˜3\.2](https://arxiv.org/html/2605.11181#S3.SS2)\. All heatmaps are clamped at zero \(only improvements shown\); the dashed black curve marks the per\-ccoptimal learning rate and the green dotted line the configured learning rate\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x15.png)Figure I\.6:Validation loss changeΔ​f​\(c,η\)=f​\(θ−η,Dc\)−f​\(θ\)\\Delta f\(c,\\eta\)=f\(\\theta\-\\eta,D\_\{c\}\)\-f\(\\theta\), evaluated on the full validation set after applying a single synchronised update to all model layers, over a joint grid of direction interpolationc∈\[0,1\]c\\in\[0,1\]\(y\-axis\) and learning rateη\\eta\(x\-axis, logarithmic\)\. Each panel corresponds to a different training checkpoint, i\.e\. the gradientGkG\_\{k\}and momentum state used to construct the update direction are taken from that training step\. For a givencc, the per\-layer update directionDcD\_\{c\}is obtained by applying the Freon power\-mean preconditioner with exponentccto the layer’s gradient:c=0c=0recovers the vanilla gradient\-descent direction \(D=G/\|G\|FD=G/\|G\|\_\{F\}\) andc=1c=1is the fully preconditionedFreondirection; intermediate values interpolate smoothly between the two\. The sameccandη\\etaare applied to every weight matrix simultaneously, soΔ​f\\Delta fmeasures the*joint*effect of a global update rather than a per\-layer optimum\. The heatmap is clamped at zero so that only loss\-decreasing regions \(negativeΔ​f\\Delta f\) are shown; the colour axis is scaled independently per panel\. The dashed black curve tracesη∗​\(c\)=arg⁡minη⁡Δ​f​\(c,η\)\\eta^\{\*\}\(c\)=\\arg\\min\_\{\\eta\}\\Delta f\(c,\\eta\), i\.e\. the empirically optimal learning rate for each direction interpolationcc\. The green dotted vertical line marks the actual learning rate used during training, allowing direct comparison of where training operates relative to the optimal ridge\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x16.png)Figure I\.7:Layer\-wise range of near\-optimalccat step 800\.Each cell shows the range of valuesc∈\[0,1\]c\\in\[0,1\]whose best\-over\-η\\etaone\-step loss decrease is within 10% of the global optimum for that layer, found by a grid search over learning rates andccusing real forward passes on the full validation set \(actual mode\)\. Cell colour indicates the optimalcc\. Rows correspond to the four weight matrices per transformer block \(WQ​K​VW\_\{QKV\},WOW\_\{O\},WMLP,1W\_\{\\mathrm\{MLP\},1\},WMLP,2W\_\{\\mathrm\{MLP\},2\}\); columns index the 12 transformer blocks\.c=0c=0recovers gradient descent,c=12c=\\tfrac\{1\}\{2\}recovers Muon, andc=1c=1recovers the pseudoinverse\-like endpoint\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x17.png)Figure I\.8:Φk\\Phi\_\{k\}at training step 3594 \(≈25%\\approx\\\!25\\%of training\) for all 48 weight matrices, arranged as 12 transformer blocks \(rows\)×\\times4 sublayer types \(columns:WQ​K​VW\_\{QKV\},WOW\_\{O\},WFC1W\_\{\\mathrm\{FC1\}\},WFC2W\_\{\\mathrm\{FC2\}\}\)\. Each bar shows the mean±\\pmstd over validation batches for a given update direction\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x18.png)Figure I\.9:Batch\-gradient alignmentγk=⟨G,D⟩/⟨gb,D⟩\\gamma\_\{k\}=\\langle G,D\\rangle/\\langle g\_\{b\},D\\rangleat training step 3594 for all 48 weight matrices \(same layout as Fig\.[I\.8](https://arxiv.org/html/2605.11181#A9.F8)\)\. A value ofγk=1\\gamma\_\{k\}=1indicates the directionDDaligns perfectly with the full gradient; values below11reflect mini\-batch noise or directional mismatch\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x19.png)Figure I\.10:Loss changeΔ​f\\Delta fover the joint\(c,η\)\(c,\\eta\)grid at training step 800, evaluated per layer on thefull validation set\. Rows correspond to transformer blocksh​\.0h\.0–h​\.11h\.11; columns correspond to the four weight matrices within each block\. The dashed black curve marks the optimal learning rate percc; the green dotted line marks the actual training learning rate\.![Refer to caption](https://arxiv.org/html/2605.11181v1/x20.png)Figure I\.11:Loss changeΔ​f\\Delta fover the joint\(c,η\)\(c,\\eta\)grid at training step 800, evaluated per layer on thesame training batch as gradient\. Rows correspond to transformer blocksh​\.0h\.0–h​\.11h\.11; columns correspond to the four weight matrices within each block\. The dashed black curve marks the optimal learning rate percc; the green dotted line marks the actual training learning rate\.