Reassessing Muon for Matrix Factorization

# Reassessing Muon for Matrix Factorization
Source: [https://arxiv.org/html/2607.13246](https://arxiv.org/html/2607.13246)
\\theorembodyfont\\theoremheaderfont\\theorempostheader

:\\theoremsep \\jmlrvolume334\\jmlryear2026\\jmlrworkshopTopology, Algebra, and Geometry in Data Science

\\NameAli Parviz222Halicioğlu Data Science Institute, UC San Diego\.\\Emailalparviz@ucsd\.edu \\NameGal Mishne222Halicioğlu Data Science Institute, UC San Diego\.\\Emailgmishne@ucsd\.edu \\NameAlex Cloninger111Department of Mathematics, UC San Diego\.222Halicioğlu Data Science Institute, UC San Diego\.\\Emailacloninger@ucsd\.edu

###### Abstract

Muonhas recently emerged as a strong optimizer for large\-scale deep learning, where it reshapes gradient updates through approximate orthogonalization and has been reported to outperformAdamandAdamWon large language model training\. Its empirical success has motivated a growing theoretical literature that interpretsMuonas steepest descent under the spectral norm\. Yet it remains unclear which ofMuon’s advantages stem from its update rule itself and which are artifacts of the scale, architecture, and data of modern deep networks\. In this work we isolate the optimizer from these confounders by studyingMuonon a simple, well\-understood, and spectrally structured problem: low\-rank matrix factorization\. Through a controlled and systematically tuned comparison against adaptive baselines, we find thatMuondoes*not*consistently outperformAdamWin this setting, and that several previously reported advantages are sensitive to hyperparameter choices\. Our results give a more nuanced picture of when spectrum\-aware orthogonalization helps, and argue for evaluating modern optimizers on controlled problems in addition to end\-to\-end benchmarks\.

## 1Introduction

Recent advances in large\-scale optimization have introduced a class of optimizers that explicitly exploit the matrix structure of gradient updates\. Among these,Muon\(MomentUm Orthogonalized by Newton–Schulz\) has emerged as a promising alternative to standard first\-order methods\(jordan2024muon\)\. Rather than applying a momentum update directly,Muontransforms the update direction through an approximate orthogonalization step based on Newton–Schulz iterations, effectively reshaping the spectrum of the gradient\.

Empirically,Muonhas demonstrated strong performance on modern deep learning workloads, including GPT\-style models, where it improves training efficiency and in some cases outperforms widely used optimizers such asAdamandAdamW\(jordan2024muon;liu2025muonscalablellmtraining;shah2025practicalefficiencymuon\)\. Notably, large\-scale studies report up to a twofold speedup overAdamWin multi\-billion\-parameter language model training\(shah2025practicalefficiencymuon\)\. These successes have motivated a growing body of theoretical work that studiesMuonthrough the lenses of spectral\-norm steepest descent and constrained optimization\. Existing analyses establish convergence guarantees under various simplifying assumptions, including exact orthogonalization, simplified momentum dynamics, and locally quadratic models\(pethick2025normconstrainedlmo;li2025noteconvergencemuon;shen2025convergencemuon;chen2025muonspectralnorm;kovalev2025understandinggo;riabinin2025gluon;gruntkowska2025ef21muon;sato2025muonconvergence;nagashima2026improvedconvergenceratesmuon\), while more recent work studies inexact Newton–Schulz iterations and how approximation errors propagate into convergence guarantees\(shulgin2025beyondideal;anon2025newtonschulz;lau2025polargrad\)\.

Despite this progress, it remains unclear which aspects ofMuon’s advantage are intrinsic to the optimizer and which arise from the complexity of large\-scale deep learning, where architectural, scale, and data\-dependent effects are hard to disentangle\. Existing analyses rarely show*when*and*why*Muonshould outperform classical optimizers in concrete, well\-specified problems, tending to concern idealized settings or purely local properties\(davis2025whenspectral;su2025isotropic\)\. We take a complementary perspective and studyMuonon low\-rank matrix factorization: a simple yet fundamental problem that allows systematic, exhaustive evaluation and whereMuon’s spectral nature might be expected to help, since optimization is governed by well\-characterized curvature and singular\-value structure\. Our goal is to characterize both the strengths and limitations ofMuonand to identify regimes where it does or does not outperform standard optimizers such asAdamandAdamW\.

This focus exposes a gap between expectation and behavior, which we organize around three questions\.

1. Q1:DoesMuonoffer any advantage in structured problems such as matrix factorization?Matrix factorization has explicit low\-rank structure and is associated with a well\-behaved optimization landscape\. IfMuongenuinely exploits spectral properties of the gradient, this should translate into faster convergence\. If it does not, what limits its effectiveness relative to adaptive methods such asAdamW?
2. Q2:AreMuon’s benefits tied to specific problem structures?Which properties of an optimization problem make orthogonalized updates advantageous? Do their benefits stem from specific geometric or spectral features of the objective, or do they generalize across matrix factorization tasks with varying conditioning, over\-parameterization, and constraints?
3. Q3:How doesMuoncompare to adaptive methods under ill\-conditioning?Matrix factorization is sensitive to conditioning, especially when singular values decay rapidly\. Optimizers such asAdamWimplicitly adapt to coordinate\-wise scaling\. DoesMuon’s spectral normalization compete with or conflict with such adaptivity, and under which regimes does one dominate?

Findings\.Contrary to its strong performance in large\-scale models, Muon’s advantage in these canonical settings is problem\-dependent\. On low\-rank factorization and matrix completion, its reported gains disappear under equal tuning, with AdamW and GD matching or exceeding its performance\. On NMF, however, Muon retains a consistent advantage, likely because its orthogonalized updates discourage redundant factors and promote more diverse representations\. This dependence on problem structure suggests that Muon’s broader success may also arise from properties specific to deep learning and underscores the value of controlled optimization benchmarks\.

## 2Matrix factorization: symmetric and nonnegative variants

We begin by fixing notation used throughout the paper\. For a matrix𝑴\{\\bm\{M\}\}, letσi​\(𝑴\)\\sigma\_\{i\}\(\{\\bm\{M\}\}\)denote itsii\-th largest singular value andσmin​\(𝑴\)\\sigma\_\{\\min\}\(\{\\bm\{M\}\}\)its smallest\. We write‖𝑴‖op\\left\\lVert\{\\bm\{M\}\}\\right\\rVert\_\{\\mathrm\{op\}\}and‖𝑴‖F\\left\\lVert\{\\bm\{M\}\}\\right\\rVert\_\{\\mathrm\{F\}\}for the spectral and Frobenius norms, respectively\. Fork≤dk\\leq d, let𝒪d×k\\mathcal\{O\}\_\{d\\times k\}be the set of matrices inℝd×k\\mathbb\{R\}^\{d\\times k\}with orthonormal columns\. Given scalarsa1,…,ada\_\{1\},\\dots,a\_\{d\}, we write𝖽𝗂𝖺𝗀​\{a1,…,ad\}\\mathsf\{diag\}\\\{a\_\{1\},\\dots,a\_\{d\}\\\}for the diagonal matrix with these entries\. We discussMuon\(jordan2024muon\)in[Appendix˜A](https://arxiv.org/html/2607.13246#A1)and defer a detailed discussion of related work to Appendix[C](https://arxiv.org/html/2607.13246#A3)\.

##### Symmetric matrix factorization

We consider the symmetric matrix factorization problem

min𝑼∈ℝd×kf​\(𝑼\)=14​‖𝑼​𝑼⊤−𝑴⋆‖F2,\\min\_\{\{\\bm\{U\}\}\\in\\mathbb\{R\}^\{d\\times k\}\}\\quad f\(\{\\bm\{U\}\}\)\\;=\\;\\frac\{1\}\{4\}\\,\\left\\lVert\{\\bm\{U\}\}\{\\bm\{U\}\}^\{\\top\}\-\{\\bm\{M\}\}^\{\\star\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\},\(1\)where𝑴⋆∈ℝd×d\{\\bm\{M\}\}^\{\\star\}\\in\\mathbb\{R\}^\{d\\times d\}is a rank\-rrpositive semidefinite matrix and𝑼∈ℝd×k\{\\bm\{U\}\}\\in\\mathbb\{R\}^\{d\\times k\}is a factor withk≥rk\\geq rcolumns \(k=rk=ris the exactly parameterized case andk\>rk\>ris over\-parameterized\)\. The goal is to recover𝑴⋆\{\\bm\{M\}\}^\{\\star\}through a low\-rank factorization𝑼​𝑼⊤\{\\bm\{U\}\}\{\\bm\{U\}\}^\{\\top\}by solving equation[1](https://arxiv.org/html/2607.13246#S2.E1)\. We assume𝑴⋆\{\\bm\{M\}\}^\{\\star\}admits the eigendecomposition𝑴⋆=𝑽⋆​𝚲⋆​𝑽⋆⊤\{\\bm\{M\}\}^\{\\star\}=\{\\bm\{V\}\}^\{\\star\}\{\\bm\{\\Lambda\}\}^\{\\star\}\{\\bm\{V\}\}^\{\\star\\top\}, where𝚲⋆=𝖽𝗂𝖺𝗀​\{λ1⋆,…,λr⋆\}\{\\bm\{\\Lambda\}\}^\{\\star\}=\\mathsf\{diag\}\\\{\\lambda\_\{1\}^\{\\star\},\\dots,\\lambda\_\{r\}^\{\\star\}\\\}collects the nonzero eigenvalues withλ1⋆≥⋯≥λr⋆\>0\\lambda\_\{1\}^\{\\star\}\\geq\\cdots\\geq\\lambda\_\{r\}^\{\\star\}\>0and𝑽⋆∈ℝd×r\{\\bm\{V\}\}^\{\\star\}\\in\\mathbb\{R\}^\{d\\times r\}is orthonormal\. The condition number of𝑴⋆\{\\bm\{M\}\}^\{\\star\}is

κ≔λ1⋆/λr⋆\.\\kappa\\;\\coloneqq\\;\\lambda\_\{1\}^\{\\star\}/\\lambda\_\{r\}^\{\\star\}\.\(2\)The symmetric problem equation[1](https://arxiv.org/html/2607.13246#S2.E1)is the positive semidefinite, factor\-tied instance of a broader family of low\-rank factorization problems; we defer the general bilinear factorization and the matrix completion variant to Appendix[B](https://arxiv.org/html/2607.13246#A2)\.

##### Nonnegative matrix factorization\.

A closely related variant constrains both the target and the factors to be entrywise nonnegative\. For a rank\-rrtarget𝑴⋆∈ℝd1×d2\{\\bm\{M\}\}^\{\\star\}\\in\\mathbb\{R\}^\{d\_\{1\}\\times d\_\{2\}\}, nonnegative matrix factorization \(NMF\) solves

min𝑼∈ℝ≥0d1×k𝑽∈ℝ≥0d2×k12​‖𝑼​𝑽⊤−𝑴⋆‖F2,\\min\_\{\\begin\{subarray\}\{c\}\{\\bm\{U\}\}\\in\\mathbb\{R\}\_\{\\geq 0\}^\{d\_\{1\}\\times k\}\\\\\[2\.0pt\] \{\\bm\{V\}\}\\in\\mathbb\{R\}\_\{\\geq 0\}^\{d\_\{2\}\\times k\}\\end\{subarray\}\}\\quad\\frac\{1\}\{2\}\\,\\left\\lVert\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}\-\{\\bm\{M\}\}^\{\\star\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\},\(3\)whereℝ≥0\\mathbb\{R\}\_\{\\geq 0\}denotes the nonnegative reals andk≥rk\\geq r\. NMF is the nonnegativity\-constrained special case of the general bilinear factorization equation[12](https://arxiv.org/html/2607.13246#A2.E12)\(Appendix[B](https://arxiv.org/html/2607.13246#A2)\)\. Moreover, We analyze the exact softplus and its Taylor surrogate as nonlinear matrix factorizations and show the truncation imposes a rank cap that explains the high\-rank gap in Fig\.[12](https://arxiv.org/html/2607.13246#A11.F12)\(Appendix\.[K](https://arxiv.org/html/2607.13246#A11)\)\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x1.png)Figure 1:Low\-rank factorization, conditioning sweep\(target rankr=15r\{=\}15, matched search rankk=15k\{=\}15, dimensiond=100d\{=\}100\)\. Tuned final loss vs\. learning rate for each condition numberκ\\kappa, from well\-conditioned \(κ=1\\kappa\{=\}1\) to strongly ill\-conditioned \(κ=625\\kappa\{=\}625\) \(geometric mean over33seeds,±\\pmlog\-std band\)\. The target’s1515singular values are spaced linearly fromκ\\kappadown to11, soκ\\kappais exactly the condition number\. When tuned, AdamW and GD reach near machine precision; the Muon configuration plateaus several orders higher\.![Refer to caption](https://arxiv.org/html/2607.13246v1/x2.png)Figure 2:Non\-negative factorization\.Final loss \(yy\-axis, log scale; geometric mean over33seeds, shaded±\\pmone\-log\-standard\-deviation band\) versus learning rate \(xx\-axis, log scale\)\. Panels\(a–c\)use the uniform spectrum and\(d–f\)the decayed spectrum, at factor ranksr=10,50,100r=10,50,100respectively\. In both spectra, tuned Muon is the only method to fit the target across ranks; AdamW trails by several orders and GD/SignGD fail\.![Refer to caption](https://arxiv.org/html/2607.13246v1/images/muon-spectrum.png)Figure 3:AdamW outperforms Muon variants across spectrum shapes\(matrix factorization, targetm×n=150×150m\{\\times\}n=150\{\\times\}150, rankr=5r\{=\}5\)\. Tuned final loss for each of the seven spectral profiles at fixedκ=104\\kappa=10^\{4\}\(mean over33seeds\)\. Both Muon and Muon without Nesterov\-style momentum exhibit degraded relative performance on most spectral shapes, whereas AdamW remains effective despite the extreme ill\-conditioning\.![Refer to caption](https://arxiv.org/html/2607.13246v1/x3.png)Figure 4:Spectral subspace recovery\(Gaussian\-kernel targetK∈ℝ100×100K\\in\\mathbb\{R\}^\{100\\times 100\}, full rankk=100k\{=\}100, tuned LR,1212inits\)\. AlignmentMi​j=\|⟨ui,ej⟩\|M\_\{ij\}=\\lvert\\langle u\_\{i\},e\_\{j\}\\rangle\\rvert\(Eq\.[16](https://arxiv.org/html/2607.13246#A10.E16)\) between learned singular vectors and true eigenvectors ofKK\.\(a\)AdamW and\(b\)Muon both recover the top eigenvectors, diverging in the near\-degenerate mid\-band\.\(c\)The difference isolates it: AdamW holds a tighter ordered diagonal \(blue\), Muon disperses onto neighbours \(red\) — lower diagonal energy \(0\.5860\.586vs\.0\.6430\.643\)\.![Refer to caption](https://arxiv.org/html/2607.13246v1/x4.png)Figure 5:Tensor\-Train Noiseless regime\(true rankr⋆=4r^\{\\star\}=4\)\. Final loss versus learning rate for each optimizer, across search ranksr∈\{2,4,8,20\}r\\in\\\{2,4,8,20\\\}spanning under\-, exactly\- \(r=r⋆r=r^\{\\star\}, bold panel\), and over\-parameterized settings\. Solid lines are the median over33seeds; shaded bands the min–max\.

## 3Experimental Design and Evaluation Protocol

We test the claim that optimizer comparisons in matrix factorization are highly sensitive to hyperparameters, so that conclusions drawn from a single \(often default\) configuration are not representative\. Going beyond prior work\(ma2026preconditioningbenefitsspectralorthogonalization\), we evaluate five matrix\-recovery problems under a controlled protocol that holds model capacity and parameterization fixed and varies only optimization\-related hyperparameters, spanning well\-conditioned to severely*ill\-conditioned*targets\. Two further analyses isolate*how*conditioning acts: the sensitivity of each optimizer to the*shape*of the eigenvalue distribution at a fixed, extreme condition number \(Section[4\.1](https://arxiv.org/html/2607.13246#S4.SS1)\), and the dynamics with which each recovers the target’s spectral subspace \(Section[J](https://arxiv.org/html/2607.13246#A10)\)\.

### 3\.1Problem settings

All problems use ambient dimensiond=100d=100and a factorized parameterizationM^=U​V⊤\\hat\{M\}=UV^\{\\top\}withU,V∈ℝd×rU,V\\in\\mathbb\{R\}^\{d\\times r\}, optimized from a small random initialization\. We consider:

- •Low\-rank factorization \(conditioning\)\.A fully observed targetM⋆=U⋆​diag​\(σ\)​U⋆⊤M^\{\\star\}=U^\{\\star\}\\mathrm\{diag\}\(\\sigma\)U^\{\\star\\top\}of rankr=15r=15\(i\.e\., we solve for the symmetric factorization\) whose singular values are spaced linearly betweenκ\\kappaand11, soκ∈\{1,5,25,125,625\}\\kappa\\in\\\{1,5,25,125,625\\\}is exactly the condition number\. Objective14​∥U​U⊤−M⋆∥F2\\tfrac\{1\}\{4\}\\lVert UU^\{\\top\}\-M^\{\\star\}\\rVert\_\{F\}^\{2\}\(Figure[1](https://arxiv.org/html/2607.13246#S2.F1)\)\.
- •Matrix completion \(conditioning\)\.A rank\-44symmetric PSD target with the same conditioning sweep, observed on a uniformly random20%20\\%of entries; loss12​∥𝒫Ω​\(U​V⊤−M⋆\)∥F2\\tfrac\{1\}\{2\}\\lVert\\mathcal\{P\}\_\{\\Omega\}\(UV^\{\\top\}\-M^\{\\star\}\)\\rVert\_\{F\}^\{2\}, the squared error on the observed support \(Figure[7](https://arxiv.org/html/2607.13246#A3.F7)\)\.
- •Matrix completion \(over\-parameterization\)\.The same completion task at fixedκ=5\\kappa=5and true rank44, sweeping the search rankk∈\{4,5,100\}k\\in\\\{4,5,100\\\}\(matched, mildly, and heavily over\-parameterized\) \(Figure[8](https://arxiv.org/html/2607.13246#A3.F8)\)\.
- •Non\-negative factorization\.A non\-negative target with a projected \(⋅≥0\\,\\cdot\\,\\geq 0\) parameterization, at factor rankk∈\{10,50,100\}k\\in\\\{10,50,100\\\}, under two target spectra: a*uniform*spectrum and a*decayed*\(exponentially graded\) spectrum, the latter removing the flat “DC\-component” degeneracy of the uniform case \(Figure[2](https://arxiv.org/html/2607.13246#S2.F2)\)\.
- •Tensor\-train factorization \(over\-parameterization, noise\)\.A symmetric positive semidefinite targetM⋆=U⋆​U⋆⊤M^\{\\star\}=U^\{\\star\}U^\{\\star\\top\}recovered with untied coresU,VU,Vin a noiseless \(r⋆=4r^\{\\star\}=4\) and a noisy \(r⋆=30r^\{\\star\}=30\) regime, sweeping the search rank across under\-, exactly\-, and over\-parameterized values\. Objective12​∥U​V⊤−M⋆∥F2\\tfrac\{1\}\{2\}\\lVert UV^\{\\top\}\-M^\{\\star\}\\rVert\_\{F\}^\{2\}; see Section[4\.3](https://arxiv.org/html/2607.13246#S4.SS3)and Figures[5](https://arxiv.org/html/2607.13246#S2.F5)and[9](https://arxiv.org/html/2607.13246#A3.F9)\.

We compareMuon\(momentum0\.950\.95, Nesterov,55Newton–Schulz orthogonalization steps\) againstAdamW, plain gradient descent \(GD\), and sign gradient descent \(SignGD\)\.

### 3\.2Evaluation protocol

We evaluate recovery with the*normalized mean\-squared error*\(NMSE\), the squared reconstruction error relative to‖𝑴⋆‖F2\\left\\lVert\{\\bm\{M\}\}^\{\\star\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\}, which makes the metric scale\-free and comparable across problems and conditioning regimes\. For every \(problem, condition, optimizer\) triple we sweep the learning rate, the dominant axis and the usual source of misleading comparisons, over a logarithmically spaced grid \(2525points in\[10−4,5×10−1\]\[10^\{\-4\},5\\times 10^\{\-1\}\]for the factorization and completion tasks;2020points in\[10−4,3\.2\]\[10^\{\-4\},3\.2\]for NMF\), holding all other hyperparameters \(momentum, and forMuonthe Newton–Schulz countJJand the orthogonalization coefficients\) at standard values\. Each configuration runs for up to30003000iterations under an identical patience\-based learning\-rate\-decay schedule and is repeated over33seeds\. Because the resulting NMSE values span roughly thirty orders of magnitude, we summarize each configuration by the*geometric mean*over seeds with a±\\pmone\-log\-standard\-deviation band, and compare optimizers by their*stable learning\-rate range*, the band of step sizes over which each converges to the NMSE floor\.

## 4Experimental Results

A single learning rate is not a meaningful operating point\.Within a fixed problem and optimizer, the final loss varies by a median of about1111orders of magnitude across the learning\-rate grid for the factorization and completion tasks, and about66orders for both NMF variants \(Figures[1](https://arxiv.org/html/2607.13246#S2.F1)–[2](https://arxiv.org/html/2607.13246#S2.F2)\)\. Performance is therefore dominated by the learning rate, and any comparison made at one pre\-selected value reflects that choice as much as the optimizer itself\.

No shared learning rate is fair to all methods\.The loss\-minimizing learning rate differs*across optimizers*by a median of1\.51\.5–3\.13\.1decades \(largest for NMF\), so evaluating all methods at a common learning rate necessarily places at least one of them far from its own optimum\. This appears as a horizontal shift between each method’s basin in every panel\.

Tuning reorders the methods, and the winner is problem\-dependent\.Table[1](https://arxiv.org/html/2607.13246#A0.T1)reports each optimizer’s best \(learning\-rate\-tuned\) final loss\. In nearly all nineteen settings the ranking induced by a fixed “default” learning rate \(≈10−3\\approx 10^\{\-3\}\) differs from the ranking obtained once each method is tuned individually, and methods that look worst at the default frequently become best after tuning\. Under per\-optimizer tuning no method dominates universally: AdamW and GD win on plain low\-rank factorization, Muon and AdamW are comparable on matrix completion, and Muon wins clearly on both NMF variants\. The detailed per\-setting comparisons, including the specific default\-vs\-tuned reorderings, are deferred to Appendix[D](https://arxiv.org/html/2607.13246#A4)\.

Ill\-conditioning widens the gaps and shifts the winner\.Conditioning is itself a hyperparameter of the problem, and it interacts strongly with the optimizer\. Asκ\\kappagrows, every method’s attainable loss degrades and its optimal learning rate migrates upward \(Figures[1](https://arxiv.org/html/2607.13246#S2.F1),[7](https://arxiv.org/html/2607.13246#A3.F7)\)\. The effect is most dramatic for the non\-adaptive methods on factorization, where GD goes from the best\-tuned method in the well\-conditioned regime to failing under extreme ill\-conditioning while AdamW becomes best \(see Appendix[D](https://arxiv.org/html/2607.13246#A4)for the precise values\)\. Thus ill\-conditioning does not merely raise the loss floor, it reorders the methods, so that no fixed configuration summarizes behavior across the conditioning range\.

### 4\.1Sensitivity to spectrum shape under extreme ill\-conditioning

Conditioning alone does not determine difficulty: at a*fixed*condition number, the convergence of first\-order optimizers also depends on the*distribution*of eigenvalues within\[smin,smax\]\[s\_\{\\min\},s\_\{\\max\}\]\. We therefore hold the conditioning fixed at an extreme value and vary only the interior shape of the spectrum, evaluating seven spectral profiles \(defined in Appendix[F](https://arxiv.org/html/2607.13246#A6)\)\. All seven families share the same endpoints\(smin,smax\)=\(10−3,10\)\(s\_\{\\min\},s\_\{\\max\}\)=\(10^\{\-3\},10\)and hence the same, deliberately extreme, condition numberκ=smax/smin=104\\kappa=s\_\{\\max\}/s\_\{\\min\}=10^\{4\}; only the interior shape differs\. Despite this identical \(and severe\) ill\-conditioning, the distributions pose markedly different preconditioning challenges\. As shown in Figure[3](https://arxiv.org/html/2607.13246#S2.F3), AdamW adapts effectively to heavy\-tailed and clustered spectra, whereas Muon underperforms across the majority of these shapes, evidence that its spectral updates are less robust to specific eigenvalue densities even when the overall condition number is held fixed\. Ill\-conditioning is therefore not a single axis of difficulty: the sameκ\\kappacan be easy or hard depending on how the eigenvalues are arranged, and the two optimizers respond to that arrangement differently\.

### 4\.2Gaussian\-kernel factorization: a well\-behaved landscape

To complement the synthetic NMF problems, we study the Gaussian radial\-basis\-function kernel\. We sampleN=100N=100points inℝ5\\mathbb\{R\}^\{5\}and formKi​j=exp⁡\(−∥xi−xj∥2/2​σ2\)K\_\{ij\}=\\exp\(\-\\lVert x\_\{i\}\-x\_\{j\}\\rVert^\{2\}/2\\sigma^\{2\}\)withσ=2\\sigma=2, giving a symmetric positive\-semidefinite matrix with strictly positive entries and a rapidly decaying spectrum\. We factorize it asK^=U​V⊤\\hat\{K\}=UV^\{\\top\}under a non\-negativity constraint \(U,V≥0U,V\\geq 0, enforced by projection after each step\), minimizing12​∥U​V⊤−K∥F2\\tfrac\{1\}\{2\}\\lVert UV^\{\\top\}\-K\\rVert\_\{F\}^\{2\}\.

Rather than fix the factor rank arbitrarily, we set it by a*principled bottleneck*: the smallest rankrrwhose leading singular values capture at least95%95\\%ofKK’s spectral energy, which yieldsr=24r=24for this target\. Each optimizer is tuned by a learning\-rate sweep \(3030points in\[10−4,100\.5\]\[10^\{\-4\},10^\{0\.5\}\]\), and every configuration is repeated over three seeds that vary the initialization while holding the target fixed\.

The result is shown in Figure[10](https://arxiv.org/html/2607.13246#A3.F10)\. On this well\-behaved landscape, Muon’s advanced orthogonalization is unnecessary: standard AdamW performs just as effectively, converging seamlessly to the same global floor as Muon\. This reinforces the study’s central theme from the opposite direction, on the adversarial NMF spectra Muon’s structured updates help, but on a smooth, well\-conditioned kernel target that advantage disappears and a well\-tuned adaptive baseline is equally good\.

### 4\.3Tensor\-train factorization

We consider the tensor\-train \(TT\) factorization of a matrix target\. For an order\-22tensor𝑴⋆∈ℝd×d\{\\bm\{M\}\}^\{\\star\}\\in\\mathbb\{R\}^\{d\\times d\}, the TT decomposition contracts two cores𝒢1∈ℝd×r\\mathcal\{G\}\_\{1\}\\in\\mathbb\{R\}^\{d\\times r\}and𝒢2∈ℝr×d\\mathcal\{G\}\_\{2\}\\in\\mathbb\{R\}^\{r\\times d\}with boundary ranksr0=r2=1r\_\{0\}=r\_\{2\}=1and a single internal rankr1=rr\_\{1\}=r,

𝑴⋆​\(i,j\)=∑α=1r𝒢1​\(i,α\)​𝒢2​\(α,j\),⟺𝑴⋆=𝒢1​𝒢2,\{\\bm\{M\}\}^\{\\star\}\(i,j\)\\;=\\;\\sum\_\{\\alpha=1\}^\{r\}\\mathcal\{G\}\_\{1\}\(i,\\alpha\)\\,\\mathcal\{G\}\_\{2\}\(\\alpha,j\),\\qquad\\Longleftrightarrow\\qquad\{\\bm\{M\}\}^\{\\star\}\\;=\\;\\mathcal\{G\}\_\{1\}\\mathcal\{G\}\_\{2\},\(4\)so the TT factorization coincides with the bilinear factorization of Appendix[B](https://arxiv.org/html/2607.13246#A2)under𝑼=𝒢1\{\\bm\{U\}\}=\\mathcal\{G\}\_\{1\},𝑽⊤=𝒢2\{\\bm\{V\}\}^\{\\top\}=\\mathcal\{G\}\_\{2\}, with the lone TT\-rankrrplaying the role of the factorization rank\. In our setting the target𝑴⋆=𝑼⋆​𝑼⋆⊤\{\\bm\{M\}\}^\{\\star\}=\{\\bm\{U\}\}^\{\\star\}\{\\bm\{U\}\}^\{\\star\\top\}is symmetric positive semidefinite, but we do*not*tie the cores: recovery uses untied factors𝑼,𝑽\{\\bm\{U\}\},\{\\bm\{V\}\}, so the problem is the general bilinear factorization equation[12](https://arxiv.org/html/2607.13246#A2.E12)applied to a symmetric PSD target, retaining theGLr\\mathrm\{GL\}\_\{r\}imbalance invariance of equation[12](https://arxiv.org/html/2607.13246#A2.E12)\. Recovery minimizes

f​\(𝑼,𝑽\)=12​‖𝑼​𝑽⊤−𝑴⋆‖F2,𝑼∈ℝd×r,𝑽∈ℝd×r,f\(\{\\bm\{U\}\},\{\\bm\{V\}\}\)\\;=\\;\\frac\{1\}\{2\}\\,\\left\\lVert\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}\-\{\\bm\{M\}\}^\{\\star\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\},\\qquad\{\\bm\{U\}\}\\in\\mathbb\{R\}^\{d\\times r\},\{\\bm\{V\}\}\\in\\mathbb\{R\}^\{d\\times r\},\(5\)the objective of equation[12](https://arxiv.org/html/2607.13246#A2.E12)with search rankrras the TT\-rank\.

We study equation[5](https://arxiv.org/html/2607.13246#S4.E5)in two regimes, a*noiseless*regime \(true rankr⋆=4r^\{\\star\}=4, recovered exactly, loss limited only by convergence\) and a*noisy*regime \(true rankr⋆=30r^\{\\star\}=30with additive observation noise, loss floored at the statistical noise level\), and within each we sweep the search rank across under\-, exactly\-, and over\-parameterized values; full specifications are given in Appendix[G](https://arxiv.org/html/2607.13246#A7)\. The quantity of interest is the*stable learning\-rate range*of each optimizer, the band of step sizes over which equation[5](https://arxiv.org/html/2607.13246#S4.E5)converges to its regime\-dependent floor, and how that band widens, narrows, or shifts as the parameterization regime and the noise level change \(Figures[5](https://arxiv.org/html/2607.13246#S2.F5)and[9](https://arxiv.org/html/2607.13246#A3.F9)\)\. Beyond low\-rank recovery, this factorized objective serves as a minimal model of the product parameterizations ubiquitous in deep learning, a connection we develop in Appendix[H](https://arxiv.org/html/2607.13246#A8)\. We also study the effect of factorization depth in Appendix[I](https://arxiv.org/html/2607.13246#A9)\.

## 5Conclusion

In this work, we revisit optimizer comparisons for matrix factorization and show that conclusions are highly sensitive to hyperparameter selection and problem conditioning\. Systematic learning\-rate sweeps reveal that rankings reported under single configurations often do not hold: well\-tuned AdamW and gradient descent frequently match or outperform Muon on standard low\-rank factorization tasks, while Muon’s advantages are mostly limited to NMF and some highly ill\-conditioned completion problems\. We further find that optimizer performance depends on both conditioning severity and spectral structure, altering method rankings across problem instances\. These results underscore the need to evaluate optimizers across tuned hyperparameter ranges and conditioning regimes rather than relying on default settings or isolated benchmarks\.

## Acknowledgements

We thank Tianhao Wang for providing helpful references and for the discussions that assisted in developing the experimental setup\. Ali Parviz was funded by NSF CIF\-2403452\.

## References

Table 1:Learning\-rate\-tuned final loss \(geometric mean over33seeds; lower is better\)\. The best optimizer per row is inbold\. The tuned ranking varies across problems and, within each problem, differs from the ranking at a fixed default learning rate\.## Appendix ABackground: The Muon Optimizer

We briefly recall theMuonupdate and the spectral viewpoint that motivates it\. Consider a weight matrix𝑾∈ℝm×n\{\\bm\{W\}\}\\in\\mathbb\{R\}^\{m\\times n\}with loss gradient𝑮t=∇𝑾L​\(𝑾t\)\{\\bm\{G\}\}\_\{t\}=\\nabla\_\{\{\\bm\{W\}\}\}L\(\{\\bm\{W\}\}\_\{t\}\)at iterationtt\.Muonmaintains a momentum buffer𝑩t\{\\bm\{B\}\}\_\{t\}and applies an*orthogonalized*update:

𝑩t\\displaystyle\{\\bm\{B\}\}\_\{t\}=μ​𝑩t−1\+𝑮t,\\displaystyle=\\mu\\,\{\\bm\{B\}\}\_\{t\-1\}\+\{\\bm\{G\}\}\_\{t\},\(6\)𝑶t\\displaystyle\{\\bm\{O\}\}\_\{t\}=Newton​\-​Schulz⁡\(𝑩t\)≈polar⁡\(𝑩t\),\\displaystyle=\\operatorname\{Newton\\text\{\-\}Schulz\}\(\{\\bm\{B\}\}\_\{t\}\)\\;\\approx\\;\\operatorname\{polar\}\(\{\\bm\{B\}\}\_\{t\}\),\(7\)𝑾t\+1\\displaystyle\{\\bm\{W\}\}\_\{t\+1\}=𝑾t−α​𝑶t\.\\displaystyle=\{\\bm\{W\}\}\_\{t\}\-\\alpha\\,\{\\bm\{O\}\}\_\{t\}\.\(8\)Hereα\>0\\alpha\>0is the step size andμ∈\[0,1\)\\mu\\in\[0,1\)the momentum coefficient\. The orthogonalization step equation[7](https://arxiv.org/html/2607.13246#A1.E7)replaces the momentum direction by \(an approximation of\) its*polar factor*\.

Polar factor and spectral steepest descent\.Let𝑩t=𝑷​𝚺​𝑸⊤\{\\bm\{B\}\}\_\{t\}=\{\\bm\{P\}\}\{\\bm\{\\Sigma\}\}\{\\bm\{Q\}\}^\{\\top\}be a singular value decomposition, with𝑷∈𝒪m×ρ\{\\bm\{P\}\}\\in\\mathcal\{O\}\_\{m\\times\\rho\},𝑸∈𝒪n×ρ\{\\bm\{Q\}\}\\in\\mathcal\{O\}\_\{n\\times\\rho\}, and𝚺=𝖽𝗂𝖺𝗀​\{σ1,…,σρ\}\{\\bm\{\\Sigma\}\}=\\mathsf\{diag\}\\\{\\sigma\_\{1\},\\dots,\\sigma\_\{\\rho\}\\\}, whereρ=rank​\(𝑩t\)\\rho=\\mathrm\{rank\}\(\{\\bm\{B\}\}\_\{t\}\)\. The polar factor is

polar⁡\(𝑩t\)=𝑷​𝑸⊤,\\operatorname\{polar\}\(\{\\bm\{B\}\}\_\{t\}\)\\;=\\;\{\\bm\{P\}\}\{\\bm\{Q\}\}^\{\\top\},\(9\)which sets every nonzero singular value to one while preserving the singular vectors\. Equivalently,𝑷​𝑸⊤\{\\bm\{P\}\}\{\\bm\{Q\}\}^\{\\top\}is the solution of the spectral\-norm linear minimization oracle,

𝑷​𝑸⊤=arg​max‖𝑿‖op≤1⁡⟨𝑩t,𝑿⟩,\{\\bm\{P\}\}\{\\bm\{Q\}\}^\{\\top\}\\;=\\;\\operatorname\*\{arg\\,max\}\_\{\\left\\lVert\{\\bm\{X\}\}\\right\\rVert\_\{\\mathrm\{op\}\}\\leq 1\}\\left\\langle\{\\bm\{B\}\}\_\{t\},\{\\bm\{X\}\}\\right\\rangle,\(10\)so theMuondirection equation[7](https://arxiv.org/html/2607.13246#A1.E7)is precisely the steepest\-descent direction with respect to the spectral norm\. This is the sense in whichMuonis “spectrum aware”: it discards the magnitude of the momentum’s singular values and acts only along their directions\.

Newton–Schulz orthogonalization\.Computing the SVD at every step is expensive, soMuonapproximates the polar factor with a fixed number of matrix\-multiplication\-only Newton–Schulz iterations\. Starting from the normalized matrix𝑿0=𝑩t/‖𝑩t‖F\{\\bm\{X\}\}\_\{0\}=\{\\bm\{B\}\}\_\{t\}/\\left\\lVert\{\\bm\{B\}\}\_\{t\}\\right\\rVert\_\{\\mathrm\{F\}\}, one iterates a fixed odd polynomial

𝑿j\+1=a​𝑿j\+b​𝑿j​𝑿j⊤​𝑿j\+c​\(𝑿j​𝑿j⊤\)2​𝑿j,j=0,1,…,J−1,\{\\bm\{X\}\}\_\{j\+1\}\\;=\\;a\\,\{\\bm\{X\}\}\_\{j\}\\;\+\\;b\\,\{\\bm\{X\}\}\_\{j\}\{\\bm\{X\}\}\_\{j\}^\{\\top\}\{\\bm\{X\}\}\_\{j\}\\;\+\\;c\\,\\big\(\{\\bm\{X\}\}\_\{j\}\{\\bm\{X\}\}\_\{j\}^\{\\top\}\\big\)^\{2\}\{\\bm\{X\}\}\_\{j\},\\qquad j=0,1,\\dots,J\-1,\(11\)with coefficients\(a,b,c\)\(a,b,c\)chosen so that the induced scalar mapσ↦a​σ\+b​σ3\+c​σ5\\sigma\\mapsto a\\sigma\+b\\sigma^\{3\}\+c\\sigma^\{5\}pushes the singular values of𝑿j\{\\bm\{X\}\}\_\{j\}toward11\. The quintic coefficients\(a,b,c\)=\(3\.4445,−4\.7750,2\.0315\)\(a,b,c\)=\(3\.4445,\-4\.7750,2\.0315\)proposed byjordan2024muonreach a usable approximation in roughlyJ=5J=5iterations; subsequent work designs improved polynomials for this orthogonalization\(amsel2025polarexpress;grishina2025chebyshevns;cesista2025muonoptcoeffs;boumal2025polarpoly\)\. AsJ→∞J\\to\\infty\(with an exact map\) the iteration converges topolar⁡\(𝑩t\)\\operatorname\{polar\}\(\{\\bm\{B\}\}\_\{t\}\), recovering the idealized update; in practiceJJis small and𝑶t\{\\bm\{O\}\}\_\{t\}is only an approximate polar factor, which several recent analyses account for explicitly\(shulgin2025beyondideal;anon2025newtonschulz;lau2025polargrad\)\.

## Appendix BMatrix Factorization Variants

This appendix collects two generalizations of the symmetric matrix factorization problem equation[1](https://arxiv.org/html/2607.13246#S2.E1): the general bilinear \(asymmetric\) factorization and the matrix completion setting in which only a subset of the entries of the target is observed\.

##### General bilinear factorization\.

The symmetric problem equation[1](https://arxiv.org/html/2607.13246#S2.E1)is a special case of the general*bilinear*\(asymmetric\) factorization, in which a possibly rectangular target𝑴⋆∈ℝd1×d2\{\\bm\{M\}\}^\{\\star\}\\in\\mathbb\{R\}^\{d\_\{1\}\\times d\_\{2\}\}of rankrris recovered from the product of two factors,

min𝑼∈ℝd1×k𝑽∈ℝd2×kf​\(𝑼,𝑽\)=12​‖𝑼​𝑽⊤−𝑴⋆‖F2,\\min\_\{\\begin\{subarray\}\{c\}\{\\bm\{U\}\}\\in\\mathbb\{R\}^\{d\_\{1\}\\times k\}\\\\\[2\.0pt\] \{\\bm\{V\}\}\\in\\mathbb\{R\}^\{d\_\{2\}\\times k\}\\end\{subarray\}\}\\quad f\(\{\\bm\{U\}\},\{\\bm\{V\}\}\)\\;=\\;\\frac\{1\}\{2\}\\,\\left\\lVert\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}\-\{\\bm\{M\}\}^\{\\star\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\},\(12\)withk≥rk\\geq r\. When𝑴⋆\{\\bm\{M\}\}^\{\\star\}is symmetric positive semidefinite and the factors are tied as𝑽=𝑼\{\\bm\{V\}\}=\{\\bm\{U\}\}, equation[12](https://arxiv.org/html/2607.13246#A2.E12)reduces to equation[1](https://arxiv.org/html/2607.13246#S2.E1)up to a constant rescaling of the objective\. Unlike the symmetric case, the asymmetric parameterization is invariant under the larger group of invertible transformations\(𝑼,𝑽\)↦\(𝑼​𝑹,𝑽​𝑹−⊤\)\(\{\\bm\{U\}\},\{\\bm\{V\}\}\)\\mapsto\(\{\\bm\{U\}\}\{\\bm\{R\}\},\{\\bm\{V\}\}\{\\bm\{R\}\}^\{\-\\top\}\)for any𝑹∈GLk​\(ℝ\)\{\\bm\{R\}\}\\in\\mathrm\{GL\}\_\{k\}\(\\mathbb\{R\}\), since𝑼​𝑽⊤=\(𝑼​𝑹\)​\(𝑽​𝑹−⊤\)⊤\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}=\(\{\\bm\{U\}\}\{\\bm\{R\}\}\)\(\{\\bm\{V\}\}\{\\bm\{R\}\}^\{\-\\top\}\)^\{\\top\}\. This invariance permits an arbitrary imbalance between the factor norms and is commonly controlled by a balancing regularizer12​\(‖𝑼‖F2\+‖𝑽‖F2\)\\tfrac\{1\}\{2\}\\bigl\(\\left\\lVert\{\\bm\{U\}\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\}\+\\left\\lVert\{\\bm\{V\}\}\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\}\\bigr\)or by enforcing𝑼⊤​𝑼=𝑽⊤​𝑽\{\\bm\{U\}\}^\{\\top\}\{\\bm\{U\}\}=\{\\bm\{V\}\}^\{\\top\}\{\\bm\{V\}\}along the trajectory\.

##### Matrix completion\.

In many applications only a subset of the entries of𝑴⋆\{\\bm\{M\}\}^\{\\star\}is observed\. Let\[n\]≔\{1,…,n\}\[n\]\\coloneqq\\\{1,\\dots,n\\\}and letΩ⊆\[d1\]×\[d2\]\\Omega\\subseteq\[d\_\{1\}\]\\times\[d\_\{2\}\]denote the set of observed indices\. Define the sampling operator𝒫Ω:ℝd1×d2→ℝd1×d2\\mathcal\{P\}\_\{\\Omega\}:\\mathbb\{R\}^\{d\_\{1\}\\times d\_\{2\}\}\\to\\mathbb\{R\}^\{d\_\{1\}\\times d\_\{2\}\}by

\[𝒫Ω​\(𝑿\)\]i​j=\{Xi​j,\(i,j\)∈Ω,0,otherwise,\\bigl\[\\mathcal\{P\}\_\{\\Omega\}\(\{\\bm\{X\}\}\)\\bigr\]\_\{ij\}\\;=\\;\\begin\{cases\}X\_\{ij\},&\(i,j\)\\in\\Omega,\\\\\[2\.0pt\] 0,&\\text\{otherwise,\}\\end\{cases\}\(13\)which retains the observed entries and zeros out the rest\. The factored matrix completion problem is

min𝑼∈ℝd1×k𝑽∈ℝd2×k12​‖𝒫Ω​\(𝑼​𝑽⊤−𝑴⋆\)‖F2\.\\min\_\{\\begin\{subarray\}\{c\}\{\\bm\{U\}\}\\in\\mathbb\{R\}^\{d\_\{1\}\\times k\}\\\\\[2\.0pt\] \{\\bm\{V\}\}\\in\\mathbb\{R\}^\{d\_\{2\}\\times k\}\\end\{subarray\}\}\\quad\\frac\{1\}\{2\}\\,\\left\\lVert\\mathcal\{P\}\_\{\\Omega\}\\\!\\bigl\(\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}\-\{\\bm\{M\}\}^\{\\star\}\\bigr\)\\right\\rVert\_\{\\mathrm\{F\}\}^\{2\}\.\(14\)TakingΩ=\[d1\]×\[d2\]\\Omega=\[d\_\{1\}\]\\times\[d\_\{2\}\]recovers the fully observed problem equation[12](https://arxiv.org/html/2607.13246#A2.E12), while the symmetric variant tied to equation[1](https://arxiv.org/html/2607.13246#S2.E1)additionally imposes𝑽=𝑼\{\\bm\{V\}\}=\{\\bm\{U\}\}\. Exact recovery in this setting requires the ground\-truth factors to be incoherent and the sample size\|Ω\|\\lvert\\Omega\\rvertto be sufficiently large relative to the degrees of freedomr​\(d1\+d2−r\)r\(d\_\{1\}\+d\_\{2\}\-r\)\. Under the Bernoulli model in which each entry is observed independently with probabilitypp, the rescaled operatorp−1​𝒫Ωp^\{\-1\}\\mathcal\{P\}\_\{\\Omega\}is an unbiased estimator of the identity, i\.e\.𝔼​\[p−1​𝒫Ω​\(𝑿\)\]=𝑿\\mathbb\{E\}\\bigl\[p^\{\-1\}\\mathcal\{P\}\_\{\\Omega\}\(\{\\bm\{X\}\}\)\\bigr\]=\{\\bm\{X\}\}\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x5.png)Figure 6:Spectral subspace recovery: Muon vs\. AdamW\.Absolute alignmentMi​j=\|⟨ui,ej⟩\|M\_\{ij\}=\\lvert\\langle u\_\{i\},e\_\{j\}\\rangle\\rvert\(Eq\.[16](https://arxiv.org/html/2607.13246#A10.E16)\) between learned singular vectorsuiu\_\{i\}and true eigenvectorseje\_\{j\}ofKK, every500500steps \(AdamW top, Muon bottom\)\. A bright, correctly\-ordered diagonal indicates faithful recovery; off\-diagonal mass reflects mixing\. AdamW sharpens its diagonal steadily, while Muon stays near\-random for∼1500\{\\sim\}1500steps, then snaps into a more diffuse alignment\.

## Appendix CRelated Work

Spectrum\-aware optimization andMuon\.Much of the theoretical analysis ofMuonbuilds on a perspective that interprets modern optimizers such asAdamandShampooas instances of steepest descent under non\-Euclidean or norm\-constrained geometries\(bernstein2024old\)\. Within this view,Muonapproximates a spectral\-norm steepest\-descent direction via polar factorization, connecting it to a broader class of methods that leverage matrix structure in gradients\(pethick2025normconstrainedlmo;li2025noteconvergencemuon;shen2025convergencemuon;chen2025muonspectralnorm\)\. Spectral transformations and orthogonalization in optimization predateMuon, appearing in earlier work on stochastic and preconditioned methods\(carlson2015stochasticB;carlson2015stochastic;carlson2015preconditioned;tuddenham2022orthogonalising\)\. A complementary line of work asks more directly when such spectral updates are beneficial in deep learning\(davis2025whenspectral;su2025isotropic\)\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x6.png)Figure 7:Matrix completion, conditioning sweep\.Each panel plots the final reconstruction loss \(yy\-axis, log scale; geometric mean over33seeds with a shaded±\\pmone\-log\-standard\-deviation band\) against the learning rate \(xx\-axis, log scale\), for condition numberκ\\kappa\. Muon and AdamW are indistinguishable in the well\-conditioned regime; all methods degrade and converge toward one another asκ\\kappagrows\.![Refer to caption](https://arxiv.org/html/2607.13246v1/x7.png)Figure 8:Matrix completion, search\-rank sweep\(true rank44,κ=5\\kappa\{=\}5\)\. Over\-parameterization \(rank100100\) widens the band of effective learning rates and lets every method except GD reach machine precision\.![Refer to caption](https://arxiv.org/html/2607.13246v1/x8.png)Figure 9:Tensor\-Train Noisy regime\(true rankr⋆=30r^\{\\star\}=30, additive observation noise\)\. Final loss versus learning rate for each optimizer, across search ranksr∈\{5,15,30,60\}r\\in\\\{5,15,30,60\\\}spanning under\-, exactly\- \(r=r⋆r=r^\{\\star\}, bold panel\), and over\-parameterized settings\. Solid lines are the median over33seeds; shaded bands the min–max; the loss floors at the noise level rather than at zero\.Implicit bias and generalization\.A second line of work studies the implicit bias induced byMuon\.fan2025implicitshow that in linear classification, idealizedMuonconverges to a solution that maximizes margin with respect to the spectral norm, contrasting with the Euclidean and coordinate\-wise biases ofSGDandAdam\. Complementary empirical studies suggest that spectrum\-aware updates can improve generalization, particularly in imbalanced or long\-tailed settings, by promoting more uniform learning across principal components rather than focusing on dominant directions\(vasudeva2025muon;wang2025muon\)\. Related analyses indicate thatMuonyields more isotropic singular\-value spectra thanAdamand may accelerate phenomena such as grokking in long\-horizon training\(zhang2025concurrence;wang2025muon;vasudeva2025muon;tveit2025muon\)\. Connections have also been drawn betweenMuonand second\-order or preconditioned methods, with several works interpreting its updates as approximations toShampoo\(gupta2018shampoo;jordan2024muon;shah2025practical\)\. Alternative derivations and closely related formulations, some predatingMuon, have appeared from different theoretical viewpoints\(pethick2025training;carlson2015preconditioned;lau2025polargrad;bernstein2024modular;bernstein2024old;an2025asgo\)\.

Preconditioning for matrix factorization\.A separate body of work shows that preconditioning can dramatically accelerate optimization in matrix factorization\. In particular,ScaledGDand its variants achieve linear convergence rates independent of the condition number under appropriate initialization, in both exactly parameterized and over\-parameterized regimes\(tong2021accelerating;tong2021low;zhang2021preconditioned;zhang2023preconditioned\), and these guarantees extend to small random initialization\(xu2023power\)\. Closest to our motivation, recent work analyzes the preconditioning effect ofMuon’s spectral orthogonalization directly\(ma2026preconditioningbenefitsspectralorthogonalization\)\.

Summary\.Taken together, these results show that spectrum\-aware and preconditioned methods can offer strong theoretical and empirical benefits, but that their behavior depends heavily on problem structure and initialization\. This motivates a closer examination ofMuonin controlled settings such as matrix factorization, where powerful alternatives already enjoy strong guarantees and where the advantage ofMuonover adaptive methods likeAdamWis not immediately evident\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x9.png)Figure 10:Gaussian\-kernel NMF at a principled bottleneck\.Non\-negative factorization of a fixed Gaussian RBF kernelK∈ℝ100×100K\\in\\mathbb\{R\}^\{100\\times 100\}at search rankr=24r\{=\}24\(the smallest rank capturing95%95\\%ofKK’s spectral energy\)\. Tuned final loss vs\. learning rate \(geometric mean over33seeds,±\\pmlog\-std band; patience\-based decay schedule\)\. On this smooth, well\-behaved landscape Muon’s orthogonalized updates give no advantage: AdamW reaches the same global floor as Muon\.
## Appendix DDetailed experimental results

This appendix expands the summary claims of Section[4](https://arxiv.org/html/2607.13246#S4)with the specific numerical comparisons\. All values refer to Table[1](https://arxiv.org/html/2607.13246#A0.T1)\.

##### Horizontal basin shifts\.

Evaluating all methods at a common learning rate necessarily places at least one of them far from its own optimum\. This is visible as a horizontal shift between each method’s basin in every panel: for instance, GD’s basin lies near the high end of the grid while Muon’s and AdamW’s lie one to three decades lower\.

##### Default\-vs\-tuned reordering\.

The ranking induced by a fixed default learning rate \(≈10−3\\approx 10^\{\-3\}\) differs from the tuned ranking in nearly all nineteen settings, and methods that look worst at the default frequently become best after tuning\. On low\-rank factorization GD is among the worst at the default rate yet attains the lowest loss of all methods once tuned \(Figure[1](https://arxiv.org/html/2607.13246#S2.F1)\); symmetrically, on both NMF variants Muon is*last*at the default rate but*first*once tuned \(Figures[2](https://arxiv.org/html/2607.13246#S2.F2)\)\.

##### Effect of ill\-conditioning\.

The effect of ill\-conditioning is most dramatic for the non\-adaptive methods on factorization: GD, the best\-tuned method in the well\-conditioned regime \(1\.6×10−131\.6\\times 10^\{\-13\}atκ=1\\kappa\{=\}1\), collapses to1\.3×10−11\.3\\times 10^\{\-1\}atκ=625\\kappa\{=\}625, where AdamW becomes best \(2\.6×10−72\.6\\times 10^\{\-7\}\), Muon second \(9\.4×10−59\.4\\times 10^\{\-5\}\), and SignGD diverges\. Matrix completion shows the same upward drift of the loss floor: atκ=625\\kappa\{=\}625all four methods are confined toO​\(10−2\)O\(10^\{\-2\}\)and become nearly indistinguishable\.

##### Problem\-dependence of the tuned winner\.

Under per\-optimizer tuning there is no universally dominant method\. On*plain low\-rank factorization*, AdamW and GD drive the loss to10−1310^\{\-13\}–10−910^\{\-9\}across the whole conditioning range, while the Muon configuration plateaus several orders higher \(10−810^\{\-8\}–10−510^\{\-5\}\)\. On*matrix completion*, Muon and AdamW are both reaching10−1610^\{\-16\}–10−1310^\{\-13\}in the well\-conditioned regime, and over\-parameterization \(search rank100100\) lets all methods except GD reach machine precision \(Figure[8](https://arxiv.org/html/2607.13246#A3.F8)\)\. On*NMF*, the picture reverses: Muon is the clear winner under both spectra, reaching10−910^\{\-9\}–10−610^\{\-6\}while AdamW stalls near10−310^\{\-3\}\(with the single exception of the uniform\-spectrum rank\-5050case, where AdamW matches Muon at∼4×10−9\\sim 4\\times 10^\{\-9\}\) and GD/SignGD fail to fit\. Whichever optimizer a single\-configuration study would crown thus depends entirely on the problem and the learning rate chosen\.

## Appendix EDynamics of spectral subspace recovery

The aggregate losses above tell us*whether*a method fits the target but not*how*it does so\. To probe the mechanism, we factorize a fixed Gaussian\-kernel targetK∈ℝ100×100K\\in\\mathbb\{R\}^\{100\\times 100\}, symmetric positive semidefinite, with eigendecompositionK=Q​Λ​Q⊤K=Q\\Lambda Q^\{\\top\}and eigenvectors ordered by descending eigenvalue, asK^=U​V⊤\\hat\{K\}=UV^\{\\top\}at full search rankk=100k=100, training AdamW and Muon each at its own tuned learning rate\. Every500500iterations we take the left singular vectors\{u^i\}\\\{\\hat\{u\}\_\{i\}\\\}of the current reconstructionU​V⊤UV^\{\\top\}and form the*alignment matrix*

Mi​j=\|⟨u^i,qj⟩\|,M\_\{ij\}\\;=\\;\\bigl\|\\langle\\hat\{u\}\_\{i\},\\;q\_\{j\}\\rangle\\bigr\|,\(15\)the absolute overlap between theii\-th learned singular direction and thejj\-th true eigenvector\. A faithful, correctly ordered recovery of the eigenbasis yieldsM→IM\\to I\(a sharp diagonal\), whereas off\-diagonal or diffuse mass indicates that the learned subspace mixes or reorders eigendirections\. Figure[6](https://arxiv.org/html/2607.13246#A2.F6)tracksMMover training\. AdamW sharpens toward a clean diagonal quickly and monotonically, recovering the eigenvectors faithfully and in order, while Muon’s alignment remains comparatively diffuse, with persistent off\-diagonal mass and slower diagonalization\. This provides a mechanistic explanation for the loss\-level gap of Sections[4](https://arxiv.org/html/2607.13246#S4)–[4\.1](https://arxiv.org/html/2607.13246#S4.SS1): on these structured and ill\-conditioned spectra, AdamW’s per\-coordinate adaptation locks onto the dominant eigendirections in the right order, whereas Muon’s orthogonalized updates distribute capacity more evenly across directions and recover the eigenbasis less cleanly\. The heatmaps shown are for a representative tuned run; the qualitative pattern is consistent across seeds\.

## Appendix FSpectral profile definitions

The experiments of Section[4\.1](https://arxiv.org/html/2607.13246#S4.SS1)use seven spectral profiles, defined by the diagonal entries ofSS\. All share the endpoints\(smin,smax\)=\(10−3,10\)\(s\_\{\\min\},s\_\{\\max\}\)=\(10^\{\-3\},10\)and hence the condition numberκ=smax/smin=104\\kappa=s\_\{\\max\}/s\_\{\\min\}=10^\{4\}; only the interior shape differs\.

- •Uniform:linear spacingsi=smin\+i−1n−1​\(smax−smin\)s\_\{i\}=s\_\{\\min\}\+\\frac\{i\-1\}\{n\-1\}\(s\_\{\\max\}\-s\_\{\\min\}\), a dense, even spread of curvature\.
- •Log\-uniform:geometric spacingsi=smin⋅κi−1n−1s\_\{i\}=s\_\{\\min\}\\cdot\\kappa^\{\\frac\{i\-1\}\{n\-1\}\}, typical of ill\-conditioned systems whose scales vary by orders of magnitude\.
- •Gaussian:eigenvalues drawn from a truncated normal𝒩​\(μ,σ2\)\\mathcal\{N\}\(\\mu,\\sigma^\{2\}\)and rescaled to\[smin,smax\]\[s\_\{\\min\},s\_\{\\max\}\], a bulk concentrated around a mean curvature\.
- •Spiked:a single top eigenvalue atsmaxs\_\{\\max\}with the remainingn−1n\-1clustered nearsmins\_\{\\min\}, a dominant rank\-one subspace\.
- •Bimodal:eigenvalues split into two clusters near the endpointssmins\_\{\\min\}andsmaxs\_\{\\max\}\.
- •Polynomial decay:si∝i−αs\_\{i\}\\propto i^\{\-\\alpha\}, a rapid decay in the importance of higher\-order directions\.
- •Exponential decay:si∝e−β​is\_\{i\}\\propto e^\{\-\\beta i\}, an aggressive concentration of energy in the top eigenvalues\.

## Appendix GTensor\-train regimes and parameterization sweep

This appendix gives the full specification of the regimes and parameterization sweep summarized in Section[4\.3](https://arxiv.org/html/2607.13246#S4.SS3)\.

##### Regimes\.

We study equation[5](https://arxiv.org/html/2607.13246#S4.E5)in two regimes\. In the*noiseless*regime the target has true rankr⋆=4r^\{\\star\}=4and is recovered exactly, so the attainable loss is limited only by the optimizer’s convergence; in the*noisy*regime the target has true rankr⋆=30r^\{\\star\}=30and is corrupted by additive observation noise, so the loss is floored at the statistical noise level rather than driven to zero\. The two regimes therefore probe complementary phenomena: optimization geometry in the well\-specified case and robustness to a nonzero residual in the misspecified case\.

##### Parameterization sweep\.

Within each regime we sweep the search rank, equivalently the TT\-rankrrof the factor𝑼\{\\bm\{U\}\}, across under\-, exactly\-, and over\-parameterized values:r∈\{2,4,8,20\}r\\in\\\{2,4,8,20\\\}for the noiseless target \(true rankr⋆=4r^\{\\star\}=4\) andr∈\{5,15,30,60\}r\\in\\\{5,15,30,60\\\}for the noisy target \(true rankr⋆=30r^\{\\star\}=30\), withr=r⋆r=r^\{\\star\}the exactly parameterized boundary case\. Under\-parameterization \(r<r⋆r<r^\{\\star\}\) caps the attainable loss because𝑼​𝑼⊤\{\\bm\{U\}\}\{\\bm\{U\}\}^\{\\top\}cannot represent𝑴⋆\{\\bm\{M\}\}^\{\\star\}, while over\-parameterization \(r\>r⋆r\>r^\{\\star\}\) enlarges the factor space and introduces additional flat directions in the landscape\.

## Appendix HRelation to deep learning models

The factorized objective equation[5](https://arxiv.org/html/2607.13246#S4.E5)is not only a tensor\-train problem in its own right; it is the minimal instance of the product parameterizations that pervade deep learning, which is what makes the optimizer behavior we study here relevant beyond low\-rank recovery\.

##### Factorization as a shallow linear network\.

The map𝑼↦𝑼​𝑼⊤\{\\bm\{U\}\}\\mapsto\{\\bm\{U\}\}\{\\bm\{U\}\}^\{\\top\}\(or\(𝑼,𝑽\)↦𝑼​𝑽⊤\(\{\\bm\{U\}\},\{\\bm\{V\}\}\)\\mapsto\{\\bm\{U\}\}\{\\bm\{V\}\}^\{\\top\}in the asymmetric case of Appendix[B](https://arxiv.org/html/2607.13246#A2)\) is exactly a two\-layer*linear*network with no intervening nonlinearity: the factors are the layer weights, and the product is the end\-to\-end map\. Consequently the loss is nonconvex in the factors despite being convex in the product, and the search rank, equivalently the TT\-rank, plays the role of the network width\. The over\-parameterized regimer\>r⋆r\>r^\{\\star\}is precisely the width\-overparameterized regime of modern networks, and the flat directions it introduces in the landscape are the source of the implicit bias of gradient methods studied for matrix and deep factorizations\(gunasekar2017implicit;arora2019implicit\)\. This is why over\-parameterized factorization is a standard proxy for the optimization of deep models: it isolates the effect of redundant parameters on the trajectory while remaining analytically tractable\.

##### Depth as a longer tensor\-train chain\.

Increasing the order of the tensor\-train chain, contractingLLcores𝒢1,…,𝒢L\\mathcal\{G\}\_\{1\},\\dots,\\mathcal\{G\}\_\{L\}rather than two, is the algebraic analogue of increasing the depth of a linear network, whose end\-to\-end map factorizes as𝑾=𝑾L​𝑾L−1​⋯​𝑾1\{\\bm\{W\}\}=\{\\bm\{W\}\}\_\{L\}\{\\bm\{W\}\}\_\{L\-1\}\\cdots\{\\bm\{W\}\}\_\{1\}\. The matrix case equation[4](https://arxiv.org/html/2607.13246#S4.E4)\(L=2L=2\) is the shallowest such chain, and depth introduces exactly the ill\-conditioning and balancing phenomena that motivate our conditioning sweep: the product of many factors amplifies spectral imbalance, and the overparameterized invariance group acts on each internal bond\. Depth has been shown to act as an implicit preconditioner that accelerates gradient descent on these product objectives\(arora2018optimization\), so the sensitivity of each optimizer to conditioning atL=2L=2is the base case of a phenomenon that sharpens with depth\.

##### Tensor\-train structure inside trained networks\.

Beyond serving as a proxy, tensor\-train factorizations appear directly inside deep models\. Reshaping a dense weight matrix into a high\-order tensor and representing it in TT format compresses fully connected and embedding layers by orders of magnitude while preserving accuracy\(novikov2015tensorizing\), and specific architectures correspond to specific tensor decompositions: recurrent networks realize the tensor\-train / matrix\-product structure\(khrulkov2018expressive\), while convolutional arithmetic circuits realize the hierarchical Tucker decomposition\(cohen2016expressive\)\. In all of these the rank of the decomposition controls expressivity exactly as the TT\-rankrrcontrols the capacity of equation[5](https://arxiv.org/html/2607.13246#S4.E5), so the interaction between search rank, conditioning, and optimizer that we characterize on the factorized problem speaks directly to how these layers are trained\.

## Appendix IOptimizer Comparison Across Tensor\-Train Depth

We benchmark four optimizers—Muon, AdamW, gradient descent \(GD\), and SignGD—on non\-negative tensor\-train \(TT\) factorization\. To systematically vary problem conditioning, we use the tensor\-train depth \(equivalently, the tensor order or number of TT cores\) as the primary experimental factor\. Each TT core is represented via a matrix unfolding, a necessary implementation detail because Muon’s Newton–Schulz orthogonalization operates only on matrix\-valued gradients; retaining a core in its native third\-order form would effectively reduce Muon to momentum SGD on that factor\. We evaluate two target regimes: a clean, low\-rank, noise\-free tensor and a noisy, higher\-rank tensor\. The study spans depths \(2\)–\(6\), an overparameterized factorization rank, a five\-point learning\-rate sweep, and three random seeds, yielding 600 total runs\. Performance is reported as the best\-tuned final loss for each optimizer, defined as the lowest mean squared error \(MSE\) achieved across the learning\-rate grid and averaged over seeds\.

The results differ markedly between the two regimes\. For the clean target, optimizer rankings change substantially with depth\. AdamW achieves the lowest losses at shallow depths, while Muon performs comparatively poorly\. However, beginning at depth \(4\), the performance of AdamW, GD, and SignGD deteriorates as the factorization deepens, whereas Muon maintains consistently low error\. By depth \(6\), Muon outperforms the competing methods by approximately one to two orders of magnitude \(Figure[11](https://arxiv.org/html/2607.13246#A9.F11), left\)\. In contrast, performance on the noisy target remains relatively flat across depths, with all optimizers achieving similar losses \(Figure[11](https://arxiv.org/html/2607.13246#A9.F11), right\)\. The added noise introduces an irreducible error floor that limits the optimization gains obtainable from improved conditioning, thereby diminishing Muon’s advantage\.

These findings suggest that Muon’s benefits arise primarily from its ability to mitigate optimization difficulties associated with ill\-conditioned deep factorizations, rather than from a general robustness advantage\. At the same time, two limitations should be noted\. First, increasing tensor\-train depth changes the underlying optimization problem itself, so the observed trends reflect changes in relative optimizer performance across a family of increasingly difficult tasks rather than within a fixed loss landscape\. Second, the reported results correspond only to the overparameterized setting considered here and may not fully characterize behavior at other scales or rank regimes\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x10.png)Figure 11:Best\-tuned reconstruction loss versus tensor\-train depth\.Each curve shows the final MSE of one optimizer at its best learning rate \(minimum over the grid\) for a given depth, averaged over three seeds; shaded bands span the per\-seed min–max\.Left:clean, low\-rank target\.Right:noisy, higher\-rank target\. Search bond rank is fixed to the over\-parameterized regime \(R=12R\{=\}12clean,R=30R\{=\}30noisy\)\. Lower is better; theyy\-axis is logarithmic\.
## Appendix JDynamics of spectral subspace recovery

The aggregate losses above tell us*whether*a method fits the target but not*how*it does so\. To probe the mechanism, we factorize a fixed Gaussian\-kernel targetK∈ℝ100×100K\\in\\mathbb\{R\}^\{100\\times 100\}, symmetric positive semidefinite, with eigendecompositionK=Q​Λ​Q⊤K=Q\\Lambda Q^\{\\top\}and eigenvectors ordered by descending eigenvalue, asK^=U​V⊤\\hat\{K\}=UV^\{\\top\}at full search rankk=100k=100, training AdamW and Muon each at its own tuned learning rate\. Every500500iterations we take the left singular vectors\{u^i\}\\\{\\hat\{u\}\_\{i\}\\\}of the current reconstructionU​V⊤UV^\{\\top\}and form the*alignment matrix*

Mi​j=\|⟨u^i,qj⟩\|,M\_\{ij\}\\;=\\;\\bigl\|\\langle\\hat\{u\}\_\{i\},\\;q\_\{j\}\\rangle\\bigr\|,\(16\)the absolute overlap between theii\-th learned singular direction and thejj\-th true eigenvector\. A faithful, correctly ordered recovery of the eigenbasis yieldsM→IM\\to I\(a sharp diagonal\), whereas off\-diagonal or diffuse mass indicates that the learned subspace mixes or reorders eigendirections\.

Figure[6](https://arxiv.org/html/2607.13246#A2.F6)tracksMMover training\. AdamW sharpens toward a clean diagonal quickly and monotonically, recovering the eigenvectors faithfully and in order, while Muon’s alignment remains comparatively diffuse, with persistent off\-diagonal mass and slower diagonalization\. This provides a mechanistic explanation for the loss\-level gap of Sections[4](https://arxiv.org/html/2607.13246#S4)–[4\.1](https://arxiv.org/html/2607.13246#S4.SS1): on these structured and ill\-conditioned spectra, AdamW’s per\-coordinate adaptation locks onto the dominant eigendirections in the right order, whereas Muon’s orthogonalized updates distribute capacity more evenly across directions and recover the eigenbasis less cleanly\. The heatmaps shown are for a representative tuned run; the qualitative pattern is consistent across seeds\.

## Appendix KSoftplus, Its Taylor Surrogate, and Nonlinear Matrix Factorization

LetX∈ℝN×dX\\in\\mathbb\{R\}^\{N\\times d\}collect theNNinputs\. We fit the Gaussian kernelK∈ℝN×NK\\in\\mathbb\{R\}^\{N\\times N\},Ki​j=exp⁡\(−∥xi−xj∥2/2​σ2\)K\_\{ij\}=\\exp\(\-\\lVert x\_\{i\}\-x\_\{j\}\\rVert^\{2\}/2\\sigma^\{2\}\), with the modelK^=ϕ​\(X​W\)​V⊤\\hat\{K\}=\\phi\(XW\)\\,V^\{\\top\}, whereW∈ℝd×RW\\in\\mathbb\{R\}^\{d\\times R\},V∈ℝN×RV\\in\\mathbb\{R\}^\{N\\times R\}, the search rankRRis the hidden width, andϕ\\phiacts entrywise\. This is a*nonlinear*matrix factorization: forϕ=id\\phi=\\mathrm\{id\}it collapses to the bilinear low\-rank modelK^=X​\(W​V⊤\)\\hat\{K\}=X\(WV^\{\\top\}\)\.

##### Softplus and its expansion\.

The softplusϕ​\(z\)=log⁡\(1\+ez\)\\phi\(z\)=\\log\(1\+e^\{z\}\)has derivativeϕ′=σ\\phi^\{\\prime\}=\\sigma\(the logistic sigmoid\)\. Sinceσ−12\\sigma\-\\tfrac\{1\}\{2\}is odd, the Maclaurin series ofϕ\\phikeeps only the constant, linear, and even\-order terms,

ϕ​\(z\)=log⁡2\+12​z\+18​z2−1192​z4\+𝒪​\(z6\)\.\\phi\(z\)=\\log 2\+\\tfrac\{1\}\{2\}z\+\\tfrac\{1\}\{8\}z^\{2\}\-\\tfrac\{1\}\{192\}z^\{4\}\+\\mathcal\{O\}\(z^\{6\}\)\.\(17\)The surrogate used in the experiments,ϕ~​\(z\)=log⁡2\+12​z\+18​z2\\tilde\{\\phi\}\(z\)=\\log 2\+\\tfrac\{1\}\{2\}z\+\\tfrac\{1\}\{8\}z^\{2\}, is the truncation after the quadratic term \(the cubic coefficient vanishes\), and is accurate when the pre\-activationsX​WXWare small, as enforced by the input scale\.

##### The quadratic surrogate is an explicit feature factorization\.

WithA=ϕ~​\(X​W\)A=\\tilde\{\\phi\}\(XW\), expanding term by term gives

K^=A​V⊤=log⁡2​1N×R​V⊤⏟constant​\(rank≤1\)\+12​\(X​W\)​V⊤⏟rank≤min⁡\(R,d\)\+18​\(X​W\)⊙2​V⊤⏟degree\-2 lift,\\hat\{K\}=AV^\{\\top\}=\\underbrace\{\\log 2\\,\\mathbf\{1\}\_\{N\\times R\}V^\{\\top\}\}\_\{\\text\{constant \}\(\\mathrm\{rank\}\\leq 1\)\}\\;\+\\;\\underbrace\{\\tfrac\{1\}\{2\}\\,\(XW\)V^\{\\top\}\}\_\{\\mathrm\{rank\}\\leq\\min\(R,d\)\}\\;\+\\;\\underbrace\{\\tfrac\{1\}\{8\}\\,\(XW\)^\{\\odot 2\}V^\{\\top\}\}\_\{\\text\{degree\-2 lift\}\},\(18\)where⊙\\odotis the Hadamard product\. The linear term is exactly the classical rank\-RRbilinear factorization, whose rank is further capped byrank⁡\(X\)≤d\\operatorname\{rank\}\(X\)\\leq d\. The Hadamard\-square term injects degree\-2 monomial features: each column ofX​WXWlies in the≤d\\leq d\-dimensional column spacecol⁡\(X\)\\operatorname\{col\}\(X\), so the columns of\(X​W\)⊙2\(XW\)^\{\\odot 2\}lie in its symmetric squareSym2⁡\(col⁡\(X\)\)\\operatorname\{Sym\}^\{2\}\(\\operatorname\{col\}\(X\)\), of dimension≤\(d\+12\)\\leq\\binom\{d\+1\}\{2\}\. Hence the model factorizesKKthrough the fixed feature set\{𝟏\}∪col⁡\(X\)∪Sym2⁡\(col⁡\(X\)\)\\\{\\mathbf\{1\}\\\}\\cup\\operatorname\{col\}\(X\)\\cup\\operatorname\{Sym\}^\{2\}\(\\operatorname\{col\}\(X\)\), and

rank⁡\(K^\)≤min⁡\(R,\(d\+22\)\)\.\\operatorname\{rank\}\(\\hat\{K\}\)\\;\\leq\\;\\min\\\!\\Big\(R,\\ \\tbinom\{d\+2\}\{2\}\\Big\)\.\(19\)The quadratic model therefore*cannot exceed*rank\(d\+22\)\\binom\{d\+2\}\{2\}\(=21=21ford=5d=5\) no matter how large the search rankRRis\.

##### Consequence for the comparison\.

Equation equation[19](https://arxiv.org/html/2607.13246#A11.E19)predicts a sharp regime split, which Fig\.[12](https://arxiv.org/html/2607.13246#A11.F12)confirms\. WhenR≤\(d\+22\)R\\leq\\binom\{d\+2\}\{2\}\(hereR=10R=10\), the search rank is the binding constraint and the truncation is invisible: exact softplus and its surrogate behave almost identically \(panels \(a\) vs \(b\)\)\. WhenR\>\(d\+22\)R\>\\binom\{d\+2\}\{2\}\(hereR=100R=100\), the degree\-2 cap binds for the surrogate but not for exact softplus, so the exact activation attains a strictly lower error \(panel \(c\) below \(d\)\)\. In both regimes the activation is a mild perturbation of the same nonlinear\-factorization family, while the search rankRRgoverns how much of the kernel’s polynomial spectrum the model can represent, and, empirically, is what separates the optimizers\.

![Refer to caption](https://arxiv.org/html/2607.13246v1/x11.png)Figure 12:Learning\-rate stability across activation and capacity\.Final MSE \(log\) vs\. learning rate \(log\) for four optimizers fitting the Gaussian kernel asK^=ϕ​\(X​W\)​V⊤\\hat\{K\}=\\phi\(XW\)V^\{\\top\}; mean over three seeds, bands show min–max, divergent runs capped at10510^\{5\}\. The activationϕ\\phiis exact softplus or its second\-order Taylor surrogate, andRRis the search rank:\(a\)softplus,R=10R\{=\}10;\(b\)Taylor,R=10R\{=\}10;\(c\)softplus,R=100R\{=\}100;\(d\)Taylor,R=100R\{=\}100\. AtR=10R\{=\}10the rank bottleneck binds and the two activations match \(a\)≈\\approx\(b\); atR=100R\{=\}100the degree\-2 surrogate is rank\-capped \(Eq\.[19](https://arxiv.org/html/2607.13246#A11.E19)\) while exact softplus is not, so \(c\) reaches a lower error than \(d\)\. Muon attains the lowest loss over the widest stable range; SignGD diverges at large learning rates\.

## Appendix LAnalytical Solutions of the Matrix Factorization for AdamW and Muon

### L\.1Problem Setup

We consider the Matrix Factorization objective function:

ℒ​\(U,V\)=12​‖U​VT−R‖F2,\\mathcal\{L\}\(U,V\)=\\frac\{1\}\{2\}\\\|UV^\{T\}\-R\\\|\_\{F\}^\{2\},\(20\)whereU,V∈ℝN×rU,V\\in\\mathbb\{R\}^\{N\\times r\}\. At iterationtt, we seek the optimal update stepΔ​U\\Delta Ufor the factorUU\(the derivation forVVis symmetric\)\. LetGt=∇Uℒ=\(Ut​VtT−R\)​VtG\_\{t\}=\\nabla\_\{U\}\\mathcal\{L\}=\(U\_\{t\}V\_\{t\}^\{T\}\-R\)V\_\{t\}be the gradient at steptt\.

Since the global optimization is intractable, iterative optimizers solve alocal proxy problemat each step\. We derive the analytical solution for these proxy problems\. Throughout,vtv\_\{t\}\(lowercase\) denotes the per\-coordinate second\-moment estimate, with entriesvt,i​jv\_\{t,ij\}; this is distinct from the factorVtV\_\{t\}\.

### L\.2The AdamW Analytical Solution

###### Proposition L\.1\.

Under the decoupled–weight\-decay approximation of Loshchilov & Hutter \(2017\), the AdamW update is the analytical solution to minimizing the linearized loss subject to an adaptive Mahalanobis\-distance trust region, with weight decay applied as a separate additive step\.

##### The Optimization Problem\.

We first consider the linearized loss penalized by the local curvature estimated by the diagonal second\-moment estimatevtv\_\{t\}:

Δ​U†=arg​minΔ​U∈ℝN×r⁡\(⟨Gt,Δ​U⟩⏟Linear Descent\+12​η​∑i,jvt,i​j​\(Δ​Ui​j\)2⏟Adaptive Trust Region\)\.\\Delta U^\{\\dagger\}=\\operatorname\*\{arg\\,min\}\_\{\\Delta U\\in\\mathbb\{R\}^\{N\\times r\}\}\\left\(\\underbrace\{\\langle G\_\{t\},\\Delta U\\rangle\}\_\{\\text\{Linear Descent\}\}\+\\underbrace\{\\frac\{1\}\{2\\eta\}\\sum\_\{i,j\}\\sqrt\{v\_\{t,ij\}\}\\,\(\\Delta U\_\{ij\}\)^\{2\}\}\_\{\\text\{Adaptive Trust Region\}\}\\right\)\.\(21\)

##### Derivation of the trust\-region step\.

The objective is separable across coordinates\. Taking the derivative with respect to a single entryΔ​Ui​j\\Delta U\_\{ij\}and setting it to zero,

Gt,i​j\+1η​vt,i​j​Δ​Ui​j=0⟹Δ​Ui​j†=−η​Gt,i​jvt,i​j\.G\_\{t,ij\}\+\\frac\{1\}\{\\eta\}\\sqrt\{v\_\{t,ij\}\}\\,\\Delta U\_\{ij\}=0\\quad\\Longrightarrow\\quad\\Delta U^\{\\dagger\}\_\{ij\}=\-\\,\\eta\\,\\frac\{G\_\{t,ij\}\}\{\\sqrt\{v\_\{t,ij\}\}\}\.\(22\)

##### Decoupled weight decay\.

AdamW then applies weight decay as a*separate*step that is not passed through the preconditioner:

ΔUi​j∗=−η\(Gt,i​jvt,i​j\+λUt,i​j\)\.\\boxed\{\\;\\Delta U^\{\*\}\_\{ij\}=\-\\,\\eta\\left\(\\frac\{G\_\{t,ij\}\}\{\\sqrt\{v\_\{t,ij\}\}\}\+\\lambda\\,U\_\{t,ij\}\\right\)\.\\;\}\(23\)

##### Remark \(why decoupling matters\)\.

Had we instead folded the penaltyλ2​‖Ut\+Δ​U‖F2\\tfrac\{\\lambda\}\{2\}\\\|U\_\{t\}\+\\Delta U\\\|\_\{F\}^\{2\}directly into the proxy, the exact stationary point would be

\(1η​vt,i​j\+λ\)​Δ​Ui​j=−\(Gt,i​j\+λ​Ut,i​j\)⟹Δ​Ui​j=−η​Gt,i​j\+λ​Ut,i​jvt,i​j\+η​λ\.\\left\(\\frac\{1\}\{\\eta\}\\sqrt\{v\_\{t,ij\}\}\+\\lambda\\right\)\\Delta U\_\{ij\}=\-\\big\(G\_\{t,ij\}\+\\lambda U\_\{t,ij\}\\big\)\\;\\Longrightarrow\\;\\Delta U\_\{ij\}=\-\\,\\eta\\,\\frac\{G\_\{t,ij\}\+\\lambda U\_\{t,ij\}\}\{\\sqrt\{v\_\{t,ij\}\}\+\\eta\\lambda\}\.\(24\)This*couples*the decay to the preconditioner \(dividing the decay term byvt,i​j\\sqrt\{v\_\{t,ij\}\}\) and is precisely Adam withL2L\_\{2\}regularization,*not*AdamW\. AdamW is recovered by \(i\) neglectingη​λ\\eta\\lambdain the denominator and \(ii\) decoupling the decay so thatλ​Ut,i​j\\lambda U\_\{t,ij\}is*not*rescaled by the curvature\. Hence the boxed expression is an approximation, not the exact minimizer of the coupled proxy\.

##### Conclusion\.

AdamW rescales every coordinate independently based on the diagonal curvaturevt,i​j\\sqrt\{v\_\{t,ij\}\}, and adds an unpreconditioned weight\-decay term\.

### L\.3The Muon Analytical Solution

###### Proposition L\.2\.

The Muon update is the exact analytical solution to the Orthogonal Procrustes Problem\. It finds the orthonormal\-column update direction that aligns maximally with the*negative*gradient; this direction is the polar factor of−Gt\-G\_\{t\}\.

##### The Optimization Problem\.

We seek the matrixO∈ℝN×rO\\in\\mathbb\{R\}^\{N\\times r\}closest to the negative gradient−Gt\-G\_\{t\}, constrained to have orthonormal columns:

O∗=arg​minO:OT​O=Ir⁡‖\(−Gt\)−O‖F2\.O^\{\*\}=\\operatorname\*\{arg\\,min\}\_\{O\\,:\\,O^\{T\}O=I\_\{r\}\}\\;\\\|\(\-G\_\{t\}\)\-O\\\|\_\{F\}^\{2\}\.\(25\)The update is thenΔ​U=η​O∗\\Delta U=\\eta\\,O^\{\*\}\.

#### L\.3\.1Step\-by\-Step Derivation

##### Step 1: Expand the objective\.

‖−Gt−O‖F2\\displaystyle\\\|\-G\_\{t\}\-O\\\|\_\{F\}^\{2\}=Tr⁡\(\(Gt\+O\)T​\(Gt\+O\)\)\\displaystyle=\\operatorname\{Tr\}\\\!\\left\(\(G\_\{t\}\+O\)^\{T\}\(G\_\{t\}\+O\)\\right\)\(26\)=Tr⁡\(GtT​Gt\)\+Tr⁡\(OT​O\)\+2​Tr⁡\(GtT​O\)\.\\displaystyle=\\operatorname\{Tr\}\(G\_\{t\}^\{T\}G\_\{t\}\)\+\\operatorname\{Tr\}\(O^\{T\}O\)\+2\\,\\operatorname\{Tr\}\(G\_\{t\}^\{T\}O\)\.\(27\)The termTr⁡\(GtT​Gt\)\\operatorname\{Tr\}\(G\_\{t\}^\{T\}G\_\{t\}\)is constant, andTr⁡\(OT​O\)=Tr⁡\(Ir\)=r\\operatorname\{Tr\}\(O^\{T\}O\)=\\operatorname\{Tr\}\(I\_\{r\}\)=ris constant\. The cross term carries aplussign, so minimizing the norm is equivalent tominimizingTr⁡\(GtT​O\)\\operatorname\{Tr\}\(G\_\{t\}^\{T\}O\), equivalently maximizing the alignment with the negative gradient,Tr⁡\(\(−Gt\)T​O\)\\operatorname\{Tr\}\\\!\\big\(\(\-G\_\{t\}\)^\{T\}O\\big\):

minO:OT​O=Ir⁡Tr⁡\(GtT​O\)\.\\min\_\{O\\,:\\,O^\{T\}O=I\_\{r\}\}\\;\\operatorname\{Tr\}\(G\_\{t\}^\{T\}O\)\.\(28\)

##### Step 2: Singular Value Decomposition\.

Let the \(thin\) SVD of the gradient beGt=P​Σ​QTG\_\{t\}=P\\Sigma Q^\{T\}, whereP∈ℝN×rP\\in\\mathbb\{R\}^\{N\\times r\}has orthonormal columns \(PT​P=IrP^\{T\}P=I\_\{r\}\),Σ∈ℝr×r\\Sigma\\in\\mathbb\{R\}^\{r\\times r\}is diagonal with singular valuesσi≥0\\sigma\_\{i\}\\geq 0, andQ∈ℝr×rQ\\in\\mathbb\{R\}^\{r\\times r\}is orthogonal\. Then

Tr⁡\(GtT​O\)=Tr⁡\(Q​Σ​PT​O\)=Tr⁡\(Σ​\(PT​O​Q\)\),\\operatorname\{Tr\}\(G\_\{t\}^\{T\}O\)=\\operatorname\{Tr\}\\\!\\big\(Q\\Sigma P^\{T\}O\\big\)=\\operatorname\{Tr\}\\\!\\big\(\\Sigma\\,\(P^\{T\}OQ\)\\big\),\(29\)using the cyclic property of the trace\.

##### Step 3: Bound the matrixZZ\.

DefineZ=PT​O​Q∈ℝr×rZ=P^\{T\}OQ\\in\\mathbb\{R\}^\{r\\times r\}\. Each ofPP,OOhas orthonormal columns andQQis orthogonal, so each is a \(semi\-\)isometry with spectral norm11\. Henceσmax​\(Z\)≤‖PT‖2​‖O‖2​‖Q‖2=1\\sigma\_\{\\max\}\(Z\)\\leq\\\|P^\{T\}\\\|\_\{2\}\\,\\\|O\\\|\_\{2\}\\,\\\|Q\\\|\_\{2\}=1\. \(NoteZZneed*not*be orthogonal whenN\>rN\>r, sinceO​OT≠INOO^\{T\}\\neq I\_\{N\}; we only require the norm bound\.\) In particular, every diagonal entry satisfies\|Zi​i\|≤σmax​\(Z\)≤1\|Z\_\{ii\}\|\\leq\\sigma\_\{\\max\}\(Z\)\\leq 1\. The objective becomes

Tr⁡\(Σ​Z\)=∑i=1rσi​Zi​i\.\\operatorname\{Tr\}\(\\Sigma Z\)=\\sum\_\{i=1\}^\{r\}\\sigma\_\{i\}Z\_\{ii\}\.\(30\)

##### Step 4: Minimize\.

Sinceσi≥0\\sigma\_\{i\}\\geq 0and\|Zi​i\|≤1\|Z\_\{ii\}\|\\leq 1, the sum∑iσi​Zi​i\\sum\_\{i\}\\sigma\_\{i\}Z\_\{ii\}isminimizedwhenZi​i=−1Z\_\{ii\}=\-1for allii, i\.e\.Z=−IrZ=\-I\_\{r\}\(attained by the admissible choiceO=−P​QTO=\-PQ^\{T\}\)\. Solving back,

PT​O​Q\\displaystyle P^\{T\}OQ=−Ir\\displaystyle=\-I\_\{r\}\(31\)O∗\\displaystyle O^\{\*\}=−P​QT\.\\displaystyle=\-\\,P\\,Q^\{T\}\.\(32\)O∗=−PQT,ΔU=ηO∗=−ηPQT\.\\boxed\{\\;O^\{\*\}=\-\\,PQ^\{T\},\\qquad\\Delta U=\\eta\\,O^\{\*\}=\-\\,\\eta\\,PQ^\{T\}\.\\;\}\(33\)Equivalently,O∗O^\{\*\}is the polar factor of the*negative*gradient−Gt\-G\_\{t\}, and the updateΔ​U=−η​P​QT\\Delta U=\-\\eta\\,PQ^\{T\}is a genuine descent step \(it removes the singular\-value magnitudes ofGtG\_\{t\}, keeping only itsP​QTPQ^\{T\}“direction”\)\.

#### L\.3\.2Connection to Newton–Schulz

The matrixP​QTPQ^\{T\}is thepolar factorofGtG\_\{t\}\. Computing the SVD at every step is expensive, so Muon computesP​QTPQ^\{T\}via the Newton–Schulz iteration

Xk\+1=12​Xk​\(3​I−XkT​Xk\),X0=Gt‖Gt‖2,X\_\{k\+1\}=\\frac\{1\}\{2\}X\_\{k\}\\big\(3I\-X\_\{k\}^\{T\}X\_\{k\}\\big\),\\qquad X\_\{0\}=\\frac\{G\_\{t\}\}\{\\\|G\_\{t\}\\\|\_\{2\}\},\(34\)which converges quadratically toP​QTPQ^\{T\}\. The descent update is thenΔ​U=−η​X∞=−η​P​QT\\Delta U=\-\\eta\\,X\_\{\\infty\}=\-\\eta\\,PQ^\{T\}, consistent with the boxed solution above\.

### L\.4Implications for Matrix Factorization

This derivation highlights the fundamental difference in update scaling\. Both updates act on matrices inℝN×r\\mathbb\{R\}^\{N\\times r\}\.

##### 1\. AdamW scaling \(dimension dependent\)\.

Assuming active gradients \(Gi​j/vi​j≈±1G\_\{ij\}/\\sqrt\{v\_\{ij\}\}\\approx\\pm 1\) and negligible weight decay, the squared norm of the update scales with the number of entriesN×rN\\times r:

‖Δ​WAdam‖F2≈∑i,jη2=η2​\(N⋅r\)\.\\\|\\Delta W\_\{\\text\{Adam\}\}\\\|\_\{F\}^\{2\}\\approx\\sum\_\{i,j\}\\eta^\{2\}=\\eta^\{2\}\(N\\cdot r\)\.\(35\)

##### 2\. Muon scaling \(rank dependent\)\.

SinceO∗=−P​QTO^\{\*\}=\-PQ^\{T\}has orthonormal columns,\(O∗\)T​O∗=Q​PT​P​QT=Ir\(O^\{\*\}\)^\{T\}O^\{\*\}=QP^\{T\}PQ^\{T\}=I\_\{r\}, so the sign is immaterial to the norm:

‖Δ​WMuon‖F2=η2​Tr⁡\(\(O∗\)T​O∗\)=η2​Tr⁡\(Ir\)=η2​r\.\\\|\\Delta W\_\{\\text\{Muon\}\}\\\|\_\{F\}^\{2\}=\\eta^\{2\}\\,\\operatorname\{Tr\}\\\!\\big\(\(O^\{\*\}\)^\{T\}O^\{\*\}\\big\)=\\eta^\{2\}\\,\\operatorname\{Tr\}\(I\_\{r\}\)=\\eta^\{2\}r\.\(36\)

##### Conclusion\.

‖Δ​WAdam‖F‖Δ​WMuon‖F≈N\.\\frac\{\\\|\\Delta W\_\{\\text\{Adam\}\}\\\|\_\{F\}\}\{\\\|\\Delta W\_\{\\text\{Muon\}\}\\\|\_\{F\}\}\\approx\\sqrt\{N\}\.\(37\)In our experiments withN=1000N=1000, AdamW naturally takes steps≈31\.6\\approx 31\.6times larger than Muon for the same learning rateη\\eta\. This explains the necessity of scalingηMuon≈N⋅ηAdamW\\eta\_\{\\text\{Muon\}\}\\approx\\sqrt\{N\}\\,\\cdot\\,\\eta\_\{\\text\{AdamW\}\}to achieve comparable convergence rates\.