Paper deep dive
The Newton-Muon Optimizer
Zhehang Du, Weijie Su
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 97%
Last extracted: 4/2/2026, 11:51:37 PM
Summary
The paper introduces the 'Newton-Muon' optimizer, a novel optimization method for training large language models. By deriving a triplet quadratic surrogate model that incorporates gradient, output-space curvature, and input data matrices, the authors show that standard Muon can be interpreted as an implicit Newton-type method. Newton-Muon improves upon Muon by incorporating right-preconditioning based on the inverse second moment of input data, resulting in faster convergence and reduced wall-clock training time.
Entities (6)
Relation Signals (3)
Zhehang Du → authored → Newton-Muon
confidence 100% · The Newton–Muon Optimizer Zhehang Du Weijie Su
Newton-Muon → improvesupon → Muon
confidence 95% · Newton-Muon reaches the target validation loss in 6% fewer iteration steps and reduces wall-clock training time by about 4%.
Newton-Muon → utilizes → Triplet Quadratic Surrogate Model
confidence 95% · In this paper, we introduce a surrogate model... leads to a new optimizer... Newton-Muon
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The Muon optimizer has received considerable attention for its strong performance in training large language models, yet the design principle behind its matrix-gradient orthogonalization remains largely elusive. In this paper, we introduce a surrogate model that not only sheds new light on the design of Muon, but more importantly leads to a new optimizer. In the same spirit as the derivation of Newton's method, the surrogate approximates the loss as a quadratic function of the perturbation to a weight matrix $W$ using only three matrices: the gradient $G$, an output-space curvature matrix $H$, and the data matrix $Z$ that stacks the layer inputs. By minimizing this surrogate in one step and adopting a certain isotropic assumption on the weights, we obtain the closed-form update rule (up to momentum and weight decay) $W \leftarrow W - \eta \cdot \mathrm{msgn}(G(ZZ^\top)^{-1})$, where $\eta$ is the learning rate and $\mathrm{msgn}(X)=UV^\top$ if $X=USV^\top$ is a compact singular value decomposition. This new optimization method, which we refer to as Newton-Muon, shows that standard Muon can be interpreted as an implicit Newton-type method that neglects the right preconditioning induced by the input second moment. Empirically, on a reproduction of the earliest publicly released Modded-NanoGPT speedrun configuration using Muon for GPT-2 pretraining, Newton-Muon reaches the target validation loss in 6\% fewer iteration steps and reduces wall-clock training time by about 4\%.
Tags
Links
- Source: https://arxiv.org/abs/2604.01472v1
- Canonical: https://arxiv.org/abs/2604.01472v1
Trouble viewing inline? Open PDF directly →
Full Text
123,931 characters extracted from source content.
Expand or collapse full text
The Newton–Muon Optimizer Zhehang Du Weijie Su (University of Pennsylvania ) Abstract The Muon optimizer has received considerable attention for its strong performance in training large language models, yet the design principle behind its matrix-gradient orthogonalization remains largely elusive. In this paper, we introduce a surrogate model that not only sheds new light on the design of Muon, but more importantly leads to a new optimizer. In the same spirit as the derivation of Newton’s method, the surrogate approximates the loss as a quadratic function of the perturbation to a weight matrix W using only three matrices: the gradient G, an output-space curvature matrix H, and the data matrix Z that stacks the layer inputs. By minimizing this surrogate in one step and adopting a certain isotropic assumption on the weights, we obtain the closed-form update rule (up to momentum and weight decay) W←W−η⋅msgn(G(ZZ⊤)−1),W← W-η·msgn(G(Z )^-1), where η is the learning rate and msgn(X)=UV⊤msgn(X)=UV if X=USV⊤X=USV is a compact singular value decomposition. This new optimization method, which we refer to as Newton–Muon, shows that standard Muon can be interpreted as an implicit Newton-type method that neglects the right preconditioning induced by the input second moment. Empirically, on a reproduction of the earliest publicly released Modded-NanoGPT speedrun configuration using Muon for GPT-2 pretraining, Newton–Muon reaches the target validation loss in 6% fewer iteration steps and reduces wall-clock training time by about 4%. †footnotetext: Emails: duz,suw@wharton.upenn.edu. †footnotetext: Code is available at https://github.com/zhehangdu/Newton-Muon. 1 Introduction Since early 2025, there has been a surge of interest in matrix-structured optimization methods for training deep neural networks and large language models (LLMs). A prominent optimizer at the center of this flurry of research activity is Muon (Jordan et al., 2024b), closely related to spectral gradient descent (Carlson et al., 2015). Specifically, letting f(W)f(W) be the loss function and WW be a layer weight matrix, we consider the optimization problem minW∈ℝm×nf(W). _W ^m× nf(W). Let G≔∇Wf(W)∈ℝm×nG _Wf(W) ^m× n be the gradient matrix or, in practice, an approximation computed from a mini-batch. Writing the matrix sign msgn(G)≔UV⊤msgn(G) UV if G=USV⊤G=USV is a compact singular value decomposition (SVD), Muon updates the weight matrix as111In practice, Muon approximates msgn(G)msgn(G) using Newton–Schulz iterations and momentum is applied to the gradient. W←W−η⋅msgn(G)W← W-η·msgn(G) for some learning rate η. Compared to AdamW (Kingma and Ba, 2014; Loshchilov and Hutter, 2017), Muon has been reported to offer a faster convergence rate and lower wall-clock time to reach the same level of loss across a broad range of model sizes (Liu et al., 2025; Shah et al., 2025; Wen et al., 2025). Furthermore, it has been used to train state-of-the-art open-source models (Team et al., 2026; Zeng et al., 2026), and many extensions of Muon have been proposed (Li et al., 2025; Pethick et al., 2025; Ahn et al., 2025; He et al., 2025; Xu et al., 2026; Qi et al., 2026; Gu and Xie, 2026). While matrix-based optimizers have demonstrated highly effective empirical performance, the theoretical mechanisms underlying Muon remain largely mysterious. For instance, it is natural to ask why preserving the matrix structure of the gradient is beneficial and why simply discarding its singular values is empirically effective. This is, however, not entirely surprising, as the nonconvex nature of deep learning optimization makes it notoriously difficult to analyze. Given that a rigorous theoretical understanding remains out of reach, a practical approach is to establish an intuitive yet principled framework for designing deep learning optimizers (Bernstein and Newhouse, 2024; Pethick et al., 2025; Lau et al., 2025; Gong et al., 2026). Among these, Su (2025) introduced the isotropic curvature model as a surrogate for approximating the loss function by assuming isotropic curvature for preconditioning and isotropic input activations. A one-step descent analysis applied to this model suggests that the optimal update direction naturally preserves the matrix subspace of the gradient, thereby partially justifying gradient orthogonalization. However, the unspecified curvature function in the isotropic curvature model limits its utility for deriving practical optimization methods for LLM training. In this paper, we address this challenge by introducing a more tractable surrogate model for understanding these optimizers and, more importantly, for proposing a new method. Taking Q∈ℝm×nQ ^m× n to be a potential update direction and letting H denote a curvature matrix,222Here, the curvature matrix H is not the full parameter-space Hessian with respect to the vectorized weights. In particular, the curvature matrix H has size m×m× m. this paper introduces a surrogate model of the form f(W−Q)−f(W)≈−tr(QG⊤)+12Ntr(HQ(ZZ⊤)Q⊤),f(W-Q)-f(W)≈-tr(QG )+ 12Ntr (HQ(ZZ )Q ), (1) where Z=[1,…,N]Z=[ z_1,…, z_N] denotes the collection of all activation inputs from the N training data points to the layer W. Because it involves three components—namely G, H, and Z—this model is referred to as the triplet quadratic surrogate model. By design, the gradient matrix G and the update direction Q naturally retain their matrix forms. The linear term −tr(QG⊤)-tr(QG ) captures the first-order approximation, as in the isotropic curvature model, while the quadratic term 12Ntr(HQ(ZZ⊤)Q⊤) 12Ntr (HQ(ZZ )Q ) is considerably simpler than the curvature function in the isotropic curvature model. Notably, this quadratic term is the matrix-form expression induced by the Kronecker-factored curvature approximation in K-FAC (Martens and Grosse, 2015). This triplet quadratic surrogate for approximating f(W−Q)−f(W)f(W-Q)-f(W) is essentially of the form “−linear gradient term+quadratic curvature term-linear gradient term+quadratic curvature term”, which closely resembles the derivation of Newton’s method. As such, this surrogate model serves as a means of deriving a Newton-type method in a manner that fully leverages the matrix structure of the gradient via a one-step descent analysis; that is, by minimizing the right-hand side of (1) over Q. In fact, although second-order methods such as Newton’s method are rarely used in training large-scale neural networks, recent work demonstrates a 33–5×5× iteration speedup when using the Gauss–Newton method compared to common optimizers in deep learning (Abreu et al., 2025). Thus, there is strong motivation to develop an implicit Newton-type method for LLM training, provided that the per-iteration computational cost remains roughly comparable to that of AdamW or Muon. Despite the ease of minimizing the triplet quadratic surrogate model, a key difficulty arises because the optimal solution involves the unknown curvature matrix H. To circumvent this, we make an assumption on the displacement W−W⋆W-W , which is the difference between the current weight matrix W and an optimal weight matrix W⋆W .333Since there are generally exponentially many optimal weight matrices, we simply consider one that is locally closest to W. We assume that this displacement is isotropic, in the sense that no direction is favored over another on average. Surprisingly, under this least-informative assumption, we obtain a closed-form expression for the Newton-type update: Q⋆∝msgn(G(ZZ⊤)−1)Q (G(Z )^-1). Accordingly, our new method updates the weight matrix according to the rule444In practice, momentum and other common implementation tricks are also used. W←W−η⋅msgn(G(ZZ⊤)−1).W← W-η·msgn(G(Z )^-1). This method differs from standard Muon in that the gradient matrix G is right-preconditioned by the inverse second moment of the input data. This can be interpreted as incorporating the geometry of the data distribution to yield an update direction distinct from that of Muon. To implement Newton–Muon efficiently in practice, we maintain a running estimate of the second moment matrix ZZ⊤Z and recompute both this estimate and its inverse only periodically, rather than at every optimization step. Between refreshes, the cached inverse is reused. We compute the damped inverse (ZZ⊤+γIn)−1(Z +γ I_n)^-1 via a Cholesky factorization followed by triangular solves. Even for very large matrices, this computational cost is no more than about 10×10× that of a single matrix multiplication on a modern GPU. This step is then followed by a Newton–Schulz polynomial approximation to the matrix sign, such as the Polar Express method (Amsel et al., 2025) or Gram Newton–Schulz (Zhang et al., 2026). Because this optimization method is derived in the same spirit as Newton’s method, we call it Newton–Muon. Standard Muon can be recovered as a special case when the input data’s second moment is approximately isotropic, that is, ZZ⊤∝InZZ I_n. Put differently, Muon can be viewed as an implicit Newton-type method that does not account for the geometry of the input data. However, this simplification is not consistent with our observation that ZZ⊤Z is highly anisotropic in practice. The right-preconditioning by the input data distribution in the update rule endows Newton–Muon with its empirical advantage over standard Muon. Indeed, compared to our reproduction of the earliest publicly logged Modded-NanoGPT speedrun (Jordan et al., 2024a) configuration using Muon, Newton–Muon achieves a notable 6%6\% reduction in the number of iterations required to reach the target validation loss, as shown in Figure 1. In terms of wall-clock time, replacing standard Muon with Newton–Muon yields a reduction of about 4%4\%. Figure 1: Top: short track Record #4 validation loss comparison on the Modded-NanoGPT speedrun benchmark. Record #4 is the earliest publicly released configuration using Muon, and our reproduction on a single H100 GPU is denoted Muon. Newton–Muon adds the activation right-preconditioner before the Newton–Schulz iterations. Newton–Muon reaches the Muon baseline final validation loss in 6%6\% fewer steps; despite a 1.8%1.8\% higher per-step cost from right-preconditioning, it reduces wall-clock time to that loss by about 4%4\%. Bottom: CIFAR-10 experiments (Appendix C) on a 32-layer residual MLP show that Newton–Muon outperforms both Muon and AdamW in both per-step efficiency and overall wall-clock time. 1.1 Structure of the Paper Section 2 introduces the triplet quadratic surrogate model and derives the Newton–Muon update. Section 3 analyzes Newton–Muon on a simple quadratic case study. Section 4 develops a one-step analysis with numerical experiments under spiked activation. Section 5 presents our experiments on LLMs. Section 6 discusses limitations and open directions. Appendix A provides the LLM training configurations. Appendix B details efficient computation of the activation second moment inverse. Appendix C provides the CIFAR-10 training configurations. Appendix D derives the Kronecker-factored curvature. Appendix E gives quadratic score formulas under isotropic activation. Appendix F discusses the non-isotropic assumption. 1.2 Related Work Matrix-based optimizers. Existing matrix-based optimizers can be grouped into several broad directions. One line uses matrix structure for preconditioning, curvature approximation, or adaptive updates, including spectral descent (Carlson et al., 2015), K-FAC (Martens and Grosse, 2015; George et al., 2018), Shampoo (Gupta et al., 2018; Morwani et al., 2024), and SOAP (Vyas et al., 2025), which performs Adafactor-style (Shazeer and Stern, 2018) updates. Another line improves training efficiency through approximations, for example by combining preconditioning with variance reduction, as in MARS (Yuan et al., 2024), or by using low-rank gradient projections, as in GaLore (Zhao et al., 2024; Su et al., 2025), or by accelerating matrix-function evaluation, as in PRISM (Yang et al., 2026). More recent work studies matrix geometry itself as the design principle, including LMO-based methods such as Scion (Pethick et al., 2025), cheaper Muon variants such as Dion2 (Ahn et al., 2025), adaptive choices of learning rate such as PolarGrad (Lau et al., 2025), and rotated coordinate updates such as ARO (Gong et al., 2026). Understanding Muon. A common view is that Muon implements normalized steepest descent under the spectral norm (Bernstein and Newhouse, 2024; Crawshaw et al., 2025; Riabinin et al., 2025), which is further generalized by Lion-K (Chen et al., 2025). We instead connect Muon to the classical Euclidean Newton step under a local quadratic surrogate, providing a complementary explanation for its strong performance. Related theory studies Muon’s implicit bias and geometric interpretation (Fan et al., 2025; Pethick et al., 2025), including simplicity bias (Dragutinović and Ranganath, 2026), separates curvature and gradient anisotropy (Lau et al., 2025), and analyzes convergence in stochastic nonconvex settings and structured models (Li and Hong, 2025; Sato et al., 2025; Ma et al., 2026), including under heavy-tailed gradient noise (Yu et al., 2026). Beyond asymptotic results, local one-step analysis has also been studied (Su, 2025; Davis and Drusvyatskiy, 2025; Gonon et al., 2026). Kronecker-factored curvature. Layerwise curvature factorizations and matrix preconditioners have been widely studied, notably in K-FAC (Martens and Grosse, 2015) and Shampoo (Gupta et al., 2018). We start from the same structural observation: for a given layer, the parameter-space curvature of a local quadratic model can be approximated by the Kronecker product (ZZ⊤/N)⊗H(ZZ /N) H, where ZZ⊤/NZZ /N is the activation second moment and HH is an output-space curvature factor. K-FAC explicitly estimates HH as a generalized Gauss–Newton/Fisher factor from gradient second moments and then applies the two-sided preconditioner H−1G(ZZ⊤)−1H^-1G(ZZ )^-1. In contrast, we show that Muon’s matrix sign can be interpreted as an implicit left preconditioner that approximates the effect of H−1H^-1. 2 Derivation of Newton–Muon In this section, we first propose an amenable surrogate function for approximating the objective function we wish to minimize, followed by the derivation of the Newton–Muon method via a one-step minimization of this surrogate. We then draw a connection between Newton’s method and standard Muon to provide a principled interpretation of the latter. 2.1 Triplet Quadratic Surrogate Starting from the current iterate, we estimate the loss change induced by a candidate update direction using a local quadratic model, in the same spirit as the derivation of Newton’s method. Assume the loss function f:ℝm×n→ℝf:R^m× n is twice continuously differentiable, and let G≔∇Wf(W)∈ℝm×nG _Wf(W) ^m× n be the gradient matrix at the current iterate W∈ℝm×nW ^m× n. For each sample i, let i∈ℝn z_i ^n denote the layer input, which is the activated output from the previous layer, and write Z=[1,…,N]∈ℝn×NZ=[ z_1,…, z_N] ^n× N. For a candidate update direction Q∈ℝm×nQ ^m× n, the layer output for sample i changes by −Qi∈ℝm-Q z_i ^m. To motivate the quadratic term, we write the averaged loss as f(W)=(1/N)∑i=1NLi(Wi)f(W)=(1/N) _i=1^NL_i(W z_i). Under the perturbation W↦W−QW W-Q, the i-th output changes from WiW z_i to Wi−QiW z_i-Q z_i. By the integral remainder form of Taylor’s theorem, Li(Wi−Qi)=Li(Wi)−∇Li(Wi)⊤(Qi)+h(Qi,Wi),L_i(W z_i-Q z_i)=L_i(W z_i)-∇ L_i(W z_i) (Q z_i)+h(Q z_i,W z_i), where h(Qi,Wi)≔(Qi)⊤[∫01(1−t)∇2Li(Wi−tQi)dt](Qi).h(Q z_i,W z_i) (Q z_i) [ _0^1(1-t)∇^2L_i(W z_i-tQ z_i)dt ](Q z_i). We then approximate this path-dependent integral ∫01(1−t)∇2Li(Wi−tQi)dt≈(1/2)H _0^1(1-t)∇^2L_i(W z_i-tQ z_i)dt≈(1/2)H by a fixed average curvature matrix H shared across samples, so that the contribution of sample i to the quadratic term is approximated by 1/(2N)(Qi)⊤H(Qi)1/(2N)(Q z_i) H(Q z_i). Summing over the N samples gives the surrogate objective minQ∈ℝm×nJ(Q)=−tr(QG⊤)+12N∑i=1N(Qi)⊤H(Qi)≔−tr(QG⊤)+12Ntr(HQZZ⊤Q⊤). _Q ^m× nJ(Q)=-tr(QG )+ 12N _i=1^N(Q z_i) H(Q z_i) -tr(QG )+ 12Ntr (HQZZ Q ). (2) This objective is closely related to the isotropic curvature model (Su, 2025): minQ∈ℝm×n−tr(QG⊤)+ζh(‖Qζ‖), _Q ^m× n-tr(QG )+E_ζh(\|Qζ\|), where ζ is sampled uniformly from the unit sphere and h is a univariate curvature function. The linear term is the same as in (2), but the higher-order term is modeled isotropically through a radial curvature function h, assumed to have super-quadratic growth. In contrast, our triplet surrogate replaces the isotropic sampling of ζ with the empirical distribution from the columns of Z, removes the isotropy assumption on the curvature, and explicitly captures the interaction among H, Z, and G in (2). 2.2 Minimization of the Triplet Model Minimizing (2) over Q is straightforward, but it becomes practically useful only if we can establish a relationship between the curvature matrix H and the gradient matrix G. To this end, we first establish a connection between (2) and the full parameter-space second-order expansion. Define the parameter-space Hessian ℋW≔∇vec(W)2f(W)∈ℝ(mn)×(mn)H_W ∇^2_vec(W)f(W) ^(mn)×(mn), where vec(W)∈ℝmnvec(W) ^mn stacks WW into a vector using a standard convention. The second-order Taylor expansion around vec(W)vec(W) gives f(W−Q)≈f(W)−tr(QG⊤)+12vec(Q)⊤ℋWvec(Q).f(W-Q)≈ f(W)-tr(QG )+ 12vec(Q) H_Wvec(Q). (3) For the quadratic proxy of (3), we employ the approximation vec(G)≈ℋWvec(W−W⋆)vec(G)≈H_Wvec(W-W ), where W⋆W denotes a nearby local minimizer, or more generally, a nearby reference point at which ∇Wf(W⋆) _Wf(W ) is assumed to be negligibly small. To relate ℋWH_W to the curvature matrix H in the surrogate model (2), we use the Kronecker-factored approximation (see its application in K-FAC (Martens and Grosse, 2015) and Appendix D for details): ℋW≈(ZZ⊤/N)⊗H,H_W≈(ZZ /N) H, (4) where ⊗ denotes the Kronecker product. We then obtain vec(G)≈((ZZ⊤/N)⊗H)vec(W−W⋆)=vec(H(W−W⋆)(ZZ⊤/N)).vec(G)≈ ((ZZ /N) H )vec(W-W )=vec (H(W-W )(ZZ /N) ). This immediately implies G≈H(W−W⋆)ZZ⊤/N.G≈ H(W-W )ZZ /N. In light of the above, we make the following assumption. Assumption 1. We assume G=H(W−W⋆)ZZ⊤/N,H≻0,ZZ⊤≻0G=H(W-W )ZZ /N, H 0, ZZ 0. The last two conditions ensure that the quadratic surrogate in (2) is strictly convex in QQ and hence has a unique minimizer. If the curvature matrix H is not positive definite but merely non-degenerate, we can instead solve for a stationary point of the surrogate model. Under this assumption, the optimal solution to the surrogate objective (2) takes a closed-form expression that does not explicitly involve the unknown curvature matrix H. Proposition 1. Denote by ΣW≔(W−W⋆)(W−W⋆)⊤∈ℝm×m _W (W-W )(W-W ) ^m× m the displacement second moment matrix and write ΣW1/2 _W^1/2 for its unique positive semidefinite square root. Then under Assumption 1, the unique minimizer Q⋆Q of (2) takes the form Q⋆=ΣW1/2msgn(ΣW1/2G(ZZ⊤)−1).Q = _W^1/2msgn ( _W^1/2G(ZZ )^-1 ). (5) Proof of Proposition 1. Define the compact SVD W−W⋆=UQSQVQ⊤W-W =U_QS_QV_Q , where UQ∈ℝm×rU_Q ^m× r and VQ∈ℝn×rV_Q ^n× r have orthonormal columns, and SQ∈ℝr×rS_Q ^r× r is a diagonal matrix with positive entries. Then ΣW=UQSQ2UQ⊤ _W=U_QS_Q^2U_Q and ΣW1/2=UQSQUQ⊤ _W^1/2=U_QS_QU_Q . From Assumption 1, we have G(ZZ⊤)−1=H(W−W⋆)/NG(ZZ )^-1=H(W-W )/N, so ΣW1/2G(ZZ⊤)−1=1NΣW1/2H(W−W⋆)=1NUQSQ(UQ⊤HUQ)SQVQ⊤. _W^1/2G(ZZ )^-1= 1N _W^1/2H(W-W )= 1NU_QS_Q (U_Q HU_Q )S_QV_Q . Since H≻0H 0 and UQU_Q has full column rank, the r×r× r matrix UQ⊤HUQU_Q HU_Q is symmetric positive definite. Hence, SQ(UQ⊤HUQ)SQS_Q (U_Q HU_Q )S_Q is also symmetric positive definite. Therefore, msgn(ΣW1/2G(ZZ⊤)−1)=UQVQ⊤.msgn ( _W^1/2G(ZZ )^-1 )=U_QV_Q . Multiplying on the left by ΣW1/2 _W^1/2 gives ΣW1/2msgn(ΣW1/2G(ZZ⊤)−1)=UQSQUQ⊤UQVQ⊤=UQSQVQ⊤=W−W⋆. _W^1/2msgn ( _W^1/2G(ZZ )^-1 )=U_QS_QU_Q U_QV_Q =U_QS_QV_Q =W-W . Under Assumption 1, the quadratic term in (2) is equivalent to the vectorized quadratic term in (3). Thus, the solution Q⋆Q to (2), that is, the Newton direction, satisfies vec(Q⋆)=ℋW−1vec(G)=ℋW−1ℋWvec(W−W⋆)=vec(W−W⋆)vec(Q )=H_W^-1vec(G)=H_W^-1H_Wvec(W-W )=vec(W-W ). Hence Q⋆=W−W⋆Q =W-W , from which (5) follows. ∎ 2.3 Newton–Muon The Newton-type update (5) can be readily used as long as the unknown displacement second moment ΣW=(W−W⋆)(W−W⋆)⊤ _W=(W-W )(W-W ) can be estimated. While the proof of Proposition 1 shows that Q⋆=W−W⋆Q =W-W , meaning that the right-hand side of (5) inherently depends on Q⋆Q , a closer look reveals that ΣW _W is a much coarser object to approximate than Q⋆Q itself. Specifically, approximating ΣW _W does not require knowledge of the right singular spaces of W−W⋆W-W . Perhaps the least-informative approximation is to assume that ΣW _W is a multiple of the identity matrix. This isotropic proxy, ΣW∝Im _W I_m, is plausible, especially at initialization and during the early stages of training when the columns of W−W⋆W-W are approximately independent and no specific directional structure has yet emerged. Since ΣW=(W−W⋆)(W−W⋆)⊤ _W=(W-W )(W-W ) aggregates the second moments of the columns of W−W⋆W-W , applying this isotropic proxy amounts to treating these unknown column directions as having no preferred orientation in aggregate. This yields a closed-form update direction depending only on observable quantities, as shown in the following result. Theorem 1. Under Assumption 1 and the isotropic proxy ΣW∝Im _W I_m, we have Q⋆∝msgn(G(ZZ⊤)−1).Q (G(ZZ )^-1 ). Although this is a corollary of Proposition 1, we choose to present it as a theorem because it formally establishes the core method introduced in this paper. At iteration t, we update the weight matrix according to the rule: Wt+1=Wt−ηt⋅msgn(Gt(ZtZt⊤)−1),W_t+1=W_t- _t·msgn (G_t(Z_tZ_t )^-1 ), where ηt _t denotes the learning rate. We call this method Newton–Muon because it implicitly acts as a Newton-type method, derived by minimizing a quadratic surrogate. Strictly speaking, when m>nm>n, it mathematically cannot hold that ΣW∝Im _W I_m, since ΣW=(W−W⋆)(W−W⋆)⊤ _W=(W-W )(W-W ) has rank at most n. Thus, the relation ΣW∝Im _W I_m should be understood purely as an isotropic proxy for the unknown displacement second moment. The following result shows that the Newton–Muon update is always a descent direction. Proposition 2. For the Newton–Muon direction Q=msgn(G(ZZ⊤)−1)Q=msgn (G(Z )^-1 ), we have tr(G⊤Q)≥0tr(G Q)≥ 0. Proof of Proposition 2. Let Gr≔G(ZZ⊤)−1G_r G(Z )^-1, and write its compact SVD as Gr=UGSGVG⊤G_r=U_GS_GV_G . Then Q=msgn(Gr)=UGVG⊤Q=msgn(G_r)=U_GV_G . Since G=GrZZ⊤G=G_rZ , we have tr(G⊤Q)=tr((GrZZ⊤)⊤Q)=tr(ZZ⊤Gr⊤Q)=tr(ZZ⊤VGSGVG⊤).tr(G Q)=tr ((G_rZ ) Q )=tr (Z G_r Q )=tr (Z V_GS_GV_G ). Now ZZ⊤Z and VGSGVG⊤V_GS_GV_G are both positive semidefinite. Therefore, the trace of their product is nonnegative, so tr(G⊤Q)≥0tr(G Q)≥ 0. ∎ Making a connection to the standard Muon optimizer, we establish the following result. Corollary 1 (Isotropic activations recover Muon). Under the assumptions of Theorem 1, if additionally ZZ⊤∝InZZ I_n, then Q⋆∝msgn(G).Q (G). That is, the triplet quadratic surrogate model recovers the standard Muon update in this case. However, as we will demonstrate in Section 5, ZZ⊤ZZ is highly anisotropic in practice. Therefore, the right preconditioner (ZZ⊤)−1(ZZ )^-1, which is readily available during training, can significantly alter the singular spaces of the gradient matrix G. Furthermore, the added computational cost of computing (G(ZZ⊤)−1)(G(Z )^-1) is insignificant compared to the empirical performance gains of Newton–Muon over standard Muon. It is also instructive to interpret Newton–Muon through the lens of orthogonal equivariance. Standard Muon is basis-free in the sense that for any orthogonal matrices Om∈ℝm×mO_m ^m× m and On∈ℝn×nO_n ^n× n, we have msgn(OmGOn)=Ommsgn(G)On.msgn(O_mGO_n)=O_mmsgn(G)O_n. Thus, rotating the gradient matrix on the left or right correspondingly rotates the Muon update in the exact same manner. For Newton–Muon, the analogous commutative diagram holds if the right rotation of the gradient matrix G is accompanied by the coordinate transformation Z↦On⊤Z O_n Z. Consequently, msgn(OmGOn((On⊤Z)(On⊤Z)⊤)−1)=Ommsgn(G(ZZ⊤)−1)On.msgn (O_mGO_n ((O_n Z)(O_n Z) )^-1 )=O_mmsgn (G(ZZ )^-1 )O_n. In contrast, if one rotates GG on the right while keeping ZZ fixed, right equivariance generally fails unless OnO_n commutes with ZZ⊤ZZ , as illustrated in Figure 2. GGmsgn(G)msgn(G)OmGOnO_mGO_nOmmsgn(G)OnO_mmsgn(G)O_n (G,Z)(G,Z)msgn(G(ZZ⊤)−1)msgn (G(Z )^-1 )(OmGOn,On⊤Z)(O_mGO_n,O_n Z)Ommsgn(G(ZZ⊤)−1)OnO_mmsgn (G(Z )^-1 )O_n Figure 2: Left: standard Muon is orthogonally equivariant. Right: Newton–Muon takes the pair (G,Z)(G,Z) as input, and the diagram commutes if the right rotation of G is accompanied by the transformation Z↦On⊤Z O_n Z. 3 Convergence Analysis of Newton–Muon: Case Study This section presents a case study of Newton–Muon in a simple quadratic model under a single spike assumption on ZZ⊤ZZ , with one spiked eigendirection and an isotropic bulk. A fully general convergence analysis is difficult, so we study this simple model in which explicit rates can be obtained. The main finding is that Newton–Muon achieves a contraction rate independent of the spike condition number κ, whereas the rates of gradient descent and Muon scale as 1/κ1/κ as κ increases. We consider the iterative optimization of WW under full-batch gradient descent, Muon, and Newton–Muon, using the learning rate sequence ηtt≥0\ _t\_t≥ 0. We consider the objective minWf(W)≔12‖WZ−W⋆Z‖F2, _Wf(W) 12\|WZ-W Z\|_F^2, where Z∈ℝn×NZ ^n× N is fixed, W⋆∈ℝm×nW ^m× n is the ground-truth matrix, and ∥⋅∥F\|·\|_F denotes the Frobenius norm. The exact gradient is ∇Wf(W)=(W−W⋆)ZZ⊤. _Wf(W)=(W-W )ZZ . (6) Denote WtW_t as the weight matrix at iteration t. The resulting updates are GD:Wt+1 : W_t+1 =Wt−ηt∇Wf(Wt)=Wt−ηt(Wt−W⋆)ZZ⊤, =W_t- _t _Wf(W_t)=W_t- _t(W_t-W )Z , (7) Muon:Wt+1 : W_t+1 =Wt−ηtmsgn(∇Wf(Wt))=Wt−ηtmsgn((Wt−W⋆)ZZ⊤), =W_t- _tmsgn ( _Wf(W_t) )=W_t- _tmsgn ((W_t-W )Z ), (8) Newton–Muon:Wt+1 --Muon: W_t+1 =Wt−ηtmsgn(∇Wf(Wt)(ZZ⊤)−1)=Wt−ηtmsgn(Wt−W⋆). =W_t- _tmsgn ( _Wf(W_t)(Z )^-1 )=W_t- _tmsgn(W_t-W ). (9) Before going into the theoretical details, we can already see the potential advantage of Newton–Muon over gradient descent and Muon. From the update equations (7)–(9), gradient descent and Muon apply the right-multiplication factor ZZ⊤Z to the displacement Wt−W⋆W_t-W . When ZZ⊤Z is highly anisotropic, as in LLM training, where the activation matrix Z has bounded stable rank (Davis and Drusvyatskiy, 2025) and a very large condition number (Table 4), this anisotropy can distort the effective update, over-emphasizing some directions toward W⋆W while suppressing others. In contrast, Newton–Muon explicitly cancels this right-side anisotropy by multiplying by (ZZ⊤)−1(Z )^-1. Single spike model. We now introduce the single-spike model and restrict to a low-dimensional invariant subspace spanned by a single spiked eigendirection of ZZ⊤ZZ and its isotropic complement. We initialize W0=0W_0=0 and assume that the singular values of ZZ satisfy σ12=κ>1 _1^2=κ>1 and σ22=⋯=σn2=1 _2^2=·s= _n^2=1, so that ZZ⊤=UZdiag(κ,1,…,1)UZ⊤,κ>1,ZZ =U_Zdiag(κ,1,…,1)U_Z , κ>1, (10) for some orthogonal matrix UZU_Z. Let 1∈ℝn e_1 ^n denote the first column of UZU_Z. We assume W⋆W has the following structured rank-r form W⋆=−1(α1,01⊤+β1,01⊤)−∑i=2ri(βi,0i⊤),W =- u_1 ( _1,0 e_1 + _1,0 b_1 )- _i=2^r u_i ( _i,0 b_i ), (11) where 1,…,r∈ℝm u_1,…, u_r ^m are orthonormal, and 1,…,r∈ℝn b_1,…, b_r ^n are orthonormal with i⟂1 b_i e_1 for all i. In particular, this requires r≤mr≤ m and r≤n−1r≤ n-1. We assume (α1,0,β1,0)≠(0,0)( _1,0, _1,0)≠(0,0) and βi,0≠0 _i,0≠ 0 for i≥2i≥ 2. This assumption ensures that (i) W⋆W is rank-r, and (i) only the first mode mixes the spiked direction 1 e_1 with a non-spike direction; all other modes lie entirely in the isotropic subspace 1⟂ e_1 . Since W0=0W_0=0, the initial residual W0−W⋆=−W⋆W_0-W =-W lies in the same family. We will show that, for all t, Wt−W⋆=1(α1,t1⊤+β1,t1⊤)+∑i=2ri(βi,ti⊤).W_t-W = u_1 ( _1,t e_1 + _1,t b_1 )+ _i=2^r u_i ( _i,t b_i ). (12) The main consequence of (10) and (11) is that the dynamics of gradient descent, Muon, and Newton–Muon decouple across the modes i=1,…,ri=1,…,r, yielding explicit one-dimensional recursions for the coefficients. We adopt the convention that msgn(0)=0msgn(0)=0, β/|β|=0β/|β|=0 when β=0β=0, and (α,β)/α2+β2=(0,0)(α,β)/ α^2+β^2=(0,0) when (α,β)=(0,0)(α,β)=(0,0). Lemma 1 (Dynamics under the single spike model). Under (10) and (11), gradient descent (7), Muon (8), and Newton–Muon (9) preserve the decomposition (12). In particular, the dynamics decouple across modes i=1,…,ri=1,…,r, and the coefficients satisfy the following recursions. • For gradient descent, α1,t+1=(1−ηtκ)α1,t,βi,t+1=(1−ηt)βi,t,i=1,…,r. _1,t+1=(1- _tκ) _1,t, _i,t+1=(1- _t) _i,t, i=1,…,r. (13) • For Muon, α1,t+1=α1,t−ηtκα1,tκ2α1,t2+β1,t2,βi,t+1=β1,t−ηt(β1,t/κ2α1,t2+β1,t2),i=1,βi,t−ηt(βi,t/|βi,t|),i=2,…,r. _1,t+1= _1,t- _t κ _1,t κ^2 _1,t^2+ _1,t^2, _i,t+1= cases _1,t- _t ( _1,t/ κ^2 _1,t^2+ _1,t^2 ),&i=1,\\ _i,t- _t ( _i,t/| _i,t| ),&i=2,…,r. cases (14) • For Newton–Muon, α1,t+1=α1,t−ηtα1,tα1,t2+β1,t2,βi,t+1=β1,t−ηt(β1,t/α1,t2+β1,t2),i=1,βi,t−ηt(βi,t/|βi,t|),i=2,…,r. _1,t+1= _1,t- _t _1,t _1,t^2+ _1,t^2, _i,t+1= cases _1,t- _t ( _1,t/ _1,t^2+ _1,t^2 ),&i=1,\\ _i,t- _t ( _i,t/| _i,t| ),&i=2,…,r. cases (15) Proof of Lemma 1. Define 1,t≔α1,t1+β1,t1 v_1,t _1,t e_1+ _1,t b_1 and, for i≥2i≥ 2, i,t≔βi,ti v_i,t _i,t b_i, so that Wt−W⋆=∑i=1rii,t⊤.W_t-W = _i=1^r u_i v_i,t . (16) Since ii=1r\ u_i\_i=1^r are orthonormal and ii=1r\ b_i\_i=1^r are orthonormal with i⟂1 b_i e_1, we have i,t⟂j,t v_i,t v_j,t for i≠ji≠ j for every t. Under (10), ZZ⊤1=κ1Z e_1=κ e_1, and the eigenspace with eigenvalue 11 is 1⟂ e_1 , so ZZ⊤i=iZZ b_i= b_i for all i=1,…,ri=1,…,r. Hence ZZ⊤1,t=κα1,t1+β1,t1∈span1,1,ZZ⊤i,t=i,t∈spani(i≥2).Z v_1,t=κ _1,t e_1+ _1,t b_1 \ e_1, b_1\, Z v_i,t= v_i,t \ b_i\\ \ (i≥ 2). (17) Moreover, the transformed vectors ZZ⊤i,ti=1r\Z v_i,t\_i=1^r remain mutually orthogonal across i. Using (6) and (16), ∇Wf(Wt)=(Wt−W⋆)ZZ⊤=∑i=1ri((ZZ⊤i,t)⊤). _Wf(W_t)=(W_t-W )Z = _i=1^r u_i ((Z v_i,t) ). (18) In particular, by (17), (ZZ⊤1,t)⊤=(κα1,t)1⊤+β1,t1⊤,(ZZ⊤i,t)⊤=βi,ti⊤(i≥2).(Z v_1,t) =(κ _1,t) e_1 + _1,t b_1 , (Z v_i,t) = _i,t b_i \ (i≥ 2). For gradient descent, substituting (18) into (7) gives Wt+1−W⋆=(Wt−W⋆)−ηt(Wt−W⋆)ZZ⊤,W_t+1-W =(W_t-W )- _t(W_t-W )Z , and matching coefficients in the basis 11⊤,11⊤,ii⊤\ u_1 e_1 , u_1 b_1 , u_i b_i \ yields α1,t+1=(1−ηtκ)α1,t _1,t+1=(1- _tκ) _1,t and βi,t+1=(1−ηt)βi,t _i,t+1=(1- _t) _i,t for all i=1,…,ri=1,…,r, i.e., (13). For Muon, since the right factors ZZ⊤i,ti=1r\Z v_i,t\_i=1^r in (18) are mutually orthogonal and the left factors ii=1r\ u_i\_i=1^r are orthonormal, the matrix sign of (Wt−W⋆)ZZ⊤(W_t-W )Z is given by msgn((Wt−W⋆)ZZ⊤)=∑i:‖ZZ⊤i,t‖2>0i((ZZ⊤i,t)⊤‖ZZ⊤i,t‖2).msgn ((W_t-W )Z )= _i:\|Z v_i,t\|_2>0 u_i ( (Z v_i,t) \|Z v_i,t\|_2 ). For i≥2i≥ 2, ZZ⊤i,t=i,tZZ v_i,t= v_i,t and ‖ZZ⊤i,t‖2=|βi,t|\|Z v_i,t\|_2=| _i,t|, yielding the i≥2i≥ 2 case in (14). For i=1i=1, we have ZZ⊤1,t=κα1,t1+β1,t1Z v_1,t=κ _1,t e_1+ _1,t b_1 and ‖ZZ⊤1,t‖2=κ2α1,t2+β1,t2\|Z v_1,t\|_2= κ^2 _1,t^2+ _1,t^2, so substitution into (8) yields the i=1i=1 case in (14). This also preserves (12). For Newton–Muon, since the right factors i,ti=1r\ v_i,t\_i=1^r are mutually orthogonal and the left factors ii=1r\ u_i\_i=1^r are orthonormal, the matrix sign of Wt−W⋆W_t-W is obtained mode-by-mode: msgn(Wt−W⋆)=∑i:‖i,t‖2>0i(i,t‖i,t‖2)⊤.msgn(W_t-W )= _i:\| v_i,t\|_2>0 u_i ( v_i,t\| v_i,t\|_2 ) . Substituting into (9) gives Wt+1−W⋆=(Wt−W⋆)−ηtmsgn(Wt−W⋆),W_t+1-W =(W_t-W )- _tmsgn(W_t-W ), and matching coefficients yields the i=1i=1 and i≥2i≥ 2 cases in (15). This preserves (12) and completes the proof. ∎ Lemma 1 reduces the matrix iterates to scalar recursions for the coefficients α1,t,βi,t\ _1,t, _i,t\. In particular, convergence of Wt→W⋆W_t→ W is equivalent to α1,t→0 _1,t→ 0 and βi,t→0 _i,t→ 0 for all i=1,…,ri=1,…,r, and the effect of the spike κ is isolated to the mixed (α1,t,β1,t)( _1,t, _1,t) mode. Corollary 2 (Convergence rates). Assume (10) and (11), and let r0>0r_0>0 satisfy |α1,0|≤r0| _1,0|≤ r_0 and |βi,0|≤r0| _i,0|≤ r_0 for all i=1,…,ri=1,…,r. Define rtr_t recursively so that |α1,t|≤rt| _1,t|≤ r_t and |βi,t|≤rt| _i,t|≤ r_t for all i=1,…,ri=1,…,r, and at each step choose ηt _t greedily to minimize the worst-case next bound rt+1r_t+1. Then, for any target 0<ε<r00< <r_0, gradient descent, Newton–Muon, and Muon reach |α1,t|≤ε| _1,t|≤ and |βi,t|≤ε| _i,t|≤ for all i=1,…,ri=1,…,r after at most TGD(ε)=O(κlog(r0/ε))T_GD( )=O(κ (r_0/ )), TNM(ε)=O(log(r0/ε))T_NM( )=O( (r_0/ )), and TM(ε)=O(κlog(r0/ε))T_M( )=O(κ (r_0/ )) iterations, respectively. In particular, Newton–Muon converges faster than both gradient descent and Muon by a factor of order κ. Proof of Corollary 2. For gradient descent, (13) gives |α1,t+1|=|1−ηtκ||α1,t|| _1,t+1|=|1- _tκ|| _1,t| and |βi,t+1|=|1−ηt||βi,t|| _i,t+1|=|1- _t|| _i,t| for i=1,…,ri=1,…,r. Hence, if |α1,t|,|βi,t|≤rt| _1,t|,| _i,t|≤ r_t, then rt+1≤max|1−ηtκ|,|1−ηt|rtr_t+1≤ \|1- _tκ|,\ |1- _t|\r_t. The greedy choice minimizes the right-hand side by solving |1−ηtκ|=|1−ηt||1- _tκ|=|1- _t|, which gives ηt=2/(κ+1) _t=2/(κ+1) and yields rt+1=(κ−1)/(κ+1)rtr_t+1=(κ-1)/(κ+1)r_t. For Newton–Muon, let st≔α1,t2+β1,t2s_t _1,t^2+ _1,t^2. By (15), α1,t+1=(1−ηt/st)α1,t _1,t+1= (1- _t/s_t ) _1,t and β1,t+1=(1−ηt/st)β1,t _1,t+1= (1- _t/s_t ) _1,t. Since |α1,t|,|β1,t|≤rt| _1,t|,| _1,t|≤ r_t, we have st≤2rts_t≤ 2\,r_t. We first bound α1,t+1 _1,t+1 and β1,t+1 _1,t+1. If ηt≤st _t≤ s_t, we obtain 1−ηt/st≤1−ηt/(2rt)1- _t/s_t≤ 1- _t/( 2\,r_t), so |α1,t+1|≤(1−ηt/(2rt))|α1,t|≤rt−ηt/2| _1,t+1|≤ (1- _t/( 2\,r_t) )| _1,t|≤ r_t- _t/ 2, and similarly |β1,t+1|≤rt−ηt/2| _1,t+1|≤ r_t- _t/ 2. If instead ηt≥st _t≥ s_t, then |1−ηt/st|=ηt/st−1|1- _t/s_t|= _t/s_t-1, and therefore |α1,t+1|=(ηt/st−1)|α1,t|≤(ηt/st)|α1,t|≤ηt| _1,t+1|=( _t/s_t-1)| _1,t|≤( _t/s_t)| _1,t|≤ _t. The same argument gives |β1,t+1|≤ηt| _1,t+1|≤ _t. Combining the two cases, |α1,t+1|,|β1,t+1|≤maxηt,rt−ηt/2| _1,t+1|,| _1,t+1|≤ \ _t,\ r_t- _t/ 2\. Next, for i=2,…,ri=2,…,r, (15) gives |βi,t+1|=||βi,t|−ηt|| _i,t+1|= || _i,t|- _t |. Since |βi,t|≤rt| _i,t|≤ r_t, this implies |βi,t+1|≤maxηt,rt−ηt≤maxηt,rt−ηt/2| _i,t+1|≤ \ _t,\ r_t- _t\≤ \ _t,\ r_t- _t/ 2\. Hence rt+1≤maxηt,rt−ηt/2r_t+1≤ \ _t,\ r_t- _t/ 2\. The greedy choice minimizes the right-hand side by equalizing the two terms: ηt=rt−ηt/2 _t=r_t- _t/ 2. Solving gives ηt=(2−2)rt _t=(2- 2)\,r_t, and therefore rt+1=(2−2)rtr_t+1=(2- 2)r_t. For Muon, let st≔κ2α1,t2+β1,t2s_t κ^2 _1,t^2+ _1,t^2. By (14), we have α1,t+1=(1−ηtκ/st)α1,t _1,t+1=(1- _tκ/s_t) _1,t and β1,t+1=(1−ηt/st)β1,t _1,t+1=(1- _t/s_t) _1,t. Since |α1,t|,|β1,t|≤rt| _1,t|,| _1,t|≤ r_t, we have st≤κ2+1rts_t≤ κ^2+1\,r_t. Applying the same argument as above, we have rt+1≤maxηt,rt−ηt/κ2+1r_t+1≤ \ _t,\ r_t- _t/ κ^2+1\. The greedy choice equalizes the two terms, so ηt=rt−ηt/κ2+1 _t=r_t- _t/ κ^2+1, and therefore rt+1=κ2+1/(κ2+1+1)rtr_t+1= κ^2+1/( κ^2+1+1)r_t. Iterating the three recursions shows that gradient descent, Newton–Muon, and Muon all decrease rtr_t geometrically. Solving rt≤εr_t≤ gives TGD(ε)=O(log(r0/ε)/log((κ+1)/(κ−1)))T_GD( )=O ( (r_0/ )/ ((κ+1)/(κ-1)) ), TNM(ε)=O(log(r0/ε))T_NM( )=O( (r_0/ )), and TM(ε)=O(log(r0/ε)/log(1+1/κ2+1))T_M( )=O ( (r_0/ )/ (1+1/ κ^2+1) ). Finally, for large κ, log(κ+1κ−1)≍1κ,log(1+1κ2+1)≍1κ. ( κ+1κ-1 ) 1κ, (1+ 1 κ^2+1 ) 1κ. Therefore TGD(ε)=O(κlog(r0/ε))T_GD( )=O(κ (r_0/ )), TNM(ε)=O(log(r0/ε))T_NM( )=O( (r_0/ )), and TM(ε)=O(κlog(r0/ε))T_M( )=O(κ (r_0/ )). ∎ 4 One-Step Analysis of Newton–Muon The convergence analysis in Section 3 relies on a simple quadratic case study in order to obtain explicit global dynamics. In this section, we provide a complementary one-step analysis that can be carried out more generally through local approximation. It has two parts: (i) we introduce a scale-invariant quadratic score that measures the best achievable one-step decrease along a chosen direction after an optimal line search; and (i) we use a simple spiked activation model to obtain a concrete numerical study. The theoretical derivation for isotropic activations is deferred to Appendix E. 4.1 Quadratic Score We compare several update directions by their predicted one-step improvement under the surrogate (2). Note first that the objective (2) is a quadratic function of the update direction QQ, and we only need to choose a direction QQ and then perform a one-dimensional line search along that direction. This reduces the comparison to a direction choice followed by one-dimensional line search. Concretely, for a fixed direction Q∈ℝm×nQ ^m× n and step size η∈ℝη , consider the update W↦W−ηQW W-ηQ. Plugging ηQηQ into the quadratic surrogate (2) yields the one-dimensional model J(ηQ)=−ηtr(QG⊤)+η22Ntr(HQ(ZZ⊤)Q⊤).J(ηQ)=- (QG )+ η^22Ntr (HQ(ZZ )Q ). Under Assumption 1, the optimal step size is η⋆=tr(QG⊤)/tr(HQ(ZZ⊤/N)Q⊤)η =tr(QG )/tr (HQ(ZZ /N)Q ). Thus, we define the score of a direction QQ to be proportional to the best one-step loss decrease under the surrogate s(Q)≔−2(J(η⋆Q)−J(0))=tr(QG⊤)2tr(HQ(ZZ⊤/N)Q⊤).s(Q) -2 (J(η Q)-J(0) )= tr(QG )^2tr (HQ(ZZ /N)Q ). (19) This provides a general, scale-invariant criterion for comparing different update directions. 4.2 Numerical Study with Spiked Activation We numerically evaluate the score s(Q)s(Q) in (19) for the six directions under G=H(W−W⋆)(ZZ⊤/N)G=H(W-W )(ZZ /N). These directions are (i) Q=GQ=G, (i) QMuon-SVD=msgn(G)Q_Muon -SVD=msgn(G) (exact SVD), (i) QMuon-NSQ_Muon -NS, obtained by applying five Newton–Schulz iterations to GG, specifically, QMuon-NS=X5Q_Muon -NS=X_5, where X0=G/‖G‖FX_0=G/\|G\|_F and Xt+1=(3.4445Im−4.7750XtXt⊤+2.0315(XtXt⊤)2)XtX_t+1= (3.4445I_m-4.7750X_tX_t +2.0315(X_tX_t )^2 )X_t for t=0,…,4t=0,…,4, (iv) QNewton-Muon-SVD=msgn(G(ZZ⊤)−1)Q_Newton -Muon-SVD=msgn (G(ZZ )^-1 ) (exact SVD), (v) QNewton-Muon-NSQ_Newton -Muon-NS, obtained by applying the same five Newton–Schulz iterations to G(ZZ⊤)−1G(ZZ )^-1, and (vi) QNewton=H−1G(ZZ⊤)−1Q_Newton=H^-1G(ZZ )^-1, which is proportional to W−W⋆W-W . Throughout, we fix the square setting m=n=512m=n=512, λmax=1 _ =1, and λmin=10−4 _ =10^-4. We generate the eigenvalues of HH by the stretched-exponential rule λk=λmaxexp(−τ(k−1)p),τ=log(λmax/λmin)(m−1)p, _k= _ (-τ(k-1)^p ), τ= ( _ / _ )(m-1)^p, so that λm=λmin _m= _ , then sample a random orthogonal PP and set H=Pdiag(λ1,…,λm)P⊤H=Pdiag( _1,…, _m)P . We sample W−W⋆∈ℝm×nW-W ^m× n with i.i.d. (0,1)N(0,1) entries. For the activation matrix, we sample the columns i∈ℝn z_i ^n of ZZ independently from (0,diag(κ,1,…,1))N(0,diag(κ,1,…,1)) with κ=64κ=64, so that the population activation second moment is spiked. We use three choices of (N,p)(N,p): the baseline case (8192,0.3)(8192,0.3), a more top-uniform curvature case (8192,2.4)(8192,2.4), and a smaller-sample case (1024,0.3)(1024,0.3). For each setting, we run 10241024 independent simulations and report the mean score with the 2.5%–97.5% interval. Figure 3 shows substantially higher scores for Newton–Muon than for Muon when both activation anisotropy and curvature anisotropy are strong, with Newton–Muon being closer to the optimal Newton direction. Muon also substantially outperforms gradient descent when the curvature is anisotropic. Figure 3: Numerical study with spiked activation second moment diag(κ,1,…,1)diag(κ,1,…,1) and κ=64κ=64. Top: baseline case (N,p)=(8192,0.3)(N,p)=(8192,0.3). Middle: more uniform curvature (N,p)=(8192,2.4)(N,p)=(8192,2.4). Bottom: smaller sample size (N,p)=(1024,0.3)(N,p)=(1024,0.3). The left column shows the spectrum of HH, and the right column shows the corresponding mean absolute scores s(Q)s(Q). 5 LLM Experiments 5.1 Pretraining Benchmark Records on Modded-NanoGPT Benchmark setup. We first compare Newton–Muon, Muon, and AdamW using historical records from the Modded-NanoGPT speedrun benchmark (Jordan et al., 2024a). This benchmark is designed around a fixed objective: train a GPT-style language model on FineWeb (Penedo et al., 2024) and report the wall-clock time needed to reach a target validation loss under a fixed hardware configuration. The benchmark has two tracks: the short track targets a validation loss of 3.28, while the medium track lowers the target to 2.92. For reproducibility, several important implementation details are presented in Appendix A. Reference records. In this benchmark, each Record #k refers to a historical leaderboard submission together with its released training configuration and runtime log. In our experiments, we use these records for reproduction and comparison. We use Record #4 from the short track as a baseline; this run was submitted shortly after Muon was introduced. To benchmark against a stronger setup, for the short track we adapt a training script that is close in configuration to Record #28. For the medium track, we compare against Record #17. However, note that the later records are heavily optimized around the original Muon update, which makes them difficult to surpass without extensive hyperparameter tuning. Since we only tune the learning rate and Newton–Muon hyperparameters, the baseline comparison should be viewed as the most direct comparison of Newton–Muon, while other comparisons should be interpreted as evaluations under training configurations that were heavily tuned for the original Muon update rather than for Newton–Muon. Algorithm 1 summarizes the implementation of Newton–Muon. Input: For each Muon layer ℓ : input activations Zℓ∈ℝn×NZ_ ^n× N, gradient Gℓ∈ℝm×nG_ ^m× n, running second moment Kℓ∈ℝn×nK_ ^n× n (initialize with 10−3In10^-3I_n), stored inverse Kℓ−1K_ ^-1. Input: Hyperparameters: EWMA coefficient β, ridge scaling γ, refresh interval k. for step=0,1,…=0,1,… do foreach ℓ do if (step+1)modk=0( step+1) k=0 then Kℓ←βKℓ+(1−β)ZℓZℓ⊤/NK_ ←βK_ +(1-β)Z_ Z_ /N // Compute via a symmetric rank-N update (Appendix B.1) γℓ←γ⋅tr(Kℓ)/n _ ←γ·tr(K_ )/n Kℓ−1←(Kℓ+γℓIn)−1K_ ^-1←(K_ + _ I_n)^-1 // Compute via Cholesky inverse (Appendix B.2) or a polynomial iteration with a custom symmetric-output matmul (Appendix B.3) Gℓ←GℓKℓ−1G_ ←G_ K_ ^-1 // Kℓ−1K_ ^-1 is applied to the raw layer gradient, before momentum, weight decay, and the remaining Muon pipeline Apply the standard Muon update using the right-preconditioned gradient GℓG_ Algorithm 1 Newton–Muon Baseline. We begin with the short track Record #4 with a single NVIDIA H100 GPU. This run trains a 124M GPT-2 architecture with 3.25B training tokens. The Muon baseline learning rate is 0.00360.0036. We modified the code to run on a single GPU without touching any training pipeline. For Newton–Muon, we use learning rate 0.00400.0040, EWMA coefficient β=0.95β=0.95, ridge scaling γ=0.2γ=0.2, and refresh interval k=32k=32. For AdamW, we assign the transformer block parameters a reduced AdamW learning rate 0.0005760.000576, selected from a learning-rate sweep. All wall-clock numbers reported below are measured in our environment under this configuration. Table 1 summarizes the final validation losses and total training times. Figure 1 compares the validation loss trajectories. Method Loss Time (s) AdamW 3.3801 7228.4 Muon 3.2793 7314.1 Newton–Muon 3.2611 7443.3 Table 1: Short track Record #4 setting (single H100). Ablation for baseline. Unless otherwise stated, we fix the Newton–Muon settings to learning rate 0.00400.0040, EWMA coefficient β=0.95β=0.95, ridge scaling γ=0.2γ=0.2, and refresh interval k=32k=32, and vary only the ablation parameters. Figure 4 shows a two-dimensional sweep over refresh interval k and EWMA coefficient β. Figure 5 sweeps ridge scaling γ and the learning rate. Figure 4: Refresh ablation for Newton–Muon on short track Record #4. Grouped bar plot over refresh interval k, with one bar per EWMA coefficient β (ridge scaling fixed at γ=0.2γ=0.2 and learning rate fixed at 0.00400.0040). Figure 5: Left: Ridge-scaling ablation for Newton–Muon on short track Record #4. Bar plot over ridge scaling γ with k=32k=32, β=0.95β=0.95, and learning rate 0.00400.0040 fixed. Right: Learning-rate ablation for Newton–Muon on short track Record #4. Bar plot over learning rate with k=32k=32, β=0.95β=0.95, and γ=0.2γ=0.2 fixed. These ablations suggest three patterns. First, Figure 4 shows that very frequent refreshes (k∈1,2k∈\1,2\) consistently underperform, while moderate refresh intervals (k∈16,32,64k∈\16,32,64\) yield the best final losses. Second, ridge scaling is essential to make the right preconditioner numerically well behaved. Figure 5 indicates that overly small ridge values can degrade training slightly since the right-preconditioner may be ill-conditioned. In contrast, a wide mid-range (roughly γ∈[0.02,0.5]γ∈[0.02,0.5] here) performs similarly. Very large ridge values eventually move the update back toward standard Muon, which corresponds to the loss rising toward the baseline 3.2793 as γ increases. Third, the learning-rate sweep in Figure 5 is relatively flat around the best region, and the best Newton–Muon learning rate is very close to the best Muon learning rate. Short Track Record #28. We reproduce a short track configuration similar to a record submitted around the same time as Record #28 on a single NVIDIA L40S GPU. This run trains a ∼ 275M-parameter model with about 670670M training tokens. We use the training script at https://github.com/KellerJordan/modded-nanogpt/blob/9d9dc969c451c87b7ad3c84f807db2c2d9109f41/train_gpt.py. For Newton–Muon, we do not change any other configurations including the learning rate. We use EWMA coefficient β=0.8β=0.8, ridge scaling γ=0.2γ=0.2, and refresh interval k=16k=16; since the run is shorter than Record #4, we refresh more aggressively. For AdamW, we assign the transformer block matrices a reduced AdamW learning rate 0.000550.00055. Because this setting is noisy, we run experiments four times and summarize results in Table 2. Figure 6 compares the validation-loss curves of the second-best run (by final validation loss) from each group. Although Newton–Muon has slightly higher total runtime due to preconditioning overhead (Table 2), its loss-versus-time curve indicates a small advantage in time to reach comparable validation loss (Figure 6). Method Run 1 Run 2 Run 3 Run 4 Avg. loss Avg. time (s) AdamW 3.4677 3.4631 3.4640 3.4566 3.4628 4272.0 Muon 3.2758 3.2777 3.2783 3.2830 3.2787 4305.9 Newton–Muon 3.2733 3.2736 3.2740 3.2745 3.2739 4342.4 Table 2: Short track Record #28 setting (single L40S). Figure 6: Validation loss trajectories for the short track Record #28 setting (single L40S). We plot the second-best run for AdamW, Muon, and Newton–Muon, and show loss versus step (left) and loss versus wall-clock time (right). Medium Track Record #17. We next reproduce the medium track Record #17 setting. This run trains a ∼ 455M-parameter model with about 3.123.12B training tokens. Similarly, for Newton–Muon we keep other configurations fixed. We use EWMA coefficient β=0.9β=0.9, ridge scaling γ=0.2γ=0.2, and refresh interval k=24k=24. We run both Muon and Newton–Muon three times and summarize results in Table 3. AdamW is not included for the medium track comparison. The improvement in final validation loss is only marginal in this setting. Method Run 1 Run 2 Run 3 Avg. loss Muon 2.9190 2.9208 2.9190 2.9196 Newton–Muon 2.9175 2.9191 2.9181 2.9183 Table 3: Medium track Record #17 setting. Wall-clock time is not reported since we did not observe a substantial improvement. 5.2 Quadratic Score We evaluate candidate directions using the same one-dimensional surrogate, but with the exact parameter-space Hessian in the denominator. Given a parameter matrix WW and a candidate direction QQ, we define the score s(Q)=tr(G⊤Q)2vec(Q)⊤ℋWvec(Q).s(Q)= tr(G Q)^2vec(Q) H_Wvec(Q). (20) When vec(Q)⊤ℋWvec(Q)>0vec(Q) H_Wvec(Q)>0, this quantity is, up to a constant factor, the predicted best decrease of the local quadratic model after an optimal line search along QQ. We compute the curvature term vec(Q)⊤ℋWvec(Q)vec(Q) H_Wvec(Q) using exact Hessian-vector products for the batch loss, without any Kronecker approximation. To make scores comparable across methods and layers, we normalize each layer’s raw gradient to unit Frobenius norm G^≔G/‖G‖F G G/\|G\|_F. Each candidate direction is also normalized Q^≔Q/‖Q‖F Q Q/\|Q\|_F. We then report three quantities. First, the alignment term tr(G^⊤Q^)tr( G Q). Second, the curvature term (1/2)vec(Q^)⊤ℋWvec(Q^)(1/2)vec( Q) H_Wvec( Q). Third, the resulting score s(Q^)s( Q). Experimental setup. We use Pythia-70M (Biderman et al., 2023) and evaluate two checkpoints, step 1000 and step 50000. We stream the Pile (Gao et al., 2020) and form a large batch of 20482048 sequences of length 10241024, totaling 1024×16×128=2,097,1521024× 16× 128=2,097,152 tokens. We study four weight matrices in the fourth transformer block: the attention output projection, the attention QKV projection, the MLP expansion, and the MLP contraction. To stabilize computation, we replace ZZ⊤ZZ by the damped matrix Kγ≔ZZ⊤+γInK_γ ZZ +γI_n with γ>0γ>0, starting from a small damping value and multiplying γ by 1010 until KγK_γ admits a stable float64 Cholesky factorization. We compare five directions per layer. The gradient descent direction is Q=G^Q= G. Muon-NS5 and Muon-NS32 apply 5 and 32 Newton–Schulz iterations to G G, respectively. Newton–Muon-NS5 and Newton–Muon-NS32 apply 5 and 32 Newton–Schulz steps to the right-preconditioned gradient. Results. Figure 7 summarizes the results. For these four matrices at the two sampled checkpoints, the ranking is consistent. Newton–Muon achieves the highest score, Muon is intermediate, and the raw gradient is worst. By construction, the gradient direction has the largest alignment with itself, but it also typically has a large curvature term. Both Muon and Newton–Muon substantially reduce this curvature term. Relative to Muon, Newton–Muon usually reduces the curvature term more than it reduces the alignment term, which explains its higher score. The advantage of Newton–Muon over Muon is smaller at step 50000 than at step 1000, and this trend appears across all four matrices. Also, NS5 and NS32 are close, suggesting that a small number of Newton–Schulz steps often suffices. Table 4 reports diagnostics of the input activation second moment computed on the same batch. The matrices are strongly anisotropic, with large diagonal spread, substantial off-diagonal mass, and very large condition numbers. Module n κ(ZZ⊤)κ(ZZ ) dmax/dmind_ /d_ o¯/dmean o/d_mean γ/dmeanγ/d_mean Attn out 512 9.18×1039.18× 10^3 2.84×1012.84× 10^1 9.89×1019.89× 10^1 1.0×10−61.0× 10^-6 Attn QKV 512 2.04×1062.04× 10^6 3.00×1003.00× 10^0 6.59×1016.59× 10^1 1.0×10−61.0× 10^-6 MLP expansion 512 ∞ 3.01×1003.01× 10^0 6.56×1016.56× 10^1 1.0×10−41.0× 10^-4 MLP contraction 2048 – 7.23×1017.23× 10^1 1.88×1021.88× 10^2 1.0×10−61.0× 10^-6 Attn out 512 2.34×1032.34× 10^3 3.02×1013.02× 10^1 7.21×1017.21× 10^1 1.0×10−61.0× 10^-6 Attn QKV 512 9.77×1049.77× 10^4 1.66×1021.66× 10^2 4.22×1014.22× 10^1 1.0×10−61.0× 10^-6 MLP expansion 512 8.22×1068.22× 10^6 1.21×1021.21× 10^2 3.97×1013.97× 10^1 1.0×10−61.0× 10^-6 MLP contraction 2048 – 7.46×1027.46× 10^2 1.88×1021.88× 10^2 1.0×10−61.0× 10^-6 Table 4: Diagnostics of the activation second moment ZZ⊤ZZ for four modules at step 1000 (top block) and step 50000 (bottom block). Here dmind_ , dmeand_mean, and dmaxd_ are the minimum, mean, and maximum diagonal entries of ZZ⊤ZZ , o¯=(1/n)∑i∑j≠i|(ZZ⊤)ij| o=(1/n) _i _j≠ i|(ZZ )_ij|, and κ(ZZ⊤)κ(ZZ ) is the condition number of ZZ⊤ZZ . The value κ(ZZ⊤)κ(ZZ ) is reported when eigenvalues were explicitly computed; “∞” indicates numerical instability in the eigenspectrum estimate (a tiny negative eigenvalue at step 1000), and “–” indicates that the spectrum was not explicitly computed at n=2048n=2048. The last column reports the relative damping level selected by the adaptive Cholesky procedure. Figure 7: Quadratic score comparison on four layer matrices from Pythia at checkpoints step 1000 (top) and step 50000 (bottom), corresponding to the one-step analysis in (20). 6 Discussion In this work, we introduce a triplet quadratic surrogate model that offers a local second-order view of Muon and leads to a new optimizer, Newton–Muon. Under an isotropic proxy, one-step minimization yields the update msgn(G(ZZ⊤)−1)msgn(G(Z )^-1), suggesting that Muon is an implicit Newton method without right preconditioning by the input second moment. Empirically, on our reproductions of historical Modded-NanoGPT speedrun benchmark configurations, Newton–Muon reaches the target validation loss in 6% fewer steps and reduces wall-clock time to that loss by over 4% relative to the Muon baseline. However, several limitations of the current framework remain and point to important directions for future work. First, Newton–Muon relies on the isotropic proxy ΣW∝Im _W I_m. This choice is simple and robust, but it discards potentially useful information about the displacement distribution and may become inaccurate later in training. It is therefore better viewed as a practical proxy than as a faithful model throughout optimization. Indeed, weight matrices in classification models often exhibit low dimensional geometric structure, as suggested by neural collapse (Papyan et al., 2020) and minority collapse (Fang et al., 2021); related phenomena, including linguistic collapse (Wu and Papyan, 2024) and cluster formation in attention dynamics (Geshkovski et al., 2023), have also been observed in transformers. A natural direction for future work is to estimate ΣW _W from training dynamics. With a sufficiently accurate estimate of ΣW _W, one need not explicitly compute its SVD or symmetric square root in order to apply (5); it suffices to obtain a factorization ΣW=MM⊤ _W=MM , and then use MM to implement the same update direction. More details on this factorized implementation are deferred to Appendix F. However, our preliminary experiments in Appendix F with non-identity proxies, such as diagonal and factorized forms, did not improve performance, suggesting that estimating ΣW _W reliably is nontrivial and may introduce feedback bias. Developing stable, cheap estimators for non-isotropic proxies remains an important open problem. Second, the Kronecker approximation of the parameter-space Hessian omits explicit token coupling in transformer attention mechanisms. In general, for token indices t≠st≠ s, one can have Hts=∂2L/(∂t∂s)≠0H_ts=∂^2L/(∂ y_t∂ y_s)≠ 0, but these cross-token curvature terms are discarded when passing from the exact token-coupled formula (21) to the averaged Kronecker approximation (22). We also use a single shared curvature matrix H across all samples and positions, whereas the local output space curvature may vary with sample or token. As a result, the approximation ℋW≈(ZZ⊤/N)⊗HH_W≈(Z /N) H may become less accurate in late training, when the Hessian can deviate further from Kronecker structure. These may help explain why Muon’s advantage decays in very late training (Wen et al., 2025). Future work should develop tractable approximations that incorporate token-coupled curvature and relax the shared H assumption. Third, our experiments computed the activation second moment inverse using a damped Cholesky solve of Kγ=ZZ⊤+γInK_γ=ZZ +γ I_n. In most cases we used the full inverse of KγK_γ, but for MLP contraction matrices of shape d×(4d)d×(4d) we instead used a block-diagonal approximation with four d×d× d blocks, which performed well in practice (see Appendix A). This suggests that structured approximations to KγK_γ may already suffice. Therefore, future work could study structured approximations such as block matrices or low-rank factorizations, which may substantially reduce memory and computation while preserving most of the benefit of Newton–Muon. Future work could also compute the activation second-moment inverse via polynomial iteration (Appendix B.3), which may be faster. It is also natural to replace our heuristic damping γ∝tr(ZZ⊤)/nγ (Z )/n by a shifted Cholesky rule (Fukaya et al., 2020), for example γ∝‖ZZ⊤‖2γ \|Z \|_2 or γ∝‖ZZ⊤‖F/nγ \|Z \|_F/ n. Fourth, future work should evaluate Newton–Muon in distributed training, since we only test it on a single GPU. In particular, it remains important to understand how to compute and invert the activation second moment ZZ⊤ZZ in multi-GPU training. This includes both efficient implementations and approximations that preserve most of the benefit of Newton–Muon. Studying these systems issues at scale is an important direction for future work. Acknowledgments This research was supported in part by the Wharton AI fund. References N. Abreu, N. Vyas, S. Kakade, and D. Morwani (2025) The potential of second-order optimization for LLMs: a study with full Gauss-Newton. arXiv preprint arXiv:2510.09378. Cited by: §1. K. Ahn, N. Amsel, and J. Langford (2025) Dion2: a simple method to shrink matrix in Muon. arXiv preprint arXiv:2512.16928. Cited by: §1.2, §1. N. Amsel, D. Persson, C. Musco, and R. M. Gower (2025) The Polar Express: optimal matrix sign methods and their application to the Muon algorithm. arXiv preprint arXiv:2505.16932. Cited by: §1. J. Bernstein and L. Newhouse (2024) Old optimizer, new norm: an anthology. In OPT 2024: Optimization for Machine Learning, Cited by: §1.2, §1. S. Biderman, H. Schoelkopf, Q. G. Anthony, H. Bradley, K. O’Brien, E. Hallahan, M. A. Khan, S. Purohit, U. S. Prashanth, E. Raff, et al. (2023) Pythia: a suite for analyzing large language models across training and scaling. In International conference on machine learning, p. 2397–2430. Cited by: §5.2. D. E. Carlson, E. Collins, Y. Hsieh, L. Carin, and V. Cevher (2015) Preconditioned spectral descent for deep learning. Advances in neural information processing systems 28. Cited by: §1.2, §1. L. Chen, J. Li, and Q. Liu (2025) Muon optimizes under spectral norm constraints. arXiv preprint arXiv:2506.15054. Cited by: §1.2. M. Crawshaw, C. Modi, M. Liu, and R. M. Gower (2025) An exploration of non-Euclidean gradient descent: Muon and its many variants. arXiv preprint arXiv:2510.09827. Cited by: §1.2. D. Davis and D. Drusvyatskiy (2025) When do spectral gradient updates help in deep learning?. arXiv preprint arXiv:2512.04299. Cited by: §1.2, §3. S. Dragutinović and R. Ranganath (2026) To use or not to use Muon: how simplicity bias in optimizers matters. arXiv preprint arXiv:2603.00742. Cited by: §1.2. 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.2. C. Fang, H. He, Q. Long, and W. J. Su (2021) Exploring deep neural networks via layer-peeled model: minority collapse in imbalanced training. Proceedings of the National Academy of Sciences 118 (43), p. e2103091118. Cited by: §6. T. Fukaya, R. Kannan, Y. Nakatsukasa, Y. Yamamoto, and Y. Yanagisawa (2020) Shifted Cholesky QR for computing the QR factorization of ill-conditioned matrices. SIAM Journal on Scientific Computing 42 (1), p. A477–A503. Cited by: §6. L. Gao, S. Biderman, S. Black, L. Golding, T. Hoppe, C. Foster, J. Phang, H. He, A. Thite, N. Nabeshima, et al. (2020) The Pile: an 800GB dataset of diverse text for language modeling. arXiv preprint arXiv:2101.00027. Cited by: §5.2. T. George, C. Laurent, X. Bouthillier, N. Ballas, and P. Vincent (2018) Fast approximate natural gradient descent in a Kronecker factored eigenbasis. Advances in neural information processing systems 31. Cited by: §1.2. B. Geshkovski, C. Letrouit, Y. Polyanskiy, and P. Rigollet (2023) The emergence of clusters in self-attention dynamics. Advances in Neural Information Processing Systems 36, p. 57026–57037. Cited by: §6. W. Gong, J. Zazo, Q. Luo, P. Wang, J. Hensman, and C. Ma (2026) ARO: a new lens on matrix optimization for large models. arXiv preprint arXiv:2602.09006. Cited by: §1.2, §1. A. Gonon, A. Muşat, and N. Boumal (2026) Insights on Muon from simple quadratics. arXiv preprint arXiv:2602.11948. Cited by: §1.2. Y. Gu and Z. Xie (2026) Mano: restriking manifold optimization for LLM training. arXiv preprint arXiv:2601.23000. Cited by: §1. V. Gupta, T. Koren, and Y. Singer (2018) Shampoo: preconditioned stochastic tensor optimization. In International Conference on Machine Learning, p. 1842–1850. Cited by: §1.2, §1.2. C. He, Z. Deng, and Z. Lu (2025) Low-rank orthogonalization for large-scale matrix optimization with applications to foundation model training. arXiv preprint arXiv:2509.11983. Cited by: §1. K. Jordan, J. Bernstein, B. Rappazzo, @fernbear.bsky.social, B. Vlado, Y. Jiacheng, F. Cesista, B. Koszarsky, and @Grad62304977 (2024a) modded-nanogpt: speedrunning the NanoGPT baseline. External Links: Link Cited by: §B.1, §1, §5.1. K. Jordan, Y. Jin, V. Boza, J. You, F. Cesista, L. Newhouse, and J. Bernstein (2024b) Muon: an optimizer for hidden layers in neural networks. External Links: Link Cited by: §1. D. P. Kingma and J. Ba (2014) Adam: a method for stochastic optimization. arXiv preprint arXiv:1412.6980. Cited by: §1. A. Krizhevsky (2009) Learning multiple layers of features from tiny images. Technical report University of Toronto. Cited by: Appendix C. T. T. Lau, Q. Long, and W. Su (2025) PolarGrad: A class of matrix-gradient optimizers from a unifying preconditioning perspective. arXiv preprint arXiv:2505.21799. Cited by: §1.2, §1.2, §1. J. Li and M. Hong (2025) A note on the convergence of Muon. arXiv preprint arXiv:2502.02900. Cited by: §1.2. Z. Li, L. Liu, C. Liang, W. Chen, and T. Zhao (2025) NorMuon: making Muon more efficient and scalable. arXiv preprint arXiv:2510.05491. Cited by: §1. 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. arXiv preprint arXiv:2502.16982. Cited by: §1. I. Loshchilov and F. Hutter (2017) Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101. Cited by: §1. J. Ma, Y. Huang, Y. Chi, and Y. Chen (2026) Preconditioning benefits of spectral orthogonalization in Muon. arXiv preprint arXiv:2601.13474. Cited by: §1.2. J. Martens and R. Grosse (2015) Optimizing neural networks with Kronecker-factored approximate curvature. In International conference on machine learning, p. 2408–2417. Cited by: §1.2, §1.2, §1, §2.2. D. Morwani, I. Shapira, N. Vyas, E. Malach, S. Kakade, and L. Janson (2024) A new perspective on Shampoo’s preconditioner. arXiv preprint arXiv:2406.17748. Cited by: §1.2. V. Papyan, X. Han, and D. L. Donoho (2020) Prevalence of neural collapse during the terminal phase of deep learning training. Proceedings of the National Academy of Sciences 117 (40), p. 24652–24663. Cited by: §6. A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. (2019) PyTorch: an imperative style, high-performance deep learning library. Advances in neural information processing systems 32. Cited by: §B.1. G. Penedo, H. Kydlíček, L. B. allal, A. Lozhkov, M. Mitchell, C. Raffel, L. V. Werra, and T. Wolf (2024) The FineWeb datasets: decanting the web for the finest text data at scale. In The Thirty-eight Conference on Neural Information Processing Systems Datasets and Benchmarks Track, External Links: Link Cited by: §5.1. T. Pethick, W. Xie, K. Antonakopoulos, Z. Zhu, A. Silveti-Falls, and V. Cevher (2025) Training deep learning models with norm-constrained LMOs. In Proceedings of the International Conference on Machine Learning (ICML), Cited by: §1.2, §1.2, §1, §1. X. Qi, M. Chen, J. Ye, Y. He, and R. Xiao (2026) Delving into Muon and beyond: deep analysis and extensions. arXiv preprint arXiv:2602.04669. Cited by: §1. A. Riabinin, E. Shulgin, K. Gruntkowska, and P. Richtárik (2025) Gluon: making Muon & Scion great again!(bridging theory and practice of LMO-based optimizers for LLMs). arXiv preprint arXiv:2505.13416. Cited by: §1.2. N. Sato, H. Naganuma, and H. Iiduka (2025) Analysis of Muon’s convergence and critical batch size. arXiv preprint arXiv:2507.01598. Cited by: §1.2. I. Shah, A. M. Polloreno, K. Stratos, P. Monk, A. Chaluvaraju, A. Hojel, A. Ma, A. Thomas, A. Tanwer, D. J. Shah, et al. (2025) Practical efficiency of Muon for pretraining. arXiv preprint arXiv:2505.02222. Cited by: §1. N. Shazeer and M. Stern (2018) Adafactor: adaptive learning rates with sublinear memory cost. In International conference on machine learning, p. 4596–4604. Cited by: §1.2. J. W. Silverstein and Z. D. Bai (1995) On the empirical distribution of eigenvalues of a class of large dimensional random matrices. Journal of Multivariate analysis 54 (2), p. 175–192. Cited by: §E.2. D. Su, A. Gu, J. Xu, Y. Tian, and J. Zhao (2025) GaLore 2: large-scale LLM pre-training by gradient low-rank projection. arXiv preprint arXiv:2504.20437. Cited by: §1.2. W. Su (2025) Isotropic curvature model for understanding deep learning optimization: is gradient orthogonalization optimal?. arXiv preprint arXiv:2511.00674. Cited by: §1.2, §1, §2.1. K. Team, T. Bai, Y. Bai, Y. Bao, S. H. Cai, Y. Cao, Y. Charles, H. S. Che, et al. (2026) Kimi K2.5: visual agentic intelligence. External Links: 2602.02276, Link Cited by: §1. P. Tillet, H. Kung, and D. Cox (2019) Triton: an intermediate language and compiler for tiled neural network computations. In Proceedings of the 3rd ACM SIGPLAN International Workshop on Machine Learning and Programming Languages, p. 10–19. Cited by: §B.1. N. Vyas, D. Morwani, R. Zhao, I. Shapira, D. Brandfonbrener, L. Janson, and S. Kakade (2025) SOAP: Improving and stabilizing Shampoo using Adam. In International Conference on Learning Representations (ICLR), Cited by: §1.2. K. Wen, D. Hall, T. Ma, and P. Liang (2025) Fantastic pretraining optimizers and where to find them. arXiv preprint arXiv:2509.02046. Cited by: §1, §6. R. Wu and V. Papyan (2024) Linguistic collapse: neural collapse in (large) language models. Advances in Neural Information Processing Systems 37, p. 137432–137473. Cited by: §6. C. Xu, W. Yan, and Y. A. Zhang (2026) FISMO: Fisher-structured momentum-orthogonalized optimizer. arXiv preprint arXiv:2601.21750. Cited by: §1. S. Yang, Z. Wang, O. Balabanov, N. B. Erichson, and M. W. Mahoney (2026) PRISM: distribution-free adaptive computation of matrix functions for accelerating neural network training. arXiv preprint arXiv:2601.22137. Cited by: §1.2. D. Yu, H. Tao, Y. Wan, L. Luo, and L. Zhang (2026) Sign-based optimizers are effective under heavy-tailed noise. arXiv preprint arXiv:2602.07425. Cited by: §1.2. H. Yuan, Y. Liu, S. Wu, X. Zhou, and Q. Gu (2024) Mars: unleashing the power of variance reduction for training large models. arXiv preprint arXiv:2411.10438. Cited by: §1.2. A. Zeng, X. Lv, Z. Hou, Z. Du, Q. Zheng, B. Chen, D. Yin, C. Ge, C. Xie, C. Wang, et al. (2026) GLM-5: from vibe coding to agentic engineering. arXiv preprint arXiv:2602.15763. Cited by: §1. J. Zhang, N. Amsel, B. Chen, and T. Dao (2026) Gram Newton-Schulz. External Links: Link Cited by: §1. J. Zhao, Z. Zhang, B. Chen, Z. Wang, A. Anandkumar, and Y. Tian (2024) GaLore: memory-efficient LLM training by gradient low-rank projection. arXiv preprint arXiv:2403.03507. Cited by: §1.2. Appendix A LLM Experimental Details This appendix collects the implementation and benchmarking details for the LLM experiments in Section 5.1, including several details that are important for reproducing the reported results. Comparison protocol. For a fair comparison, we use several publicly logged benchmark records as configuration references and report wall-clock results from our reproductions. To keep the comparison controlled, for each selected record we modify only the optimizer and tune learning rates when needed. For Newton–Muon, we replace the Muon update with Newton–Muon. For AdamW, to keep the comparison controlled, we replace the Muon-optimized parameters by AdamW, set β=(0.9,0.999)β=(0.9,0.999), and tune the AdamW learning rate over a wide range, reporting the setting that achieves the lowest validation loss, while keeping the remaining training settings unchanged. All other training settings are kept identical to the original record; we may make implementation-level changes (e.g., custom kernels) only to improve efficiency without altering the underlying training pipeline, and we apply any such changes uniformly across methods within the same comparison. For each selected record, we reproduce the run in our own environment and report wall-clock time from these reproductions. We do not claim that our measured times match the leaderboard submission times; our focus is the within-environment wall-clock difference among Newton–Muon, Muon, and AdamW under otherwise identical training settings. Layerwise second moment computation. For the attention QKV projection and the MLP expansion, the input dimension is d, so we form a single d×d× d activation second moment matrix from their inputs. For the attention output projection, the input dimension is also d, but the activations are taken from the attention output immediately before applying the output projection. For the MLP contraction, the input dimension is 4d4d. To avoid forming a full (4d)×(4d)(4d)×(4d) second moment, we use a block-diagonal approximation: reshape each activation ∈ℝ4d z ^4d into four contiguous blocks =[(1);…;(4)] z=[ z^(1);…; z^(4)] with (b)∈ℝd z^(b) ^d, and form four separate d×d× d second moments Z(b)Z(b)⊤Z^(b)Z^(b) . The corresponding inverse is represented as four independent d×d× d inverses, and applying it is done by splitting the MLP contraction gradient into four d×d× d blocks along the input dimension and right-multiplying each block by its matching inverse. Numerical precision. There are also a few important implementation details regarding numerical precision. First, the activation second moment update ZℓZℓ⊤Z_ Z_ is relatively well behaved, and in practice the accumulation can be done with bfloat16. However, the inverse computation is much more sensitive. We compute Kℓ−1=(Kℓ+γℓIn)−1K_ ^-1=(K_ + _ I_n)^-1 using a float32 Cholesky factorization followed by a Cholesky inverse (see Appendix B.2). Second, when applying the right-preconditioner, the matrix multiply GℓKℓ−1G_ K_ ^-1 should also be performed in float32. If the gradients are stored in bfloat16, we can upcast GℓG_ to float32 before the multiplication. Appendix B Optimizing Computation B.1 Symmetric Matrix Multiplication Our implementation exploits the fact that the activation second moment ZZ⊤ZZ (and any polynomial p(ZZ⊤)p(ZZ )) is symmetric. This structure allows us to avoid redundant computation by computing only one triangle and reconstructing the other by symmetry. SYRK (symmetric rank-k update). Forming the second moment K←ZZ⊤K←ZZ is exactly a SYRK pattern: it suffices to compute either the lower or upper triangle of KK and fill the other half by mirroring across the diagonal. SYPP (Symmetric polynomial product). In the polynomial inverse, we repeatedly multiply factors that are polynomials in ZZ⊤ZZ . Since ZZ⊤ZZ is symmetric, any polynomial p(ZZ⊤)p(ZZ ) is also symmetric. Moreover, because these factors are functions of the same matrix, they commute, so products such as p(ZZ⊤)q(ZZ⊤)p(ZZ )q(ZZ ) remain symmetric. Consequently, we can compute only one triangle of the product and write the mirrored entries to complete the matrix, reducing nearly half of the work compared to a dense matrix multiplication. PyTorch (Paszke et al., 2019) does not expose a simple interface for enforcing this kind of triangular compute. Instead, we implement Triton (Tillet et al., 2019) custom kernels, adapted from the Modded-NanoGPT repository (Jordan et al., 2024a). B.2 Cholesky Inverse When Kγ=ZZ⊤+γIn≻0K_γ=ZZ +γI_n 0, Cholesky factorization gives Kγ=LL⊤K_γ=LL with LL lower triangular. We then explicitly form the inverse Kγ−1=(LL⊤)−1=L−⊤L−1K_γ^-1=(LL )^-1=L^- L^-1. We use torch.linalg.cholesky_ex followed by torch.cholesky_inverse in float32. B.3 Polynomial Iteration Inverse We also consider a polynomial iteration that explicitly constructs an approximation to Kγ−1K_γ^-1 using only SYPP. When KγK_γ has moderate condition number after damping, this approach can be cheaper than Cholesky. Principle. Let Kγ≻0K_γ 0 and choose a scalar α>0α>0 such that the scaled matrix K~≔αKγ K αK_γ has spectrum contained in (ε,1]( ,1], i.e., ε<λmin(K~)≤λmax(K~)≤1 < _ ( K)≤ _ ( K)≤ 1 for a small safety margin ε>0 >0. Define the residual matrix R0≔In−K~.R_0 I_n- K. With this scaling, the spectrum of R0R_0 lies in [0,1−ε][0,1- ]. Since K~=In−R0 K=I_n-R_0, Kγ−1=αK~−1=α(In−R0)−1.K_γ^-1=α K^-1=α(I_n-R_0)^-1. The goal is therefore to build an explicit approximation X≈(In−R0)−1X≈(I_n-R_0)^-1. We initialize the inverse estimate as X0=In.X_0=I_n. At step k=1,…,Tk=1,…,T, we pick a polynomial qk(⋅)q_k(·) and apply it to the current residual matrix Qk≔qk(Rk−1),Xk≔Xk−1Qk.Q_k q_k(R_k-1), X_k X_k-1Q_k. This induces a residual update. Using Rk−1=In−K~Xk−1R_k-1=I_n- KX_k-1, we have K~Xk−1=In−Rk−1 KX_k-1=I_n-R_k-1, hence Rk=In−K~Xk=In−K~Xk−1Qk=In−(In−Rk−1)qk(Rk−1)=ϕk(Rk−1),R_k=I_n- KX_k=I_n- KX_k-1Q_k=I_n-(I_n-R_k-1)q_k(R_k-1)= _k(R_k-1), where the induced scalar map on eigenvalues is r+=ϕk(r)≔1−(1−r)qk(r).r^+= _k(r) 1-(1-r)q_k(r). Crucially, the admissible interval for the residual is iteration dependent. If Rk−1R_k-1 has spectrum contained in an interval r∈ℐk−1r _k-1 (ℐ0=[0,1−ε]I_0=[0,1- ]), then the next spectrum is contained in r+∈ϕk(ℐk−1),r^+∈ _k(I_k-1), and the interval shrinks as k grows. As the residual contracts toward zero, the effective design interval becomes much smaller than [0,1−ε][0,1- ], which allows later polynomials to be optimized for a tighter range and achieve stronger contraction per SYPP. Convergence guarantee. Define s0≔1−εs_0 1- , so that the initial residual interval is ℐ0=[0,s0]I_0=[0,s_0]. If the first polynomial q1q_1 is chosen so that supr∈ℐ0|ϕ1(r)|≤s1<s0, _r _0| _1(r)|≤ s_1<s_0, and for each later step k≥2k≥ 2 the polynomial qkq_k is chosen so that sup|r|≤sk−1|ϕk(r)|≤sk<sk−1, _|r|≤ s_k-1| _k(r)|≤ s_k<s_k-1, then ‖Rk‖2≤sk\|R_k\|_2≤ s_k and the residual interval shrinks monotonically to 0. Repeating a short sequence of such contractive steps yields RTR_T close to zero and XTX_T as an explicit approximation to (In−R0)−1=K~−1(I_n-R_0)^-1= K^-1. Finally, Kγ−1≈αXT.K_γ^-1≈αX_T. How the polynomials are chosen. At the first step we choose q1q_1 to minimize s1=supr∈ℐ0|ϕ1(r)|s_1= _r _0| _1(r)| on the one-sided interval ℐ0=[0,s0]I_0=[0,s_0]. For each later step k≥2k≥ 2, we take the current symmetric spectral bound sk−1s_k-1 and choose qkq_k to minimize the next bound sk=sup|r|≤sk−1|ϕk(r)|s_k= _|r|≤ s_k-1| _k(r)| under coefficient-magnitude constraints and an explicit robustness margin that accounts for finite precision and modeling error. We implement this minimax design on a grid (via a linear program) and then chain multiple steps under a total SYPP budget (via dynamic programming). The first step is designed on a one-sided interval, and subsequent steps are designed on the symmetric interval to remain stable once roundoff introduces small negative eigenvalues. Numerical safeguards for precision. Two hyperparameters in the LP directly control how conservative the bound update sk−1↦sks_k-1 s_k is. First, INTERVAL_PAD_REL enlarges the design interval. For the first step, instead of optimizing ϕ1 _1 only on ℐ0=[0,s0]I_0=[0,s_0], the solver designs on ℐ0design=[0,(1+INTERVAL_PAD_REL)s0].I_0^design=[0,(1+ INTERVAL\_PAD\_REL)s_0]. For each later step k≥2k≥ 2, instead of optimizing ϕk _k only on the current spectral bound |r|≤sk−1|r|≤ s_k-1, the solver designs on |r|≤skdesign|r|≤ s_k^design with skdesign=(1+INTERVAL_PAD_REL)sk−1.s_k^design=(1+ INTERVAL\_PAD\_REL)s_k-1. This guards against underestimating the current residual bound by enforcing contraction on a slightly wider interval than the nominal bound. Second, NOISE_ABS is a worst-case absolute perturbation radius for the scalar residual map. Concretely, the LP is solved under the robust requirement that for all r in the design interval and for all perturbations δ∈[−NOISE_ABS,NOISE_ABS]δ∈[- NOISE\_ABS, NOISE\_ABS], the perturbed update remains bounded |ϕk(r)+δ|≤γk. | _k(r)+δ |≤ _k. Equivalently, the solver enforces |ϕk(r)|≤γk−NOISE_ABS| _k(r)|≤ _k- NOISE\_ABS on the grid, and then reports the robust bound as sk=γk=sup|ϕk(r)|+NOISE_ABSs_k= _k= | _k(r)|+ NOISE\_ABS. The algorithm is shown in Algorithm 2. Input: positive semidefinite matrix K∈ℝn×nK ^n× n, damping γ>0γ>0, SYPP budget B. Output: Explicit approximation Kγ−1~≈(K+γIn)−1 K_γ^-1≈(K+γI_n)^-1. Kγ←K+γInK_γ←K+γI_n Estimate any upper bound of the largest eigenvalue λ¯max(K)≥λmax(K) λ_ (K)≥ _ (K) λ¯max(Kγ)←λ¯max(K)+γ λ_ (K_γ)← λ_ (K)+γ α←1/λ¯max(Kγ)α← 1/ λ_ (K_γ), K~←αKγ K←αK_γ, ε¯←αγ ←αγ, X0←InX_0←I_n, R0←In−K~R_0←I_n- K Choose a plan (q1,…,qT)(q_1,…,q_T) from Table 5 with tabulated ε≤ε¯ ≤ and total SYPP cost ≤B≤ B for k←1k← 1 to T−1T-1 do Qk←qk(Rk−1)Q_k← q_k(R_k-1) Xk←Xk−1QkX_k←X_k-1Q_k Rk←In−(In−Rk−1)QkR_k←I_n-(I_n-R_k-1)Q_k QT←qT(RT−1)Q_T← q_T(R_T-1) XT←XT−1QTX_T←X_T-1Q_T return Kγ−1~←αXT K_γ^-1←αX_T such that ‖In−KγKγ−1~‖2≤sout \|I_n-K_γ K_γ^-1 \|_2≤ s_out, where souts_out is the certified residual bound associated with the selected plan Algorithm 2 Inverse of Kγ=K+γInK_γ=K+γI_n via polynomial approximation ε Total souts_out SYPP Polynomial qk(x)q_k(x) 0.0015 12 0.030717 3 q1(x)=1.991037−15.856588x+31.760959x2q_1(x)=1.991037-15.856588x+31.760959x^2 4 q2(x)=0.102569+0.102569x+7.383161x2+7.383161x3q_2(x)=0.102569+0.102569x+7.383161x^2+7.383161x^3 3 q3(x)=1+2.541910x+2.541910x2q_3(x)=1+2.541910x+2.541910x^2 2 q4(x)=1+1.192261x+1.192261x2q_4(x)=1+1.192261x+1.192261x^2 0.0015 13 0.004865 3 q1(x)=1.991037−15.856588x+31.760959x2q_1(x)=1.991037-15.856588x+31.760959x^2 3 q2(x)=1+3.839962x+3.839963x2q_2(x)=1+3.839962x+3.839963x^2 3 q3(x)=1+2.989700x+2.989700x2q_3(x)=1+2.989700x+2.989700x^2 2 q4(x)=1.244063+1.244063xq_4(x)=1.244063+1.244063x 2 q5(x)=1+1.047265x+1.047265x2q_5(x)=1+1.047265x+1.047265x^2 0.003 10 0.019885 3 q1(x)=1.964953−15.439061x+31.064790x2q_1(x)=1.964953-15.439061x+31.064790x^2 3 q2(x)=1+3.346712x+3.346712x2q_2(x)=1+3.346712x+3.346712x^2 2 q3(x)=1.403255+1.403255xq_3(x)=1.403255+1.403255x 2 q4(x)=1+1.140006x+1.140006x2q_4(x)=1+1.140006x+1.140006x^2 0.003 11 0.002839 3 q1(x)=1.964953−15.439061x+31.064790x2q_1(x)=1.964953-15.439061x+31.064790x^2 3 q2(x)=1+3.346712x+3.346712x2q_2(x)=1+3.346712x+3.346712x^2 3 q3(x)=1+1.757644x+1.757644x2q_3(x)=1+1.757644x+1.757644x^2 2 q4(x)=1+1.028634x+1.028634x2q_4(x)=1+1.028634x+1.028634x^2 0.006 8 0.047094 3 q1(x)=1.915935−14.653845x+29.754543x2q_1(x)=1.915935-14.653845x+29.754543x^2 3 q2(x)=1+2.716205x+2.716205x2q_2(x)=1+2.716205x+2.716205x^2 2 q3(x)=1+1.262596x+1.262596x2q_3(x)=1+1.262596x+1.262596x^2 0.006 9 0.014106 3 q1(x)=1.915935−14.653845x+29.754543x2q_1(x)=1.915935-14.653845x+29.754543x^2 4 q2(x)=0.639753+0.639753x+4.060692x2+4.060692x3q_2(x)=0.639753+0.639753x+4.060692x^2+4.060692x^3 2 q3(x)=1+1.108734x+1.108734x2q_3(x)=1+1.108734x+1.108734x^2 0.006 10 0.002087 3 q1(x)=1.915935−14.653845x+29.754543x2q_1(x)=1.915935-14.653845x+29.754543x^2 3 q2(x)=1+2.716205x+2.716205x2q_2(x)=1+2.716205x+2.716205x^2 2 q3(x)=1.160973+1.160973xq_3(x)=1.160973+1.160973x 2 q4(x)=1+1.020113x+1.020111x2q_4(x)=1+1.020113x+1.020111x^2 0.012 7 0.047594 3 q1(x)=1.828900−13.257429x+27.420730x2q_1(x)=1.828900-13.257429x+27.420730x^2 3 q2(x)=1+2.072900x+2.072900x2q_2(x)=1+2.072900x+2.072900x^2 1 q3(x)=1.046594+1.046594xq_3(x)=1.046594+1.046594x 0.012 8 0.008118 3 q1(x)=1.828900−13.257429x+27.420730x2q_1(x)=1.828900-13.257429x+27.420730x^2 3 q2(x)=1+2.072900x+2.072900x2q_2(x)=1+2.072900x+2.072900x^2 2 q3(x)=1+1.071558x+1.071558x2q_3(x)=1+1.071558x+1.071558x^2 0.025 6 0.048057 4 q1(x)=1.528164+1.400800x−12.902311x2+32x4q_1(x)=1.528164+1.400800x-12.902311x^2+32x^4 2 q2(x)=1+1.266514x+1.266514x2q_2(x)=1+1.266514x+1.266514x^2 0.025 7 0.008458 3 q1(x)=1.679044−10.844625x+23.373926x2q_1(x)=1.679044-10.844625x+23.373926x^2 2 q2(x)=1.301562+1.301562xq_2(x)=1.301562+1.301562x 2 q3(x)=1+1.073878x+1.073878x2q_3(x)=1+1.073878x+1.073878x^2 Table 5: Polynomial plans with CMAX=32 (the maximum absolute polynomial coefficient), INTERVAL_PAD_REL=0.001, and NOISE_ABS=0.001. Total denotes the total number of symmetric polynomial products (SYPP) used by the plan, and the SYPP column reports the SYPP cost of each polynomial qkq_k. The quantity souts_out is the certified final residual bound, i.e., an upper bound on ‖In−KγKγ−1~‖2 \|I_n-K_γ K_γ^-1 \|_2. Appendix C CIFAR-10 Experiment Dataset and model. This appendix describes the CIFAR-10 (Krizhevsky, 2009) experiment shown in the second row of Figure 1, where we report test accuracy versus training step and training time. The 50,000 training images are split into 45,000 training examples and 5,000 validation examples, and the standard 10,000-image test set is used for final evaluation. Training augmentation consists of random cropping with padding 4 and random horizontal flipping. All images are normalized channelwise using the standard CIFAR-10 mean and standard deviation. The model is a residual MLP with 32 hidden layers of width 512. Each hidden layer consists of a linear map followed by LayerNorm and GELU. Residual connections are used whenever the input and output dimensions match. Thus, the first layer 3072→5123072→ 512 has no skip connection, while the remaining 31 hidden layers of shape 512→512512→ 512 use residual additions. The output layer is a linear classifier from 512 to 10. Training setup. All three methods are trained for 100 epochs with batch size 4096 on an A100 GPU. The learning-rate schedule is linear warmup followed by cosine decay, with 100 warmup steps and minimum learning-rate ratio 0.1. Validation is evaluated every 24 training steps during hyperparameter tuning. For Muon and Newton–Muon, only the hidden-layer weight matrices are assigned to the matrix optimizer, while the remaining parameters are optimized by AdamW. For the pure AdamW baseline, all trainable parameters are optimized by AdamW. We tuned the hyperparameters of all methods on the validation split. We then retrained each method with its selected hyperparameters on the full 50,000-image CIFAR-10 training set and report test accuracy versus training step and training time in Figure 1. Final hyperparameters. The final AdamW baseline uses learning rate 8×10−48× 10^-4, weight decay 10−210^-2, and (β1,β2)=(0.9,0.999)( _1, _2)=(0.9,0.999). The final Muon configuration uses AdamW on the non-matrix parameters with learning rate 1.6×10−31.6× 10^-3, weight decay 10−210^-2, and (β1,β2)=(0.9,0.999)( _1, _2)=(0.9,0.999), and applies Muon to the hidden-layer weight matrices with matrix learning rate 0.160.16, matrix weight decay 10−310^-3, momentum 0.80.8. The final Newton–Muon configuration uses AdamW on the non-matrix parameters with learning rate 8×10−48× 10^-4, weight decay 10−210^-2, and (β1,β2)=(0.9,0.999)( _1, _2)=(0.9,0.999), and applies Newton–Muon to the hidden-layer weight matrices with matrix learning rate 0.160.16, matrix weight decay 3×10−43× 10^-4, momentum 0.750.75, together with EWMA β=0.95β=0.95, ridge γ=0.05γ=0.05, and refresh interval k=16k=16. Appendix D Kronecker-Factored Curvature Here we derive the Kronecker-factored approximation in (4). In a transformer with sequence length N, let t∈ℝn z_t ^n denote the input activation at token t and let t∈ℝm y_t ^m denote the output of a linear map t=Wt y_t=W z_t for t=1,…,Nt=1,…,N. For an update direction Q∈ℝm×nQ ^m× n, the output perturbation at token t is Δt=Qt y_t=Q z_t. To expose the Kronecker structure, first consider a single token loss Lt(t)L_t( y_t) with output-space curvature Ht≔∇t2Lt(t)∈ℝm×mH_t _ y_t^2L_t( y_t) ^m× m. Then the exact second-order change satisfies δ2Lt=(Δt)⊤Ht(Δt)=(Qt)⊤Ht(Qt).δ^2L_t=( y_t) H_t( y_t)=(Q z_t) H_t(Q z_t). Using vec(Qt)=(t⊤⊗Im)vec(Q)vec(Q z_t)=( z_t I_m)vec(Q), we obtain (Qt)⊤Ht(Qt)=vec(Q)⊤((tt⊤)⊗Ht)vec(Q).(Q z_t) H_t(Q z_t)=vec(Q) (( z_t z_t ) H_t )vec(Q). For token-coupled losses, stack token outputs into ≔[1⊤,…,N⊤]⊤∈ℝmN y [ y_1 ,…, y_N ] ^mN, and let the averaged scalar loss be f()≔L()/Nf( y) L( y)/N, where L()L( y) denotes the summed loss over the N tokens. The token-coupled output-space curvature of L is ℋ≔∇2L()∈ℝ(mN)×(mN). H _ y^2L( y) ^(mN)×(mN). Write ℋ H in N×N× N blocks of size m×m× m: ℋ=[H11⋯H1N⋮⋱⋮HN1⋯HNN],Hts∈ℝm×m. H= bmatrixH_11&·s&H_1N\\ & & \\ H_N1&·s&H_N bmatrix, H_ts ^m× m. Then δ2f=1N∑t=1N∑s=1N(Qt)⊤Hts(Qs).δ^2f= 1N _t=1^N _s=1^N(Q z_t) H_ts(Q z_s). Each term can be written as (Qt)⊤Hts(Qs)=vec(Q)⊤((ts⊤)⊗Hts)vec(Q),(Q z_t) H_ts(Q z_s)=vec(Q) (( z_t z_s ) H_ts )vec(Q), so the exact parameter-space Hessian for the token-coupled loss is ℋW=1N∑t=1N∑s=1N(ts⊤)⊗Hts∈ℝ(mn)×(mn).H_W= 1N _t=1^N _s=1^N( z_t z_s ) H_ts ^(mn)×(mn). (21) Define the input activation matrix Z=[1,…,N]∈ℝn×NZ=[ z_1,…, z_N] ^n× N, the diagonal-block average curvature H≔1N∑t=1NHtt.H 1N _t=1^NH_t. A simplifying approximation retains only diagonal-token curvature contributions and replaces HttH_t by their average HH, yielding ℋW≈1N∑t=1N(tt⊤)⊗H=(ZZ⊤/N)⊗H.H_W≈ 1N _t=1^N( z_t z_t ) H=(ZZ /N) H. (22) Remark. For architectures where loss decomposes across tokens/examples such as MLPs, the parameter-space Hessian (21) is block-diagonal, i.e., Hts=H_ts=0 for t≠st≠ s. In this case, (21) reduces to the diagonal-token form ℋW=1N∑t=1N(tt⊤)⊗Htt,H_W= 1N _t=1^N( z_t z_t ) H_t, so the only approximation in (22) comes from replacing the varying blocks HttH_t by their average HH. As a result, the Kronecker estimator (22) is expected to be more accurate when the loss is not token-coupled. Appendix E Theoretical Quadratic Score Under Isotropic Activation This appendix develops the theoretical quadratic score analysis under an isotropic activation second moment, deriving explicit formulas and approximations for gradient descent, Muon, Newton–Muon, and Newton, and then using them to interpret the corresponding numerical study. E.1 Score Formulas for Different Directions We use (19) to compare four update directions under the quadratic surrogate: (i) the raw GD direction Q=GQ=G, (i) the Muon direction QMuon=msgn(G)Q_Muon=msgn(G), (i) the Newton–Muon direction QNewton-Muon=msgn(G(ZZ⊤)−1)Q_Newton -Muon=msgn (G(ZZ )^-1 ), and (iv) the Newton direction QNewton=H−1G(ZZ⊤)−1Q_Newton=H^-1G(ZZ )^-1. GD direction. For Q=GQ=G, the numerator is tr(GG⊤)=‖G‖F2tr(GG )=\|G\|_F^2 and the denominator is tr(HG(ZZ⊤/N)G⊤)tr (HG(ZZ /N)G ). Using G=H(W−W⋆)(ZZ⊤/N)G=H(W-W )(ZZ /N), we have tr(GG⊤)=tr(H2(W−W⋆)(ZZ⊤/N)2(W−W⋆)⊤).tr(GG )=tr (H^2(W-W )(ZZ /N)^2(W-W ) ). Similarly, tr(HG(ZZ⊤/N)G⊤)=tr(H3(W−W⋆)(ZZ⊤/N)3(W−W⋆)⊤).tr (HG(ZZ /N)G )=tr (H^3(W-W )(ZZ /N)^3(W-W ) ). Hence the score for the GD direction has the closed form s(G)=tr(H2(W−W⋆)(ZZ⊤/N)2(W−W⋆)⊤)2tr(H3(W−W⋆)(ZZ⊤/N)3(W−W⋆)⊤).s(G)= tr (H^2(W-W )(ZZ /N)^2(W-W ) )^2tr (H^3(W-W )(ZZ /N)^3(W-W ) ). (23) Muon direction. We consider the matrix sign QMuon≔msgn(G)Q_Muon (G). For G=USV⊤G=USV , tr(QMuonG⊤)=tr(UV⊤(VSU⊤))=tr(S)=‖G‖∗=‖H(W−W⋆)(ZZ⊤/N)‖∗,tr(Q_MuonG )=tr (UV (VSU ) )=tr(S)=\|G\|_ = \|H(W-W )(ZZ /N) \|_ , and the denominator becomes tr(HQMuon(ZZ⊤/N)QMuon⊤)=tr(HUV⊤(ZZ⊤/N)VU⊤)=tr(U⊤HU⋅V⊤(ZZ⊤/N)V).tr (HQ_Muon(ZZ /N)Q_Muon )=tr (HUV (ZZ /N)VU )=tr (U HU·V (ZZ /N)V ). Therefore, the Muon-direction score is s(QMuon)=‖H(W−W⋆)(ZZ⊤/N)‖∗2tr(U⊤HU⋅V⊤(ZZ⊤/N)V).s(Q_Muon)= \|H(W-W )(ZZ /N) \|_ ^2tr (U HU·V (ZZ /N)V ). (24) The score for Newton–Muon can be derived similarly. We do not analyze it theoretically here; instead, we evaluate its predicted score numerically and compare it with Muon in our numerical study. Newton direction. The Newton step is QNewton∝H−1G(ZZ⊤)−1Q_Newton H^-1G(ZZ )^-1. Substituting into (19) gives tr(QNewtonG⊤)=tr(H−1G(ZZ⊤)−1G⊤),tr(Q_NewtonG )=tr (H^-1G(ZZ )^-1G ), and, using HQNewton(ZZ⊤/N)=G/NHQ_Newton(ZZ /N)=G/N, tr(HQNewton(ZZ⊤/N)QNewton⊤)=1Ntr(GQNewton⊤)=1Ntr(QNewtonG⊤).tr (HQ_Newton(ZZ /N)Q_Newton )= 1Ntr (GQ_Newton )= 1Ntr(Q_NewtonG ). Therefore, s(QNewton)=Ntr(QNewtonG⊤)s(Q_Newton)=Ntr(Q_NewtonG ). Under G=H(W−W⋆)(ZZ⊤/N)G=H(W-W )(ZZ /N), this can be rewritten as s(QNewton)=tr(H(W−W⋆)(ZZ⊤/N)(W−W⋆)⊤).s(Q_Newton)=tr (H(W-W )(ZZ /N)(W-W ) ). (25) E.2 Isotropic Baseline Numerical Study From (23), (24), and (25), the quadratic score (19) yields explicit expressions for the GD, Muon, and Newton directions. By using different distributions on HH, ZZ, and (W−W⋆)(W-W ), we can evaluate how anisotropy in curvature, activations, and displacement interacts to favor different update geometries, either analytically when expectations simplify or numerically by Monte Carlo simulation. Here we study the simplest baseline in which both the activation and the displacement are isotropic Gaussian; note that the same framework can be reused under more structured, non-isotropic specifications by changing only the distribution for HH, ZZ, and (W−W⋆)(W-W ). Assumptions. First, we sample ZZ with i.i.d. standard normal entries. For the theoretical analysis, we assume the activation second moment is isotropic (while in the simulations we use a finite sample size): ZZ⊤/N=In.ZZ /N=I_n. (26) Second, we assume the square case m=nm=n. Third, we further assume the displacement matrix W−W⋆W-W has i.i.d. standard normal entries. Under (26), G=H(W−W⋆)G=H(W-W ) and the score (19) simplifies to s(Q)=tr(QG⊤)2/tr(HQQ⊤)s(Q)=tr(QG )^2/tr (HQQ ). Let DH≔diag(λ1,…,λm)D_H ( _1,…, _m). We fix this diagonal spectrum (to analyze different levels of curvature anisotropy through the choice of λk\ _k\), then sample a random orthogonal matrix P∈ℝm×mP ^m× m and define H≔PDHP⊤H PD_HP . GD score. Under (26), (23) reduces to s(G)=tr(H2(W−W⋆)(W−W⋆)⊤)2tr(H3(W−W⋆)(W−W⋆)⊤).s(G)= tr (H^2(W-W )(W-W ) )^2tr (H^3(W-W )(W-W ) ). (27) Moreover, for k∈2,3k∈\2,3\, [tr(Hk(W−W⋆)(W−W⋆)⊤)]=tr(Hk[(W−W⋆)(W−W⋆)⊤])=ntr(Hk).E [tr (H^k(W-W )(W-W ) ) ]=tr (H^kE [(W-W )(W-W ) ] )=ntr(H^k). In the large-n regime, the random traces above concentrate around their means, so a standard spectrum-level approximation to (27) is s(G)≈ntr(H2)2tr(H3)=n(∑j=1mλj2)2∑j=1mλj3.s(G)≈ n tr(H^2)^2tr(H^3)=n ( _j=1^m _j^2 )^2 _j=1^m _j^3. (28) Muon score. In the square full-rank case, U∈ℝm×mU ^m× m is orthogonal, so tr(U⊤HU)=tr(H)tr(U HU)=tr(H). Thus s(QMuon)=‖H(W−W⋆)‖∗2tr(U⊤HU)=‖H(W−W⋆)‖∗2tr(H).s(Q_Muon)= \|H(W-W )\|_ ^2tr(U HU)= \|H(W-W )\|_ ^2tr(H). (29) While ‖H(W−W⋆)‖∗2E\|H(W-W )\|_ ^2 does not simplify to a closed form in general, it admits an accurate high-dimensional deterministic approximation that depends only on the spectrum of HH and can be computed numerically. Let S≔1mH(W−W⋆)(W−W⋆)⊤H∈ℝm×m.S 1mH(W-W )(W-W ) H ^m× m. The eigenvalues of SS are the squared singular values of H(W−W⋆)H(W-W ) divided by m. Consequently, ‖H(W−W⋆)‖∗=∑i=1mσi(H(W−W⋆))=m∑i=1mλi(S).\|H(W-W )\|_ = _i=1^m _i (H(W-W ) )= m _i=1^m _i(S). (30) In the regime m=n→∞m=n→∞ with HH deterministic and W−W⋆W-W having i.i.d. standard normal entries, Silverstein–Choi’s analysis of sample-covariance type matrices implies that the empirical spectral distribution of SS converges to a nonrandom limiting law characterized by a fixed-point equation for its Stieltjes transform (Silverstein and Bai, 1995). In particular, writing the eigenvalues of HH as λkk=1m\ _k\_k=1^m, the Stieltjes transform mS(z)m_S(z) of the limiting law is the unique solution in the upper half-plane to the equation mS(z)=−1z⋅1m∑k=1m11+λk2mS(z),z∈ℂ+.m_S(z)=- 1z· 1m _k=1^m 11+ _k^2m_S(z), z _+. (31) Given mS(z)m_S(z), the limiting spectral density ρS(x) _S(x) is obtained by the standard inversion formula ρS(x)=limη↓01πℑmS(x+iη),x>0. _S(x)= _η 0 1π m_S(x+iη), x>0. (32) Equations (31) and (32) suggest a concrete numerical pipeline: 1. Choose a small η>0η>0 and a grid of x values covering the support of the spectrum of SS. 2. For each z=x+iηz=x+iη, solve the fixed-point equation (31). 3. Approximate ρS(x)≈(1/π)ℑmS(x+iη) _S(x)≈(1/π) m_S(x+iη). 4. Compute the limiting mean singular value of H(W−W⋆)H(W-W ) via μ1/2≔∫0∞xρS(x)dx, _1/2 _0^∞ x _S(x)dx, using numerical quadrature on the grid. 5. Using (30), the leading-order scaling is ‖H(W−W⋆)‖∗≈m3/2μ1/2,hence‖H(W−W⋆)‖∗2≈m3μ1/22,\|H(W-W )\|_ ≈ m^3/2 _1/2, \|H(W-W )\|_ ^2≈ m^3 _1/2^2, (33) where the approximation captures the dominant m3m^3 growth and depends only on the spectrum of HH through (31). Substituting (33) into (29) yields a corresponding practical approximation for the expected Muon score, s(QMuon)≈‖H(W−W⋆)‖∗2tr(H)≈m3μ1/22tr(H).Es(Q_Muon)≈ E\|H(W-W )\|_ ^2tr(H)≈ m^3 _1/2^2tr(H). (34) We will use (31)–(34) as the basis for numerically evaluating the Muon numerator (and thus the score) from the eigenvalues of HH. Newton–Muon score. Under (26), Newton–Muon coincides exactly with Muon: QNewton-Muon=msgn(G(ZZ⊤)−1)=msgn(G)=QMuon.Q_Newton -Muon=msgn (G(ZZ )^-1 )=msgn(G)=Q_Muon. Consequently, we use the same approximation s(QNewton-Muon)≈m3μ1/22tr(H).Es(Q_Newton -Muon)≈ m^3 _1/2^2tr(H). (35) Newton score. Applying (26) to (25) yields s(QNewton)=tr(H(W−W⋆)(W−W⋆)⊤).s(Q_Newton)=tr (H(W-W )(W-W ) ). Taking expectation gives the exact identity [s(QNewton)]=tr(H⋅nIm)=ntr(H)=n∑k=1mλk.E [s(Q_Newton) ]=tr (H· nI_m )=ntr(H)=n _k=1^m _k. (36) Comparison. Under the simplifying assumptions (26) and m=nm=n, we can compare the three directions at the level of the theoretical quantities already derived for Newton, GD, Muon, and Newton–Muon. Newton has an exact expected score (36): [s(QNewton)]=ntr(H)=n∑k=1mλk.E [s(Q_Newton) ]=ntr(H)=n _k=1^m _k. For the GD direction, the spectrum-level approximation (28) gives s(G)≈ntr(H2)2tr(H3)=n(∑k=1mλk2)2∑k=1mλk3.s(G)≈ n tr(H^2)^2tr(H^3)=n ( _k=1^m _k^2 )^2 _k=1^m _k^3. For nonnegative eigenvalues λk≥0 _k≥ 0, apply Cauchy inequality to the two sequences λk3/2 _k^3/2 and λk1/2 _k^1/2: (∑kλk2)2∑kλk3≤∑kλk⟹s(G)≲n∑kλk=[s(QNewton)]. ( _k _k^2 )^2 _k _k^3≤ _k _k s(G) n _k _k=E [s(Q_Newton) ]. Equality in Cauchy holds if and only if λk _k is constant across k. Thus, GD matches Newton only when H∝ImH I_m. Newton–Muon and Muon also simplify in this equal-eigenvalues case, and its spectrum-only approximation (34) becomes fully explicit. When H=λImH=λI_m and m=nm=n, the limiting law of S=(1/m)H(W−W⋆)(W−W⋆)⊤HS=(1/m)H(W-W )(W-W ) H is λ2λ^2 times the Marchenko–Pastur law at aspect ratio 11, which yields the closed form μ1/2=83πλ. _1/2= 83πλ. Substituting into (34) gives s(QMuon)≈m3μ1/22tr(H)=m3λ2(83π)2mλ=(649π2)m2λ=(649π2)[s(QNewton)],Es(Q_Muon)≈ m^3 _1/2^2tr(H)= m^3λ^2 ( 83π )^2mλ= ( 649π^2 )m^2λ= ( 649π^2 )E [s(Q_Newton) ], so the isotropic asymptotic approximation predicts that Muon achieves a constant fraction of the Newton score (about 0.720.72), while GD coincides with Newton. Beyond H=λImH=λI_m, the Muon and Newton–Muon approximation (34) and (35) still depend only on the spectrum of HH through μ1/2 _1/2 (computed from (31) and (32)), but there is no longer a simple closed-form comparison between the Muon/Newton–Muon score (34) and the GD approximation (28), as we will show next in our numerical study. Under the isotropic-input theory (26), the Newton–Muon and Muon theoretical predictions coincide exactly; any empirical separation between them comes solely from finite-sample deviations of ZZ⊤/NZZ /N from InI_n. In all experiments below we fix λmin=10−4 _ =10^-4 and vary only (N,p)(N,p). For each setting, we run 10241024 independent simulations and report the mean score together with the 2.5%–97.5% interval. The theory marker in the plot denotes the theoretical score in the limit of infinite data and infinite dimensions m,n→∞m,n→∞. Baseline (N=8192N=8192, p=0.3p=0.3). Figure 8 shows a clearly anisotropic spectrum of HH. In this ill-conditioned regime, both Muon and Newton–Muon yield substantial score improvements over the raw GD direction, and Newton–Muon is slightly better than Muon. Figure 8: Baseline configuration (N=8192N=8192, p=0.3p=0.3): spectrum of HH (left) and mean absolute scores s(Q)s(Q) for GD, Muon, Newton–Muon, and Newton (right). Uniform curvature (N=8192N=8192, p=2.4p=2.4). Figure 9 corresponds to a spectrum that is nearly flat at the top and drops sharply only near the tail. In this more uniform-curvature setting, the score gaps among GD, Muon, and Newton–Muon narrow, indicating that Newton–Muon and Muon help when HH is strongly anisotropic. Here, Newton–Muon is still slightly better than Muon. Figure 9: More top-uniform curvature (N=8192N=8192, p=2.4p=2.4): spectrum of HH (left) and mean absolute scores (right). Smaller data (N=1024N=1024, p=0.3p=0.3). Figure 10 keeps the same curvature shape as the baseline but uses a smaller sample size N to form ZZ⊤/NZZ /N, making the empirical second moment noisier. Relative to Figure 8, the gap between Newton–Muon and Muon becomes much larger, suggesting that Newton–Muon more effectively compensates for activation anisotropy when the sample size is small. Figure 10: Smaller-N data (N=1024N=1024, p=0.3p=0.3): spectrum of HH (left) and mean absolute scores (right). Conclusion. Across these three runs, Muon and Newton–Muon achieve their largest score gains when curvature is strongly anisotropic (Figure 8); when HH becomes more uniform, the benefit shrinks or even vanishes (Figure 9). Newton–Muon outperforms Muon in these three experiments, and especially when N is reduced so that ZZ⊤/NZZ /N is less isotropic, the advantage of Newton–Muon over Muon further increases (Figure 10). Appendix F Non-Isotropic Assumption Here we report attempts to replace the isotropic assumption ΣW∝Im _W I_m by estimating ΣW _W from training dynamics. (5) motivates the use of a non-isotropic ΣW _W. When ΣW _W admits a factorization of the form ΣW=MM⊤ _W=M , the update in (5) can be implemented without explicitly constructing ΣW1/2 _W^1/2. Proposition 3 (Factorized form of (5)). Suppose ΣW=MM⊤ _W=MM for some M∈ℝm×rM ^m× r. Then Q⋆=Mmsgn(M⊤G(ZZ⊤)−1).Q =Mmsgn (M G(ZZ )^-1 ). Proof of Proposition 3. Let k≔rank(M)k (M) and take a compact SVD M=UMSMVM⊤M=U_MS_MV_M , where UM∈ℝm×kU_M ^m× k and VM∈ℝr×kV_M ^r× k have orthonormal columns and SM∈ℝk×kS_M ^k× k is diagonal with positive entries. Then ΣW=MM⊤=UMSM2UM⊤ _W=M =U_MS_M^2U_M , so ΣW1/2=UMSMUM⊤ _W^1/2=U_MS_MU_M . Hence ΣW1/2G(ZZ⊤)−1=UMSMUM⊤G(ZZ⊤)−1. _W^1/2G(Z )^-1=U_MS_MU_M G(Z )^-1. Using the compact-SVD definition of the rectangular matrix sign, if UMU_M has orthonormal columns, then msgn(UMX)=UMmsgn(X)msgn(U_MX)=U_Mmsgn(X) for any matrix X. Applying this with X=SMUM⊤G(ZZ⊤)−1X=S_MU_M G(Z )^-1 gives msgn(ΣW1/2G(ZZ⊤)−1)=UMmsgn(SMUM⊤G(ZZ⊤)−1).msgn ( _W^1/2G(Z )^-1 )=U_Mmsgn (S_MU_M G(Z )^-1 ). Substituting into (5), we obtain Q⋆=ΣW1/2msgn(ΣW1/2G(ZZ⊤)−1)=UMSMmsgn(SMUM⊤G(ZZ⊤)−1).Q = _W^1/2msgn ( _W^1/2G(Z )^-1 )=U_MS_Mmsgn (S_MU_M G(Z )^-1 ). (37) On the other hand, M⊤G(ZZ⊤)−1=VMSMUM⊤G(ZZ⊤)−1M G(Z )^-1=V_MS_MU_M G(Z )^-1. Applying the same identity with VMV_M yields msgn(M⊤G(ZZ⊤)−1)=VMmsgn(SMUM⊤G(ZZ⊤)−1).msgn (M G(Z )^-1 )=V_Mmsgn (S_MU_M G(Z )^-1 ). Therefore Mmsgn(M⊤G(ZZ⊤)−1) Mmsgn (M G(Z )^-1 ) =(UMSMVM⊤)(VMmsgn(SMUM⊤G(ZZ⊤)−1)) =(U_MS_MV_M ) (V_Mmsgn (S_MU_M G(Z )^-1 ) ) =UMSMmsgn(SMUM⊤G(ZZ⊤)−1), =U_MS_Mmsgn (S_MU_M G(Z )^-1 ), which matches (37). ∎ Motivated by this identity, we explored practical estimators of MM in benchmark training runs, including momentum-based and diagonal variants derived from update statistics. However, these non-isotropic variants did not yield consistent improvements over Newton–Muon. Keeping a momentum buffer of the actual update. We maintain a momentum buffer of the actual parameter updates, denoted by MM. Heuristically, if the update directions tend to point toward an optimum W⋆W , then this buffer can be viewed as a coarse proxy of W−W⋆W-W . In practice, we formed two candidate directions QI=msgn(G(ZZ⊤)−1),QM=Mmsgn(M⊤G(ZZ⊤)−1),Q_I=msgn (G(ZZ )^-1 ), Q_M=Mmsgn (M G(ZZ )^-1 ), and then combined them by a convex average after normalizing their scales Q=(1−λ)QI‖QI‖F+λQM‖QM‖F,λ∈[0,1].Q=(1-λ) Q_I\|Q_I\|_F+λ Q_M\|Q_M\|_F, λ∈[0,1]. Despite the motivation from Proposition 3, this averaging did not yield consistent gains, and typically underperformed the Newton–Muon update QIQ_I. Using diagonal estimation of the covariance from update statistics. We also tried a computationally cheaper approximation that estimates a diagonal ΣW _W from the second moment of the actual update directions. Concretely, for each matrix parameter W∈ℝm×nW ^m× n we maintain a state vector ∈ℝm u ^m that tracks the row-wise magnitudes of the applied Muon directions via an EWMA: ←β+(1−β)‖Q‖row, u←β u+(1-β) \|Q \|_row, where ∥⋅∥row\|·\|_row denotes the ℓ2 _2 norm along the column dimension. We then form a damped diagonal estimate of ΣW _W using the squared row magnitudes, ⊙2≈diag(ΣW),M2=diag(λ⊙2+(1−λ)⊙2¯ 1m), u 2\ ≈\ diag( _W), M^2=diag (λ u 2+(1-λ) u 2\,1_m ), where ⊙2¯ u 2 is the mean of the entries of ⊙2 u 2, λ∈[0,1]λ∈[0,1] is a damping coefficient, and m∈ℝm1_m ^m is the all-ones vector. Let MM be the positive diagonal square root of M2M^2. The resulting update takes the form QM∝Mmsgn(MG(ZZ⊤)−1).Q_M Mmsgn (MG(ZZ )^-1 ). Despite being inexpensive, this diagonal ΣW _W variant did not outperform Newton–Muon in our experiments and in some cases performed worse than standard Muon. Conclusion. Overall, if the update direction is systematically misaligned with W−W⋆W-W , then the resulting estimate of ΣW _W inherits this bias, and since the next update direction depends on this estimate, the bias can reinforce itself, leading to slower training. By using an isotropic proxy (i.e., ΣW∝Im _W I_m), we avoid imposing assumptions about the unknown distribution of displacements while still capturing the most robust geometric component of the update through the matrix sign. This helps explain both the computational simplicity and the stability of Newton–Muon.