The Polar Express: Optimal Matrix Sign Methods and their Application to the Muon Algorithm
Authors: Noah Amsel
David Persson
Christopher Musco
Robert M. Gower
Venue: ICLR 2026 Oral (2026)
Keywords: polar decomposition
matrix sign
numerical linear algebra
muon
optimization
approximation theory
Primary area: optimization
TL;DR: We introduce a GPU-friendly algorithm for computing the polar decomposition of a matrix to low accuracy that is optimal in its class. This improves Muon.

Abstract: Computing the polar decomposition and the related matrix sign function has been a well-studied problem in numerical analysis for decades. Recently, it has emerged as an important subroutine  within the Muon algorithm for training deep neural networks.
However, the requirements of this application differ sharply from classical settings: deep learning demands GPU-friendly algorithms that prioritize high throughput over high precision. We introduce *Polar Express*, a new method for computing the polar decomposition. Like Newton–Schulz and other classical polynomial methods, our approach uses only matrix-matrix multiplications, making it 
very efficient on GPUs.
Inspired by earlier work of Chen \& Chow and Nakatsukasa \& Freund, *Polar Express* adapts the update rule at each iteration by solving a minimax optimization problem.
We prove that this strategy minimizes error in a worst-case sense, allowing *Polar Express* to converge as rapidly as possible both in the early iterations and asymptotically.
We also address finite-precision issues, making it practical to use in `bfloat16`. When integrated into Muon, our method yields consistent improvements in validation loss for a GPT-2 model on one to ten billion tokens from the FineWeb dataset, outperforming recent alternatives across a range of learning rates.

Published as a conference paper at ICLR 2026
THE POLAR EXPRESS: OPTIMAL MATRIX SIGN METHODS AND THEIR APPLICATION TO THE MUON ALGORITHM
Noah Amsel
David Persson
Christopher Musco New York University emusco@nyu. edu
New York University noah. wwwellnyu.eeee.mmommmm edu
New York University
Flatiron Institute dup21 0@nyu edu
Robert M. Gower
Flatiron Institute
rgover@flatiroonnssittee
ABSTRACT
Computing the polar decomposition and the related matrix sign function has been well-studied problem numerical analysis for decades. Recently, it has emerged as an important subroutine within the Muon optimizer for training deep neural networks. However, the requirements of this application differ sharply from classical settings: deep learning demands GPU-friendly algorithms that prioritize high throughput over high precision. We introduce Polar Express, new method for computing the polar decomposition. Like Newton-Schulz and other classical polynomial methods, our approach uses only matrix-matrix multiplica- tions, making it very efficient on GPUs. Inspired by earlier work of Chen & Chow and Nakatsukasa & Freund, Polar Express s adapts the update rule at each it- eration by solving minimax optimization problem. We prove that this strategy minimizes error in a worst-case sense, allowing Pol ar Express to converge rapidly possible both in the early iterations and asymptotically. We also address fimite- precision issues, making it practical to use in b bfloat16. When integrated into Muon, our method yields consistent improvements loss for GPT-2 model trained on one to ten billion tokens from the FineWeh dataset, outperforming recent alternatives across a range of learning rates.
finite-precision
validation
INTRODUCTION
Advanced linear algebra is making its way into deep learning. Efficient algorithms for computing matrix functions have found exciting new applications in training neural networks. particular, ap- proximations to the matrix-invers are used in the full Adagrad method (Duchi et al., 201 1), the ma- trix square-root and quarter-root appear as subroutines the Shampoo and Soap optimizers (Gupta et al., 2018; Shi et al., 2023; Vyas et al., 2025), and most recently, the matrix sign function has become a key ingredient the Muon optimizer (Bernstein & Newhouse, 2024b;a; Jordan et al., 2024b). While the problem of computing these matrix functions has been studied by numerical analysts for decades, applications in deep learning come with different requirements than those in computational science. For deep learning, it is critical take maximum advantage of GPU-friendly operations like matrix-matrix products and to avoid less parallel operations. Moreover, memory overhead must be small handle large models. On the other hand, high accuracy is typically less important; the gold standard of sixteen digits of accuracy is overkill in deep learning.
Given these considerations, there is need to develop new matrix function methods that are tailor- made for deep learning applications. We take on this challenge by designing state-of-the-as GPU-friendly algorithm for computing the matrix sign function, or more generally, for computing
'https://gg thub. com NoahAmseh /Pol |Pclarxxpres. arExpress

---

Published as a conference paper at ICLR 2026
the polar decomposition of rectangular matrix. We apply (Algorithm 1, Implementation 1) to compute the descent direction in the increasingly popular Muon optimizer. In Figure 1, show that using Polar Express within Muon consistently results in lower validation loss across all learning rates when training a GPT-2 model, as compared to other matrix sign methods (Cesista et al., 2025; Jordan al., 2024h).
Polar Express method
4.4
3.46
Tupe-PdazzZSy momnshhilln
4.2-
3.44
Fanl Validalion
4.0
3.40
0.3
3.38
3.6
3.36
3.4
3.31
Lesaning Rate
Epoch
Figure 1: Training GPT-2-Large model (774M params) 1 billion tokens from the FineWeb dataset (Aroca-Quelle et al., 2023). The label muon-<names refers to implementing Muon using <name> to compute the polar factor. Left: final validation loss across learning rates. Right: valida- tion loss across epochs using the best learning rate. The hest learning rate (Ir) and final validation loss for each method were muon You (Ir 0.02): 3.399, muon Jordan. (Ir 0.02): 3.398 and muon Polaollaatee Pxpp (Ir = 0.02): 3.340.
1.1 THE MUON METHOD
The Muon optimizer has recently gained popularity for training large language models, often out- performing state-of-the-art adaptive gradient methods like Adam and AdamW (Kingma & Ba, 2015; Loshchilov & Hutter, 2019). Muon has been used to set records for the NanoGPT speedrun (Jordan et al., 2024b), expand the Pareto frontier of performance guage models (Liu et al., 2025; Shah et al., 2025), and even to train a trillion parameter frontier L.LM (Kimi Team al., 2025).
training FLOPs for large lan-
The Muon update rule (Bernstein & Newhouse, 2024b) is defined as follows. Let X, B >> 0 be the learning rate and momentum coefficient hyperparameters. (By default, 0.9.) Let W E mxxx be the weight matrix of a given neural network layer at iteration f, and let G; E mmxx be its (stochastic) gradient. Let M, Ruxn be the running momentum estimate of the gradient, where Mo = O. The Muon update is given by
M, = 3M, + (1 B)G
W = Wr Aolaar(M.)
Whereas standard stochastic gradient descent (SGD) with momentum updates the weight matrix by taking step in the direction MM, the Muon method steps in the direction polar MM), where polar (M)) denotes the closest semi-orthogonal matrix to M (Higham, 2008, Chapter 8). Con- cretely, if M UEVT T is the singular value decomposition (SVD) of M, then
(1) The matrix polar (M) can be seen as generalization the matrix sign function to rectangular ma- trices (Benzi& & Huang, 2019). Indeed, when M is square symmetric with cigendecomposition M = VAVT. polar(M) exactly coincides with the matrix sign function sign(M) V sign A V (Higham, 2008, Chapter 5). Equivalently, polar(M M) is the left orthogonal factor of the polar de- composition M (Higham, 2008, Chapter 8). The motivation for Muon that polar (M) gives the steepest-descent direction with respect to the spectral norm (instead of the Frobenius norm, as in standard SGD). For analysis and further discussion on Muon we refer the reader to (Jordan et al., 2024b; Bernstein & Newhouse, 2024b; Pethick al., 2025; Riabinin al., 2025; Carlson et al., 2015a;b). this paper, we take the Muon update rule as given and focus on the problem of efficiently computing the polar decomposition polar(M).
polar (M) == UVT
2

---

Published as a conference paper at ICLR 2026
1.2 COMPUTING THE POLAR FACTOR
Although polar(M) can be computed directly via an SVD in O min(m, n3) time, doing so is prohibitively expensive deep learning applications, especially as standard SVD algorithms fail to take full advantage of the parallelism available on GPUs. There has been significant work on highly- parallel methods for the SVD, but the most common approaches actually require computing the matrix-sign function as a subroutine (Nakatsukasa & Freund, 2016; Nakatsukasa & Higham, 2013). Numerical analysts have spent decades developing iterative methods for computing polar (M). This rich line of work includes Newton-Schulz (Higham, 2008, Chapter 8), Padé iteration (Kenney & Laub, 1991; Higham, 1986), the Newton and scaled Newton iterations (Higham, 2008, Chapter 8), the QDWH iteration (Nakatsukass al., 2010; Nakatsukasa & Higham, 2013), and Zolo-pd (Nakat- sukasa & Freund, 2016). Unfortunately, as discussed in Appendix B, most of these methods are based rational approximations to the function sign(æ) and require computing matrix inverses or QR decompositions. Such methods are ill-suited to GPU acceleration and deep learning applica- tions. In contrast, the older Newton-Schulz method is based polynomial approximation of sign(a:) and uses only matrix-matrix products. Thus, Muon initially used Newton-Schulz (Bernstein & New- house, 2024a). Indeed, Muon stands for "MomentUm Orthogonalized by Newton-Schulz" (Jordan et al., 2024b). For a more comprehensive discussion on prior work, see Appendix B.
The Newton-Schulz methods. Newton-Schulz constructs a sequence of approximations X, P polar(M) follows:
3 X, X X-- 2
1
X M MF.
x X,
At each iteration, this rule effectively applies the cubic polynomial p(x) = d 1 x to each sin- gular value of X,. The scalar fixed-point iteration #4+1 = p(xx) converges to sign( (z)) as 0o, provided |oo| < 1. As a result, the matrix iteration satisfies lim X, = UV (-00
2
polar Xo)- Higher-degree versions of Newton-Schulz follow the same principle. For example, the degree-5 polynomial p(z) = (15x- - 10x 3//)/8 converges even faster. The Newton-Schulz iterations converge super-exponentiall when X, is sufficiently close to polar(M) (M), but they suffer from slow initial convergence; when Xo is far from polar(MM M,), the approximation improves slowly over the first few iterations. Due to the slow initial convergence of Newton-Schulz, Chen & Chow (2014) 4) developed version of the Newton-Schulz iteration, which adapts the polynomial at each iteration. The resulting method achieves a faster initial convergence, while retaining super-exponential con- vergence in later iterations. Polar Express is inspired by their method.
The Jordan and You methods. In Muon, high accuracy approximations to polar( ) are usually not necessary. The primary goal is instead to compute a coarse approximation in as few iterations as possible. To accelerate convergence in the low-accuracy regime, Jordan recently proposed fixed-point iteration based on the polynomial p(z) = 3.4445x - 4.7750x3 + 2.031 5x which was found using heuristic numerical search (Jordan et al., 2024b). Unlike Newton-Schule the scheme that Jordan proposed does not converge to (M), but plateaus at an error of ~ 0.3. However, it reaches this level of accuracy rapidly and outperforms the Newton-Schulz when only small number of iterations are performed. Building on this idea, You proposed method that applies six different polynomial updates in succession, which were again found by heuristic search. This method achieves better accuracy than Jordan's 's but still fails to converge (Cesista et al., 2025).
1.3 CONTRIBUTIONS
We present Polar Express (Algorithm 1), an iterative method for approximating polar(M Our method dynamically adapts the polynomial update rule at each iteration, prioritizing rapid progress in the initial stage and high accuracy in the later stage. Polar Express constructs polynomials PI....,P7 PT so that the resulting composition is the optimal approximation to the sign function with respect to the supremum (L*) norm (Theorem 3.1). By iteratively applying these polynomials to M, Polar Express computes an approximation to polar(M) that is optimal in the warst-case. Our method converges to polar( M) super-exponentially (Theorem 3.3), and it quickly reaches a good approximation within just five to ten iterations. This carly-stage acceleration especially valuable in deep learning applications, where runtime efficiency takes precedence over
3

---

Published as a conference paper at ICLR 2026
high accuracy. In contrast, classical methods like Newton-Schull suffer from a slow initial conver- gence, while recent heuristic proposals (Jordan et al., 2024b; Cesista et al., 2025) fail to converge. Our method is efficient run GPUs, using only few matrix- matrix products per iteration. We give an explicit instantiation of Polar Express in Implementation 1, which incorporates minor modifications to make it compatible with half-precision arithmetic (see Section 3.4). Implementa- tion is very short and easy to use, with no dependencies except PyTorch. It serves as a drop-in replacement for previous methods. In numerical experiments, Polar Express outperforms pre- vious methods on synthetic matrices and gradient matrices from GPT-2 transformer (Figure 3). We demonstrate the effectiveness of using Polar Express within the Muon optimizer in Figure I, showing that it consistently improves the training of GPT-2 language models on billion tokens of the Web dataset (Aroca-Ouellette et al., 2023). Our method has been adopted into the NanoGPT speedrun (Jordan et al., 2024a), heavily optimized implementation that serves as a benchmark for LLM training efficiency.
Notation. We let || and || denote the Frobenius norm and spectral norm (largest singular value) of matrix M, respectively. We denote the spectrum (set of singular values) by o ((M). Let Pa be the set of polynomials of degree at most d. For odd d, Podd denotes the set of polynomials of degree at most d containing only odd-degree monomials. For a polynomial p. deg(p) is its degree. Let sign(x) be the scalar sign function, which satisfies sign(0)
1 if > 0 and sig)(z) = 1 if x < O. For a polynomial E podd and a matrix M with rank reduced SVD given by M = UEVT T and positive singular values 01 2 ... 2 (aank(M) > 0, we define p(M) == Up(E)VT, where p( [)) is the diagonal matrix with diagonal entries p(o:)) for i = 1.....,ank((MM
O, sign(x)
2 APPROXIMATIONS BY COMPOSITIONS OF POLYNOMIALS
To design a GPU-friendly method for computing polar(M), we limit ourselves the following GPU-friendly operations: (i) linear combinations of matrices (given scalars B, 7 E R and matri- ces B and C, compute 3B + 7C) and (ii) matrix-matrix products (compute BC). While both these computational primitives are well-suited for parallel computing environments, matrix-matrix products come at a higher computational cost than linear combinations. Therefore, our method attempts to minimize the number of matrix-matrix products. A key observation is that we can com- pute odd monomials of M = U VVV using the following formula: M2 == U M (M) M) Hence, for an odd polynomial p(z) = agæ + +. + + 22 we can compute
p(M) M + ( M) +.... a. MMM M)
It has been shown that for an arbitrary polynomial one requires ((de(p))) products to compute M(M) (Paterson & Stockmeyer, 1973); see also Jarlebring & Lorentzon (2025) for related work. This compares favorably to the naive approach that forms all monomials in and then sums them together, which requires (dde(pp) products. However, if p can be expressed as composition of T polynomials, each of degree d
== Pr o pr-1 ..... o Pl: then the degree of p is d", and p(M) can be efficiently computed recursively by
(3)
Xo = M, X, = (x., for t = 1,2.....7?
(4)
The final iterate is XT p(M), which we compute with just O(Td) matrix-matri products. Iterative methods for polar( can be seen in this light, For instance, the degree-5 Newton-Schulz I,.... .T The Pı approximates sign(z), and the approximation error goes to 0 as T grows. In this paper, we ask the following question: what choice of pr o P1 gives the best approximation to sign(x))
method uses the polynomial update pz(z) 15 10 + -52² for each composition p Pr
The method we will present is optimal in the following sense: given lower and upper bounds E and u on the singular values of M, an odd degree d E N, and the number of iterations T E N, our method computes the composition p" (M) that minimizes the worst-case error in the spectral norm. That is,
arg min pppropy-1 11 ops P.end
max [|polar(M) p(M) ME (M)
(5)
For non-symmetric matrices, e.g. rectangular matrices, cannot compute even polynomials of the singu- lar values without first explicitly computing the SVD. We therefore restricted to odd polynomials.
4

---

Published as a conference paper at ICLR 2026
22
6
8g
10g
Figure 2: The evolution of the first three optimal polynomials P1, P2, and p3 and the corresponding lower bounds l++1 (f,) and upper bounds W,+1 = 2 le+1, as described in Theorem 3.1. The horizontal black line shows y The polynomial degree
We set = 0.03 and #y 1.
Given that polar(M) p(M) U(I P(E))VT, and by the unitary invariance of the spectral norm, we have that (5) is equivalent to
pz
arg min
max
- --p(())]
(6)
pppropr 117-16600 16..pp1 2E P
In other words, the problem given in (5) reduces to that of finding a "uniform" approximation to the constant function + over the interval [E, u]. as given in (6). Uniform approximation on an interval by polynomials or rational functions of given degree theory (Trefethen, 2020). Here, seek approximation of particular form-6 composition of odd polynomials of fixed degrees. In the next section, we solve the optimization problem of (6) and the solution to create Polar Express.
central topic approximation
THE POLAR EXPRESS
3.1 GREEDY IS OPTIMAL
The key observation is that the polynomial used in each iteration can be chosen greedily, given the choice of polynomials from the previous iterations. For the first iteration, choose Pl so map the interval [é, u] as close to as possible. That is, it minimizes maxze](.a] [1 P1 (x)). The image of pı will be a new interval [22, w2]. where
f2 = min (z) 2E(,.u)
U22 = max pi (i())
=E((.u)
now pick P2 to map the interval [62, u2] as close to as possible, obtaining a new interval [63, w3]
that is the image of [e, through P2 P1- We continue this process for many iterations desired.
The following theorem guarantees that this process finds the solution to (6), and thereby also (5). The scheme is also outlined in Figure 2. which demonstrates the evolution of the lower bounds &,, the upper bounds we, and the polynomials pt across iterations. The proof is in Appendix C.
Theorem 3.1. Let be odd and define & = and u1 = u. For = 1,..., define
P: = arg min max - (-(())) l = min pr(x), Ue+1 = max Pc(x) pepodd ae(e..)
(8)
sE(f..4
=E(P...))
The resulting composition p= := pr o PT o ..... o pı is optimal and the error is given by:
max p" (xxx
min
max |1 p(x) 1 Crtr.
(9)
DE [E(
popropt por xERe,u] Prepodd
For completeness, the equivalence between (5) and (6) proven in Appendix E.

---

Published as a conference paper at ICLR 2026
Furthermore the new error, lower and upper bounds can be computed through
l 4-1 = pp(ee),
u++1 = 2 - c++1: and
max [e(
|| pt(x)) = =1-4 4411
(10)
Remark 3.2 (Why a fixed degree?). We note that choice of the degree of each pr, .2...... .r need not be the same for Theorem 3.1 to hold. More generally, one may specify sequence of degrees di....,dy dr and define each Pt as Pt. arg min pePa, maxxz](e, [p(x)) 1| for 1 = 1.... T.
However, Lee et al. (2022, Table 2) supports setting de = 5, as we do.
Fortunately, (10) shows that once Pt has been found, we can compute the new lower and upper bounds 4++1 and u4+1 simply by evaluating ps (f,). Hence, for any fixed upper and lower bounds on the singular values of M, we can precompute all the polynomials P1..... PT and the bounds [&,, ui]... ..l...... [ uT-1]- Then, applying the iterative procedure of (4), the final iterate Xr will satisfy the following error bound:
polar(M) - Xrll2 polar(M) - (M -1- 6T+1.
(11)
From the optimality guarantee of Theorem 3.1, know that our method converges least fast as the Newton-Schuil iteration of the same degree. Combining this fact with an existing analysis of Newton-Schulz, we immediately get the following convergence guarantee showing that our method enjoys faster than exponential convergence. The proof can be found in Appendix D.
Theorem 3.3. Let M be matrix normalized so that o (M) C [ 1]. Let X p" (M), where p" is the polynomial from Theorem 3.1 with 2q + 1. Then, we have
polar(M) Xr||2 < | e (1(())))
(12)
Hence, for d = 3 and d = 5 the method converges quadratically and cubically, respectively.
In fact, our method is strictly faster than Newton-Schulz, even if m (M) < l. When min = t Polar Express is about twice as fast as Newton-Schulz (cf. Chen & Chow (2014, Section 3.1)). Recent work has analyzed the stability and convergence of Muon when the polar factor is computed inexactly (Shulgin et al., 2025; Refael et al., 2025). Combining these analyses with Theorem 3.3 immediately yields convergence guarantee for Muon implemented with Polar Express.
3.2 FINDING THE OPTIMAL POLYNOMIAL FOR EACH ITERATION
Theorem 3.1 shows that we can solve (6) by greedily choosing the optimal approximation E Padd for each interval [&,, ur] for t I,... T. In this section, show how find each pt. Since are now focused on just one iteration, we drop the subscripts. Given & and u, we wish to solve the following optimization problem:
arg min max [1 - p(z)
(13)
pEpadd zE|f.a)
That is, we seek a minimax or uniform approximation of the function z + on [e, u] from the set of odd polynomials. (Equivalently, we seek minimax optimal approximation to sign(x: (r) on [-u, E] U [e, u).) Problems of this form are well-studied approximation theory and numerical analysis. The key mathematical insight underlying their solution is the Equioscillation Theorem, which we state formally for our setting in Lemma C.I. This theorem the basis of the Remez algorithm (Pachón & Trefethen, 2009; Parks & McClellan, 1972), general-purpor finds a (nearly) optimal polynomial approximation of a given degree to any function on any interval. With very minor modification to handle the constraint that p be odd, Remez can solve (13).
method that
However, the Remez algorithm is complicated and notoriously difficult to implement correctly.* Fortunately, do not need the algorithm in its full generality; seek only low-degree polynomial approximations, and the function we wish to approximate is just f(x) = 1. We use the Equioseilla- tion Theorem to derive (17), an explicit, closed-form solution (13) for the degree
case. Up
Foor implementations of the general Remez algorithm, we recommend Chehfun or o remez.

---

Published as a conference paper at ICLR 2026
to rescaling, this turns out be the same polynomial derived by different means in Chen & Chow (2014). For d = 5, we present Algorithm 2, simpler way of solving (13) that is mathematically equivalent to Remez in our setting. This algorithm is implemented in its entirety in Implementa- tion 2. For more details, we refer the reader to Appendix F.
3.3 UPPER AND LOWER BOUNDS ON THE SINGULAR VALUES
To instantiate our method, we need upper and lower bounds and on the singular values of the input matrix M. A trivial upper bound is given by M IE.. This can be quite loose in the worst case. In practice, it is off only by small constant factor because the gradient matrices of the weights of dense linear layers in neural networks tend to have small effective rank (Yang et al., 2024). We therefore rescale M by || and set 1. It is difficult efficiently find a good lower bound on minn so we are forced to guess. Fortunately, the consequences of a bad guess are not severe. The method converges for any E (0, u]. and even an order of magnitude error only delays convergence by few iterations. For matrices stored in floating point arithmetic, the singular values are usually larger than machine precision Cmach (Boutsikas et al., 2024). We work bfloar fflttatl 6, which has Gmach = 2-8 z 3.91 10-3 so we set l = 10-3 Since we use these bounds for all input matrices, we can pre-compute the optimal polynomials once and apply them to as many inputs as we want.
3.4 FINITE PRECISION CONSIDERATIONS
When working in unite-precision arithmetic, especially the half- precisrrciionn bfloatl 6 format used in deep learning, we must take some care avoid blowups and other problems due to numerical error. To this end, we make a few small but crucial changes to the method in the offline stage that stabilize it with a negligible effect on aceuracy. One issue arises when numerical round-off creates singular values that are slightly larger than our current upper bound ut. To fix it, we replace each polynomial Pt by I + (r(1),.0)), effectively increasing ut- Another issue, identified by Nakatsukasa & Higham (2013), is due to the non non-monotonicity of pe. We address it by using slightly suboptimal (but less oscillatory) polynomials in the early iterations, as suggested by Chen & Chow (2014). For detailed discussion on the finite precision considerations, refer to Appendix G.
3.5
THE ALGORITHM
We give the pseudocode of our proposed method for any degree in Algorithm 1. We give the specific Python code of the Pola ar Express with degree d 5 and = 10-3 used in our GPT experi- ments in Implementations and in Ap- pendix A. Both incorporate the finite pre- cision considerations discussed in Sec- tion 3.4. Our algorithm precomputes the polynomials PI.....PP of Theorem 3.1 in full precision using the results of Sec- tion 3.2 (or the Remez algorithm for d > 5). This stage is offline because the coef- ficients the polynomials are only com- puted and stored once. For every subse- quent call to the algorithm, these coef- ficients are reused and the offline stage is skipped. For instance, in Implemen- tation 1 these polynomials have been precomputed and stored in the variable coeffs list.
Algorithm 1 The General Polar Express
input: Matrix M, iteration count T, degree d, approxi- mate lower bound l. output: An approximation Xy to polar(M).
Offline: precompute poly nominls in float 64
&1 = L, u1 = fort = 1,2, 1,2....,T do
4
Solve using Remez (Appendix F):
Pt = arg min pEpadd x =E(max(fe,3u,/ (1, /10), uil 11 p(x)
max
Pe - Pc(111011)
l+ < (eo), Me+1 t 2 2-644 l171
end for
Online: apply precomputed poly polynomials in beloat16
10
Set Xa M/ [[| II: + 10
11
for = 2, do
12
X. (X,,))
13
end for
14
return XT.
The online stage can be performed in lower precision (bfloat16) for greater speed on a GPU. Horner's rule can be used to carry out each iteration. For instance, if pt had + cx then
7

---

Published as a conference paper at ICLR 2026
Synthetic Matrix
GPT-2 Gradient
min maax 10-0
minn/max
Layer mlp.c proj
100
1.00
Spectral Error
0.75
10-1
0.50
10-2
0.25
0.00
10
5
10
15
20
25
2
4
6
8
10
Iteration
Iteration
Newton-Schulz (d 55
You
PolarExp, 6= 10-6
Jordan
PolarExp, 8 10
PolarExp, 8= 10
Figure 3: Convergence of degree-5 polynomial methods. Polar Express outperforms other methods at every iteration when tuned properly. Left panel: synthetic matrix with max 1, min 10 Right panel: gradient from mandoml--nnitiaiize GPT-2 model on a batch of language modeling data. Shaded region shows 90% interval over 512 batches of data.
X (al + Y,- [9] + cY1-1)) where Y X, X, A simple implementation of the offline stage of Algorithm is given Implementation 2. For deep learning applications, we recommend using d = 5 and 7' = 5 or with L = 10-3 With these parameters, the offline Implementation
gives the polynomials encoded in coef fs in Implementation 1. All told, our proposal for Muon is to apply the composition of these polynomials to M/[[M|]] F + 105
stage as implemented
NUMERICAL EXPERIMENTS
4.1 CONVERGENCE OF POLAR EXPRESS
We compare Polar Express against degree-5 Newton-Schulz and the methods of Jordan et al. (2024b) and Cesista et al. (2025). We first generate random matrix whose singular values are evenly spaced on logarithmic scale between 10-6 and 1, with singular vectors chosen randomly. The left panel of Figure shows the results. Since all the methods in this plot use degree-5 poly- nomials, their computational and runtime costs are proportional to the number iterations. As expected, Newton-Schulz converges but makes almost no progress for the first 17 iterations. Jor- dan's method rapidly achieves an error of 0.3 after just iterations, but ceases to converge fur- ther. You's method, which is only defined for six iterations, converges at similar rate Jordan's method. When Polar Express instantiated with f = 7mins it dominates the other methods every iteration, achieving excellent accuracy after just 11 iterations and converging about twice as fast as Newton-Schulz to any given error. Even when É is wrong by two orders of magnitude in either direction, the method remains competitive, though it does outperform Jordan's 's method until iteration 13 or 14. We also test convergence on non-synthetic matrix: the gradient of a weight matrix from the fourth transformer block of GPT-2 model (Figure right). Again, the version of Polar ar Express outperforms the other methods, but setting & to be many orders of magnitude too small can delay convergence. Note that Figure measures error in the spectral norm. For many applications may be satisfied with looser measure of error; Appendix H.I.
'In Appendices I and J, describe two further algorithmic ideas. They not used in our Muon experi- ments but they may be heneficial in other settings, and helieve they merit further study.
8

---

Published as a conference paper at ICLR 2026
mau-Jonnin wou--You You ailonr Bum--Pllaeap
4.2
Validation Loos
Valaiation
4.0
3.8
3.8
10"3
10
0.0
0.2
0.4
0.6
0.8
Learning Rate
Epoch
Figure 4: Training GPT-2-Small (124M) model on Billion tokens of the Fine data set (Aroca- Ouellette et al., 2023). muon-<met hhod> denotes Muon with 5 iterations of <method> to compute polar(M) No weight decay used. Left: final validation loss vs. learning rate. The best final val- idation losses for each method were adamw(Ir =0.0005): 4.197, muon- Jorddnnn =0.01); 3.639, muon muon-Youdl You(lr =0.01); 3.629 and muon muon-Pola PolarExp(Ir =0.005); 3.588. Right: Validation loss vs. training iteration.
4.2 TRAINING GPT-2
We compare the performance of using Polar Express (Implementation 1) inside Muon against Jordan's 's (Jordan et al., 2024b) and You's (Cesista et al., 2025) methods. We train two architectures: 12) and GPT-2-Large (Nemba 1280, Rlaye
GPT-2-Small (Nembe
768, Maaye 12, Alaad 36, Thead 20), both with a vocabulary size of 50,257 and a context length of 1024. We train on IB tokens of the Web dataset AureleOOelleette et al., 2023) for epoch with batch size 32. All runs use mixed precision (bfloat16) on 4 H100 GPUs with the learning rate schedule proposed in Jordan al. (2024a)-a constant phase for the first 40% training steps followed by linear decay. All methods for the matrix sign computations are performed in bf bfloat] 1 6 precision and use five iterations. Following nano-gpt (Jordan al., 2024a), we assign Muon to all param- eters with at least two dimensions (e.g., excluding RMS norm parameters), except for embeddings, unembeddings, and positional encodings. These excluded parameters are optimized with A
Figures and 4 show the resulting in terms of validation loss for the GPT-La rge and GPT-Sma models, respectively.
both cases, muon-Pol Fxxxp achieves better validation loss than muon-Jordan muon-You. The advantage is remarkably consistent across all learning rates and epochs. While not shown Figures and 4, loss than the baselines, and the improvements in training loss are nearly identical the improve- ments in validation loss. Furthermore, since all three of these matrix sign methods are equally expensive (they all apply degree polynomial at each iteration), improved validation loss in terms of training steps also implies improved loss in terms of wall clock time. For figures displaying the improvements in training loss and wall-clock time, see Appendix H.2, Figure 11.
also achieves a better training
4.3 ABLATIONS
Accuracy of polar approximation We now explore how the accuracy of approximating polar(M) affects the optimization quality of Muon. Our main experiments with GPT-2 use it- erations. We trained GPT-2 Small with Muon using between 2 and 30 iterations of Polar Express instead. For comparison, we also implemented Muon with the exact polar factor, computed using torch. 1 Minalg. svd. Figure 5 shows the results. The left plot shows that when using only 2 or iterations of Polar Express, the final validation loss is worse than when using 5 or 6 iterations. However, increasing the accuracy of the polar approximation further-eves computing it exactly
*Code for our LLM training experiments is available at https::// http://ggg.t...tttteete//t sPttoppt/tree/pllel in the pola branch.
tettthu.b... com wwwmmodiiooooooooooooooooooo..oc /odi chi rag/

---

Published as a conference paper at ICLR 2026
Polar Method
Final Validation Loss
3.75
PoierExpress
Median Training Step Time
SVD
3.70
2
3.65
1
3.60
SVD
30 SVD
Iterations of Exprees
Therations Polar Express
Figure 5: Ablating the number of iterations of Polar Express used to implement Muon, and comparing to computing polar (M) exactly via SVD. Left: using > 6 iterations the SVD does not improve final validation loss. Right: Runtime of Muon is not sensitive to the number of iterations of Polar Express. s, but the SVD makes significantly slower. All runs use GPT-2-Small with Billion tokens of Wib data, learning rate 0.05, and weight decay 0.1.
3.5
3.06
TIADEL P Padniaap
3.4
3.3
Plas Valiaaiinn
3.00
Validatioo
3.2
2.98
3.1
2.96
3.0
291
2.9
10-1
10
Learning
Epoch
Figure 6: Training GPT-2-Large 10 billion tokens of FineWeb with weight decay 0.1. Best final validation losses were muon- muon-Jordan Jond (Ir = 0.002): 2.921, muen-You (Ir 0.002): 2.919 and muon-PolarEn (Ir = 0.002): 2.913.
with the SVD does improve the optimization quality. The right plot shows that changing the number of iterations does not meaningfully change the runtime of Muon; in setting, the runtime of computing polar(M MM) )) is dominated by the forward and backward passes. However, the SVD is so costly that using it doubles the runtime of each training step. These results validate the standard way of implementing Muon: using or 6 iterations an iterative approximation like Polar Express rather than computing polar (M) exactly. For further experiments supporting this conclusion. Appendix H.I, Figure 9.
Weight decay We also experimented with adding weight decay of 0.1 to the GPT-2 training runs, keeping all else the same. The results are presented in Appendix H.2, Figure 12. They are quite similar Figures and 4. We again find that PooaarExp outperforms the other methods.
Number of Training Tokens Our main experiments with GPT-2 use billion tokens of training data from FineWeb (Aroca-Ouellette et al., 2023). We now select a subset our training runs and extend them to 10 billion tokens. 10 billion tokens roughly matches the Chinchilla scaling rule for GPT-2-Large (774M params) and exceeds for GPT-2-Small, as per Table 3 in Hoffmann et al. (2022). Figure shows the results for GPT-2-Large with weight decay. (For GPT-2-Small, see Appendix H.2. Figure 13b). Polar Express still outperforms the baselines by a small but consistent margin.
10

---

Published as a conference paper at ICLR 2026
Acknowledgments This work partially supported by NSF awards 2045590 and 2234660.
Reproducibility statement A complete Pytorch implementation of our method given in Ap- pendix A. Details of experiments, including hyperparameters, are given Sections 4.1 and 4.2. Source code to reproduce our experiments is given in the supplementary materials and is available https: httpp://githuui /hittub.... com/ com/nodi chi Petg//PTTopp/// wetteeeepllloo in the polar branch. Proofs of all theoretical claims can be found in the appendices.
REFERENCES
N. I. Achieser. Theory of approximation. Dover Publications, Inc., New York, 1992. ISBN 0-486- 67129-1. Translated from the Russian and with a preface by Charles Hyman, Reprint of the 1956 English translation.
Samuel Aroca-Ouelletee Philippe Beaudoin, Guillaume Lajoie, Liam Paull, Joelle Pineau, Pascal Vincent, and Anirudh Goyal. Fineweb: Learning language models with high quality web data. In NeurlPS Datasets and Benchmarks Track, 2023. URL https: : //ta:///rxvvvvo ooooorogaaaoo. httg//b//23066 03061.
Michele Benzi and Ru Huang. Some matrix properties preserved by generalized matrix functions. Spec. Matrices, 7:27-37, 2019. ISSN 2300-7451. doi: 10. (0.1515/spma-2019900 URL https: aiddiiorgia 10. 151 wttpp/pm22222// 200---001100000
Jeremy Bernstein and Laker Newhouse. Modular duality in deep learning. arXiv preprint
arXiv:2410.21265, 2024a. URL httpttttp//// /ttpa/iii...gaa////////aaaaaa//a .12655.
Jeremy Bernstein and Laker Newhouse, Old optimizer, new norm: An anthology. arXiv preprint 249.. 20325.
arXiv:2409.20322 2024b. URL
A. Björck and C. Bowie. An iterative algorithm for computing the best estimate of an orthogonal matrix. SIAM J. Numer Anal. 8:358-364, 1971. ISSN 0029.1429. doi: 10/1137777 URL https: /oii. wtwgg/c..o.. 113 37 07/7080066
Christos Boutsikas, Petros Drineas, and Ilse C. F. Ipsen. Small singular values can increase in lower 1137
precision. SIAM J. Matrix Anal. Appl., 4518):1518-1540, 2024. ISSN 1995-4799,1095-7162 doi: 10.1137/23M1557200, URL https: http::////d..ogg11.1111////3333
David Carlson, Volkan Cevher, and Lawrence Carin. Stochastic Spectral Descent for Restricted Boltzmann Machines. In Guy Lebanon and V. N. Vishwanathan (eds.), Proceedings of the Eighteenth International Conference Artificial Intelligence and Statisties, volume 38 of Pro- ceedings of Machine Learning Research, pp. 111-119, San Diego, California, USA, 09-12 May 2015a. PMLR. URL httph:: http:///ccceeen....c.e.. mll www.pngw.....................e.e
press v3////3/////// car/sonllsn 5. html.
David E Carlson, Edo Collins, Ya-Ping Hsieh, Lawrence Carin, and Volkan Cevher. Precondi- tioned spectral descent for deep learning. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015b. URL /tttps::/pprcceeednggg.... neurip...........o.............e www.:eeewris....pppeeeeeeeeeeeee wwtp:ppp_ppppeeeeeeeeeeeeeeeeeee ffller 4ppp://055111f5a0aaac33333111333
Franz Louis Cesista, Jiacheng You, and Keller Jordan. Squeezing 1-2% efficiency gains out of muon
by optimizing the newton-schulz coefficients, 2025, URL htthtt///// httpt:::////lllunnnttt hwt.yykkuunnn.hhub....... io/
PL Chebyshev. Questions on smallest quantities connected with the approximate representation of functions (1859). Collected works, 2:151-235, 1947.
Jie Chen and Edmond Chow. A stable scaling of newton-schulz for improving the sign function
computation of hermitian matrix. Preprint ANI/MCS-P5059-01 14, 2014. URL https:// www. mcs anc..anl... www.aaan........oa wwtp/ppppee//S0SS55 1110099---1..11 pdf.
E. W. Cheney. Introduction to approximation theory. McGraw-Hin Book Co., New York- Toronto-
London, 1966.
11

---

abs/ Published as a conference paper at ICLR 2026
J. Douglas Carroll and Phipps Arabie. Chapter 3 multidimensional scaling. In Michael H. Birnbaum (ed.), Measurement, Judgment and Decision Making, Handbook of Perception and Cognition (Second Edition), pp. 179-250. Academic Press, San Diego, 1998. ISBN 978-0-12- 099975-0. doi: http://ttt.................ccaaa scienced rect
URL https:: www.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. J. Mach. Learn. Res., 12:2121-2159, 2011. ISSN 1532-4435,1533-799
Alexandre Eremenko and Peter Yuditskii. Uniform approximation of sgn æ by polynomials and entire functions. Anal. Math., 101:313-324, 2007. ISSN 0021-7670,,1665 8538. doi: 10.1007/
$11854-007-0011-- 1-3. URL https: // wted...rg/.....0/ssl.l...... /ttp:////oo...orr/110000///////0
Gene H. Golub and Charles F. Van Loan. Matrix computations. Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, Baltimore, MD, fourth edition, 2013. ISBN 978 100----414-0794444 1-4214-0794-9; 978-1-4214-085944
J. C. Gower and G. B. Dijksterhuis. Procrustes problems, volume 30 of Oxford Statistical Sci- ence Series. Oxford University Press, Oxford, 2004. ISBN 0-19-851058-6. doi: 10.1093/ 1cprof:oso/78801855 105819511111.... 10581.001..00... URL httLtttps//dee /tt.t/ttse////.....g//0/0..00/// 9780198510588 05118111111111111 001.0001. 000.00000000000000000000..
Ekaterina Grishina, Matvey Smirnov, and Maxim Rakhuba. Accelcrating newton-schulz iteration for orthogonaliiation via chebyshev-type polynomials, 2025. URL https: pettps://arxi...r
abs / 2506 10935.
Vineet Gupta, Tomer Koren, and Yoram Singer. Shampoo: Preconditioned stochastic tensor op- timization. Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 1842-1850. PMLR, 10-15 Jul 2018. guptal8a. html.
wwr.mmlr...... /
URL https: https///ccccceee.cccooo
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770-778, 2016.
Nicholas Higham. Computing the polar decomposition applications. SIAM J. Sci. Statist, Comput.,
1986. ISSN 0196-5204. doi: 10.1137/0907079. URL attps://
doi. org. 10. (tt.r////..011//911111 11117719797779
Nicholas Higham.
Functions of matrices.
SIAM, Philadelphia, PA, 2008 ISBN 978-0-
89871-646-7. doi: 9780898717778.
URL https://doi http://oii.rg/1..11111.1010 11.31/11.
Jordan Hoffmmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, Tom Hen- nigan, Eric Noland, Katie Millican, George van den Driessche, Bogdan Damoc, Aurelia Guy, Simon Osindero, Karen Simonyan, Erich Elsen, Jack W. Rae, Oriol Vinyals, and Laurent Sifre. Training compute-optimal large language models, 2022. URL 2203.15555
Elias Jarlebring and Gustaf Lorentzon. The polynomial associated with fixed number of matrix- matrix multiplications. arXiv preprint arXiv:2504.01500 2025. URL https: wwwxi.ivoo abs/2504. 01500.
Keller Jordan, Jeremy Bernstein, Brendan Rappazzo, @fernbear.bkky.socilll Boza Vlado, You Ji- acheng, Franz Cesista, Braden Koszarsky, and @Grad62304977. modded--anogpt: Speedrun- ning the nanogpt bascline, 2024a. URL https: modded-nane -andedd-nnogpt..
Keller Jordan, Yuchen Jin, Vlado Boza, Jiacheng You, Franz Cesista, Laker Newhouse. and Jeremy Bernstein. Muon: An optimizer for hidden layers in neural networks, 2024b. URL https:
//kellerjordan .iciiinnngittt github..t wtt.hb..o///aae tt/postss //uon//.
12

---

Published as a conference paper at ICLR 2026
Tetsuya Kaneko, Simone Fiori, and Toshihisa Tanaka. Empirical arithmetic averaging the pact Stiefel manifold. IEEE Trans. Signal Process. 61())8883 894, 2013. ISSN 1053-587X,1941 0476. doi: 10.1109/TSP:212.2226666 URL https: /tt::///ooiiii/ doi www.i.orgg1.... TSP 2012. 2226167.
Charles Kenney and Alan J. Laub. Rational iterative methods for the matrix sign function. SIAM Matrix Anal. Appl. 12(2):273-291 1991. ISSN 0895-4798. doi: 10.1137/0612020. URL
https:://bi..org/10..
Kimi Team, Yifan Bai, Yiping Bao, Guanduo Chen, Jiahao Chen, Ningxin Chen, Ruijue Chen, Yanru Chen, Yuankun Chen, Yutian Chen, Zhuofu Chen, et Kimi k2: Open agentic intelligence, 2025.
URL https: /twpiv/aaaii.a//aa//eaa org/
20534.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. URL 6980.
wbbp 14412.
Zdislav Kovářík. Some iterative methods for improving orthonormality. SIAM J. Numer. Anal., 7:
386-389, 1970. ISSN 0036-1429. doi: 10.1137/0707031. URL httpp:://ooi.oee 1137/0707000
Alex Krizhevsky. Learning multiple layers of features from tiny images. Technical Report TR-2009, University of Toronto, 2009. URL https: /tttp:::/www... /wwp///www/....///////////.o toronto.edu//krr....... learning- fereffeatuussse
Eunsang Lee, Joon-Woo Lee, Jong-Seon No, and Young-Sik Kim. Minimax approximation of sign function by composite polynomial for homomorphic comparison. IEEE Transactions on Dependable and Secure Computing, 19(6):371 100027771-777,, 1-7777, 2022. doi: 10. 109/TDDC..2211.1055111
R. B. Leipnik. Rapidly convergent recursive solution of quadratic operator equations. Numer. Math., 17:1-16, 1971. ISSN 0029- 599X,9999992955. doi: 10. 10.10///00035586/1 URL httpp:///doi. org/10. 000000///00196666..
Jingyuan Liu, Jianlin Su, Xingeheng Yao, Zhejun Jiang, Guokun Lai, Yulun Du, Yidao Qin, Weixin Xu, Enzhe Lu, Junjie Yan, et al. Muon is scalable for LLM training. arXiv preprint 17X1:22502.16982. 2025. URL https:ttps/..mm //t//aaxvvvvrrg.e/ee wttpgagb//./////// /ttt:s///axxd..orr/aa/////// 16982.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization, International Confer- ence on Learning Representations, 2019. URL https:ttp//... /oppereveeeeeeeeeeeeeeeeeee
htt.ht:::///peereeiiie net forum? d=
Bkg6RiCqY7.
Modula. nowton-schule algorithm - jiacheng's six-step method. https: //docs modula. systems / a gori ceoorthhmme tmms hww.://wwwwshhooeeeee / aiaccheeng-eee --ssixsseep. 2024. Ac-
cessed: 2025-05-19.
Yuji Nakatsukasa and Roland W. Freund. Computing fundamental matrix decompositions accurately via the matrix sign function in two iterations: the power of Zolotarev's 's functions. SIAM Rev,
58(3):461-493 2016. ISSN 0036-1445, 10955,000. doi: 10.1 137/ 4099733344 URL https:
//ddi. ott.e00.rg11..0 113711
Yuji Nakatsukasa and Nicholas J. Higham. Backward stability of iterations for computing the polar decomposition. SIAM J. Matrix Anal. Appl., 33(2):460 479, 2012. ISSN 1895-4798 doi: 10. 10857544. URL htt.t:::. : /ooioooo//////// wor./orgg1/.. 1137/1108177544
Yuji Nakatsukasa and Nicholas J. Higham. Stable and efficient spectral divide and conquer al- gorithms for the symmetric eigenvalue decomposition and the SVD. SIAM J. Sci. Comput.,
A3355-A1325-11349 2013. ISSN 1064-8275.1095-7197 doi: 10.1137/120876605. URL httpp:://df .org/orgggo... 1137/120876605
Yuji Nakatsukasa, Zhaojun Bai, and François Gygi. Optimizing Halley's iteration for computing the matrix polar decomposition. SIAM J. Matrix Anal. Appl., 31(5):2700-2700, 2010. ISSN
1055-449811095-711 doi: 10.1137/090774999. 090774999.
URL
https: : /ttp////.i////oo
13

---

Published as a conference paper at ICLR 2026
Herbert Neuberger. Exactly massless quarks the lattice. Phys. Lett. B, 417(1-2)111--44, 1998. ISSN 1377-2269, doi: 1010016600370226559970113336 URL https: //doi.
(97) (117))1133553 01368-3.
Ricardo Pachón and Lloyd N. Trefethen. Baryeentric-eeeer.ggrriiiii for best polynomial approx- imation the chebfun system. В/T; 4944):2-11111 2009. ISSN 1572-9983 doi: 10. URL
1007///05433000-00220000 / ott.a//aaa////ii..or///..../0//0 0000000002--00000000000000000000 0240
T Parks and James McClellan. Chebyshev approximation for nonrecursive digital filters with linear phase. IEEE Transactions on circuit theory, 19(2):18-1994 1972. doi: 10.1109/CCT.1922 1083419.
Michael S. Paterson and Larry J. Stockmeyer. On the number of nonscalar multiplications necessary to evaluate polynomials. SIAM J. Comput., 2:60-66, 1973. ISSN 0097-5397, doi: 10.1137/ 0202007. URL httpt::t///........ //t.////oii/////o httg://ooooooro///o///oo 1117000000000007..
Thomas Pethick, Wanyun Xie, Kimon Antonakopoulos, Zhenyu Zhu, Antonio Silveti-Falls, and Volkan Cevher. Training deep learning models with norm-constrained Imos, 2025. URL https:
/
07529.
Yehonathan Refael, Guy Smorodinsky, Tom Tirer, and Ofir Lindenbaum. Sumo: Subspace-
aware momeent-rtthooonlizziiii for accelerating memory-efficient Ilm training. arXiv preprint arXiv:2505.24444, 2025.
Artem Riabinin, Egor Shulgin, Kaja Gruntkowska, and Peter Richtárik. Gluon: Making muon & scion great again! (bridging theory and practice of Imo-based optimizers for Ilms), 2025, URL https: //tpiii////xvv../////// org/ /ttp////....... 13416.
Ishaan Shah, Anthony M Polloreno, Karl Stratos, Philip Monk, Adarsh Chaluvaraju, Andrew Hojel, Andrew Ma, Anil Thomas, Ashish Tanwer, Darsh J Shah, al. Practical efficiency of muon for pretraining. arXiv preprint arXiv:2505.02222 2025. URL /tttt:s///axiv.oot abs/ 2505. 02222.
Hao-Jun Michael Shi, Tsung-Hsien Lee, Shintaro Iwasaki, Jose Gallego-Posada, Zhijing Li, Kaushik Rangadurai, Dheevatsa Mudigere, and Michael Rabbat. A distributed data-paraliel Py- Torch implementation of the distributed Shampoo optimizer for training neural networks at-scale. arXiv preprint 230.002909 2023. htts::/ arxiv. org abs/abs/////aa.... htt.b///3/// 2309. 06497.
Egor Shulgin, Sultan AlRashed, Francesco Orabona, and Peter Richtárik. Beyond the ideal: Ana- lyzing the inexact muon update. arXiv preprint arXiv:2510.1993: 2025.
Attila Szabo and Neil Ostlund. Modern quantum chemistry: introduction to advanced electronic structure theory. Courier Corporation, 1996.
Lloyd N. Trefethen. Approximation theory and approximation practice. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, extended edition, 2020. ISBN 978- 178-1-611975- 93-2.
Nikhil Vyas, Depen Morwani, Rosie Zhao, Itai Shapira, David Brandfonbrenerr, Lucas Janson, and Sham M. Kakade. SOAP: Improving and stabilizing shampoo using adam for language modeling. In The Thirteenth International Conference Learning Representations, 2025. URL https: // //opereeview.o formm? d=TxxhxxppN.
Greg Yang. James B. Simon, and Jeremy Bernstein. A spectral condition for feature learning, 2024.
URL pestts::/arxiie
abs/23100 17813.
Zhenyue Zhang, Hongyuan Zha, and Wenlong Ying. Fast parallelizable methods for computing invariant subspaces of Hermitian matrices. J. Comput. Math., 25(5)5583-594, 2007. ISSN 0254- 9409,1991-71339 URL hwtphhtpt/wwee www..jstor.....oo......o........ http://ggg///oeee.ooooo httss/aabee44333330 43693395.
14

---

Published as a conference paper at ICLR 2026
CONTENTS
1
Introduction
1.1
The Muon Method
1.2
Computing the Polar Factor
1.3
Contributions
2
Approximations by Compositions of Polynomials
3
The Polar Express
3.1
Greedy is optimal
3.2
Finding the optimal polynomial for each iteration
3.3
Upper and lower bounds on the singular values
3.4
Finite precision considerations
3.5
The algorithm
Numerical Experiments
4.1
Convergence of Polar Express
8
4.2
Training GPT-2
4.3 Ablations
Code for Polar Express
16
Related Work
17
Proof of Theorem 3.1
19
D Proof of Theorem 3.3
22
E Proof of equivalence between (5) and (6)
23
F Remez algorithm
23
Finite precision considerations
26
Additional Experimental Results
26
H.l Convergence of Polar Express and Its Impact Muon
26
H.2 Training GPT-2
30
H.3 Image Classification
30
I Initialization for Matrices with Large Spectral Gaps
30
J Fast Polynomial Iteration for Rectangular Matrices
35
J.1 Application to Muon
37
15

---

Published as a conference paper at ICLR 2026
A CODE FOR POLAR EXPRESS
Implementation gives Python implementation of the online stage of Algorithm for degree 5, which in our numerical experiments. It uses hard-coded polynomials generated from Implementation and incorporates a numerical safety factor of 1.01 as described Section 3.4. This implementation is designed for ease of use. is short, it has no dependencies besides PyTorch, and it a drop-in replacement for previous implementations of matrix sign methods (Cesista et al., 2025: Jordan et al., 2024b), such Modula (2024).?
Implementation 1 Python code for Polar Express of degree = 5.
From itertcols import repeat
import torch
coeffs 1ist =
_ostfs_listtt
(8. 0t.287212018145555
0...9998865110988833
17.3003873125300933
(4. (4.107059111542203,,
2.9978499166379906
0.54484310829266011
(3. 93.9486908333822946, -2.908902115992999,
0.55181913993771337)
(3.31841965737060155
-2.488488024314874,
0.51004994012372),
(2.300652019954881,,
-1.6699039445747499,
0.41897331195256733
1.8813011407787398
-12277955221195868,, 0.37680408946524669),
11.8750018888534419
-12500011445399987,
0.375000164547444811
(1.875, -1.25, 0.375),
# subscquent coeffs equal this numerically
safety factor for numerical stability (but exclude last polynomial) coeffs_list. [{a 1.01, 1.01x*3, 1.01vx5) for (a, b, c) coeffs_list::-1 + [coeffs_list[]] [-11]
@torch.
def FolarExpress (G: Tensor, steps: int) torch. Tensor: assert n.ndim
X - 6.bflcat16() I for speed
if G.size
(-1): x X.mT this reduces FLOPS
X = X / (X.norm (dim= -1), keepdim=True) 1.01 +1e-7) hs beefssliit[:steee
list(
(coeff_ssil
steps len (coeffss__iss
for a, b, c in hs:
A = X X.mI
A A
+
<- bx*3 cx
if (- G.size X X.mT return x
Implementation 2 gives a Python implementation of the offline stage of Algorithm 1. This code was used construct the coefficients of the polynomials given Implementation 1, which in turn were used in our Muon experiments (Section 4.2). It uses f = 10-3 and by default. It incorporates Algorithm and the finite precision modifications described in Section 3.4.
Implementation 2 Polar Express, Offline Stage
from math import inf, sqrt import numpy as np
def @odimal_qqintic(]] u): assert <- 1 U if - 5e-6 <- 1 u: Above this threshold, the equioscillating polynomials
is numerically equal to...
return (15/8)/u, (-10/8) (u**3), (3/8) (u**5) This initialization becomes exact as 1 -> u
= (3+1 + u)
"Code including Implementations and can also be found at pettps://ge thub. com /NoahAmshhh/////
Pol arkxprsss
16

---

Published as a conference paper at ICLR 2026
3+u) E, old_E inf, None
while not old_ or abs (old_E - E) 1c-15:
old_ =
LAS np.array
[1, 1**3, 1x*5, 11, [q, qu*3, q -1], [r, r**3, rew5, 11, [u, u+*3, u -1l,
1)
a, b, G, np. lenallnng...... ww..inalg.sccveellle sclve(HHS) np.ones (4))
4 r up. sqrt (-33b np.array([-l, 11)
sqrt (9 - 20naxc)) (10xc))
return fa)at(a), float (b), (ccat(c)
def ootimal_comosstioonll num_iters, cushion=0. 02w00000000000077772222202777
=
coefficients []
for range (num_iters):
a, b, c apt.iml__qiintin (1, aushion<u), u)
# Due to cushioning, this may be centered around 1 with # respect to 0.024*u, u. Recenter it around 1 with respect that cLp cxp(u) 1:
u, meaning find
pl a+1 + b+l+*3 c+1++5
pu axu + beu++3 c+U++5
rescalar - (pl + +
-= rescalar; b - rescalar;
rescalar
Optionally incorporate safety factor here:
# a 1/=1, /= 1.01**3; c /= 1.01*x5
ierFiicients.ppend..
b, c))
1 = a+1 + b*1*+3 + c*1*+5
2 return coefficients
print (oopt i ma1 wwwposiiiooiiooooccccoooo.cooo.e (le-3, 10), sep-"in")
B RELATED WORK
Computing polar(M (M) an important and longstanding problem in numerical linear algebra, with applications spanning electronic structure calculations, lattice quantum chromodynamics, orthogo- nal Procrustes analysis, parallel algorithms for computing the SVD, and beyond; see e.g. (Higham, 1986: Kaneko et al., 2013; Douglas Carroll & Arabie, 1998; Gower & Dijksterhuis, 2004; Neu- berger, 1998; Szabo Ostlund, 1996).
Newton-Schulz and polynomial Padé methods. The earliest methods in the literature poly- nomial iterations like (2). Several nearly simultaneous papers introduced the family of polynomial Padé iterations, comprising Sowton-Schull and its higher-degree analogues (Kovářík, 1970: Björck & Bowie, 1971; Higham, 1986; Leipnik, 1971). These higher-degree methods are also sometimes called "Newton-Schulz"" when doing so, we will specify the degree for clarity. In these methods, each iteration refines the current approximation X, by applying a low-degree odd matrix polyno- mial, where any odd monomial x >> 220+1 defined for rectangular matrices by the formula X, X, (X," X,)* Our Polar Express method also takes this form, though unlike Newton- Schulz, it changes the polynomial at each iteration.
The polynomials used Padé methods are chosen to match the value and first few derivatives of (zg)(x) at the points x = +1. For instance, the update rule of the third method in this family is defined by p(z)
To (35x 35 + 21x 5x"), which is the unique degree-7 polynomial satisfying p(=1) = +1 and p' (=11) = p=(=1) (=11= = These methods converge so long as
all singular values of X, lie in (0, 1). condition guaranteed by the initialization of (2). Furthermore,
17

---

Published as a conference paper at ICLR 2026
the order of convergence of the degree 29 + method is q 4+1 (Björck & Bowie, 1971). In particular, the Newton-Schule method (q= = 1) converges quadratically.
Newton's method and rational Padé. In the numerical analysis literature, polynomial methods were succeeded by rational iterations like Newton's 's method (Higham, 1986), defined as follows":
1 (Xx xX 2
X M
X
(14)
Newton's 's method also converges quadratically. Like Newton-Schulz, it works because the rational function r(x) = 2/2 (z x=') has stable fixed point 1; unlike for Newton-Schulz this point is global attractor for the whole positive real line. At first glance, Newton's 's method has nothing to do with the Padé iterations discussed above. However, after change of variables Y, = X,7, it can be reinterpreted as Y++1 = 2 + Y, Y,) which is sometimes called inverse Newton Observing that r(z) I+z² satisfies r( =1 and 1=1) O, we see that (inverse) Newton also a Padé method, though a rational rather than polynomial one. In fact, given odd degree 2q,, + for the numerator and an even degree 2qa for the denominator, there is a unique rational function that matches the value and first On + Qu derivatives of sign(x)) at x #1. This directly yields Padé method for computing polar(Mi whose order of convergence is 4n + qu + 1. For instance,
2m
r(x) 3z+=" I+327 is called Halley's 's method, which converges cubically. When qa = O, we recover the polynomial Padé methods.
There are two main weakness of Newton's method and the Padé iterations: slow convergence in the initial phase and the need to compute explicit inverses. To accelerate initial convergence, Higham popularized the technique of rescaling the matrix after every Newton iteration (Higham, 1986). Intuitively, rescaling X, so that max = 1/ centers the spectrum around 1, where convergence is fastest. Several casiuy-computable choices of scaling factor exist to accomplish this approximately. Note that this rescaling scheme would fail for Newton-Schule which likewise suffers from slow initial convergence but which would diverge if Fmax > 1.
Computing matrix inverses is difficult to parallelize and implement stably in low precision arith- metic However, a trick was developed for stably computing many rational methods without explicit inverses; QR decompositions can be used instead (Nakatsukasa al., 2010; Zhang et al., 2007). Applying this trick to Halley's 's method and combining with special rescaling scheme yields the QDWH -as-d dynamically weighted Halley) method, which converges in just six iterations for any reasonably conditioned matrix (Nakatsukasa et al., 2010).
Adaptive rational methods from optimal approximations. A landmark 2016 paper introduced a new paradigm to design iterative methods for computing polar(M) (Nakatsukasa & Freund, 2016). The main insight is as follows. Padé methods choose the update rule to be an approximation to sign(x) of a given degree that is optimally accurate in the neighborhood of x = 1. Instead, we should choose the approximation to sign(x) that is optimal over an interval [E. 1] C R2o that contains the singular values. Moreover, after each step of the algorithm, the range of the singular values changes; therefore, we adapt the update rule at each iteration to match the new interval. When the range of the singular values is large, this approach ensures that the update rule shrinks it quickly as possible. As the algorithm proceeds and the interval shrinks to small neighborhood of 1, the update rule approaches that of Padé method, maintaining the same high order of convergence as has.
Within the class of odd rational functions whose numerators and denominators have degree 24 + and 2q. respectively, an explicit formula for this optimal approximation to sign(z) on any interval [e, 1] was found by Zolotarev. It was shown that these rationals have remarkable convergence prop- erties for any q (Nakatsukasa & Freund, 2016). For q = 1. this optimal approximation coincides exactly with the dynamically weighted Halley's 's method (QDWH) referenced above. For even faster convergence than QDWH, (Nakatsukasa & Freund, 2016) proposed the Zolo-pd method, which uses q = 17. Finally, these methods all admit the same -R-ba implementation trick as QDWH.
daptive polynomial methods. In this paper, we adopt the paradigm of Zolo-pd (Nakatsukasa & Freund, 2016) but with polynomials rather than rationals of degree (2g + 1,)4). This choice
*Our description of Newton's 's method and other rational methods assumes square non-singular M. Non- square problems can be reduced to the square case by initial QR decomposition, but this is not an option for purely polynomial methods like ours.
18

---

Published as a conference paper at ICLR 2026
avoids the need for QR factorizations, relying solely GPU-friendly matrix-matrix multiplications in low-precision arithmetic. While this class of methods has not been fully developed in the numer- ical analysis literature, similar ideas have been rediscovered in different guises. In an unpublished manuscript that predates Zolo-pd, Chen & Chow (2014) descrihe a rescaling strategy for Newton- Schulz. Though motivated differently, their method is equivalent to ours for degree-3 polynomials (unlike our work, they do not consider general odd degree). They also observe numerical instability that prevents the method from converging all the way machine precision. Using the insights of Nakatsukasa & Higham (2012), they propose simple mitigation for this issue that we adopt in Section 3.4. Our work gives the approach from Nakatsukasa & Higham (2012) a stronger theoretical foundation that connects to the paradigm of Zolo-pd. Concretely, we prove that choosing an optimal polynomial at each iteration leads to a composed polynomial that is globally optimal in the sense of (5).
Independently, a group of cryptographers developed a similar method for approximating the scalar function sign(z) in the context of homomorphic encryption schemes (Lee al., 2022). Their focus mainly on tuning the analogues in their setting of the polynomial degree and number iterations, whereas focus on demonstrating optimality and efficiently constructing the update polynomials for degree 3 and 5. addition, we consider matrix-valie inputs in low- low-precision arithmetic- not scalars in exact arithmetic-and we demonstrate our method's effectiveness within the Muon algorithm for training deep neural networks.
Application within Muon. The designers of Muon realized that, due to the extreme efficiency requirements and lax accuracy requirements of their setting, rational- based methods from the nu- merical analysis literature are inapplicable. However, polynomial-basse full advantage of GPUs because they use only matrix-matrix products in half-precision arithmetic, not inverses or QR decompositions. The preference for speed over accuracy motivates methods that aim to quickly produce coarse approximations, even the cost of asymptotic convergence. Examples include the proposals of Jordan (Jordan et al., 2024b) and You (Cesista et al., 2025), as discussed in Section 1.2. Like Chen & Chow (2014), Jordan found that convergence in the initial phase can be accelerated by choosing update rules that have a large derivative near zero, so as to increase the small singular values much as possible each iteration. You furthermore chose to use different update rules at each iteration, allowing extra flexibility to tune the trade-off between speed and accuracy. Both used degree-5 polynomials that were found through gradient descent on heuristic objective functions. These proposals were previously compared to Newton-Schulde but never to Nakatsukasa & Higham (2012). We find that our method (which generalizes Nakatsukasa & Higham (2012) 2)) outperforms them all.
iteration schemes can take
Finally, we remark that concurrent work of Grishina, Smirnov, and Rakhuba also proposes an adap- tive polynomial method that generalizes Nakatsukasa & Higham (2012) and applies it to accelerat- ing Muon (Grishina et al., 2025). Like Nakatsukasa & Higham (2012), this work does not establish global optimality of the composed polynomial as we do in Section 3 or address finite precision considerations.
C PROOF OF THEOREM 3.1
The aim of this section to prove Theorem 3.1. We begin with a result that provides a few essential properties for the the polynomial solving (6) when T 1. This result is known Chebyshev's theorem (Chebyshev, 1947) or the equioscillation theorem (Trefethen, 2020, Chapter 10).
Lemma C.1. Let d = 2q + 1 and u, é > O. Consider the problem
min max | p(x))-
(15)
peredd (E(M,))
"Joddn al. (2024b) actually compares 2z 3 400, whereas the true degree-5 Newton-Schuls polynomial is (15z 10z + 3/8) However, the difference in performance is negligible for the first few iterations.
19

---

Published as a conference paper at ICLR 2026
There exists a unique polynomial p" Pagd solving (15). Furthermore, p" is the unique solution to the above problem if and only if there exist q + distinct points (zo..... Tg+1} C [e, u| such that
(xt) = vt-1)
max [f
[ - p (x)), for = .,....q
for 17 = 1 or 7 = -1.
Proof. A discussion be found in Eremenko & Yuditskii (2007). Here we include formal proof for completeness.
By Chebyshev's Theorem (Achieser, 1992; Chebyshev, 1947; Cheney, 1966) it is sufficient show that podd satisfies the Haar condition: any non-zero p E podd = span [ (z...., ........ x. 20+1 can have most q roots in [e, u].
Since deg(p) = = 2q + 1 we know that can have at most 2g + 1 roots in R. However, since p(0) and p(x)
know that p has one root zero, and the remaining roots come in symmetric pairs (z, -x). Because of this, p can have at most roots in the positive orthant, and thus can have most roots in [4, u| (0, 0o). Hence, good satisfies the Haar condition, which yields the desired result.
p( (-==))
The proof of Theorem 3.1 will be by induction on T. We begin by establishing the base case, T 1, which is handled by the following result.
Lemma C.2. Let u, E >>> 0 and define
p" == arg min max p(x))- PEPj sE(t,u)
Then
p" (1) min p" (z),
max (z()) pf), and max
|1 - p" (x) = (().
=E(f.u]
=E(t.u
Proof. Throughout the proof assume d 2q 1. We begin with proving
(e)
min p" (x). 2E(t.u.
Consider the polynomial e(x) == 1 - pr (x)) The proof will contain three steps. We first rule out the trivial case that p" * 0, since p(x) ट,ッ would then be better approximation. Hence, p² cannot be the zero polynomial.
Step 1: e(x) has exactly q stationary points inside the open interval (l, u).
Note that e(x) has at most 2q stationary points in R, since its derivative e'(x) is a polynomial of degree 2q. Furthermore, since p* odd, have that e'(x) -p'(zi (x) is even of degree 2q. and thus can have at most stationary points contained in (0, +oo). Hence, there can be at most 9 stationary points of e(x) inside the interval [e, uj.
Iq+1 [f, u] where e(x) is maximized or minimized in [E, u]. These points are either stationary points or they are endpoints of the interval [e, u4). Let next be the number of stationary points and nsa be the number of endpoints in the set {ro, Since a point can be both a stationary point and an endpoint we have & + 2 S Rend + Tlstat However, Mend < 2 and Nlsta < q. which follows from the previous paragraph where we showed that there are at most y stationary points e(x) in [e, u). So nend + Nsstat S q + 2, and consequently we must have nand = 2 and nssaa q, as required.
By Lemma C.I there are q + 2 points 2o. Iq+1
xg+1)
Step 2: x is maximum of e(x) on the interval [e, u]
By Lemma C.1 and the discussion from Step 1, we know that le(x) is maximized at q + points inside [e, u] and q of these points contained inside the open interval (e, u). Hence, must
20

---

by Published as a conference paper at ICLR 2026
either be maximum minimum of e(x). We will show that l must be maximum by contradiction.
Suppose e was minimum of exx on [e, u]. First note that is trivially non-negative on [e, ud, or else p(z) = would be a better polynomial. Hence, since p" p0(0) we must have p²' > for some [0, f. else the zero polynomial p(x) would be better approximation. Hence, for some E [0, E] we have (6) < O.
We must also have el( 2 0 else é is not minimum of e(z). Since e'(6) < 0 for some E [ and e') 2 O, by the intermediate value theorem there exists point x² E [0, such that e' e'(x*) 0. However, by the discussion above know that all stationary points of are contained inside the open interval (e, u). Hence, I = l cannot be a minimum of e(x) on [e. uj. However, by
Step know that the endpoints of [e, u] must be either minima or maxima e(r). Hence, x is maximum of e(x) on [e, u].
Step 3:
Obtaining the desired equalities
Since e(r) has maximum in [e, u| at x L, have p²( min p" (x)). The other two equalities (E(f,u) are immediate consequences of the equioscillation property of p" Lemma C.I and that L is minimum of p" over the set [e, u].
With the above-mentioned result hand, we are ready to prove Theorem 3.1.
Theorem 3.1. Let be odd and define & and 4n = u. For = I,..., T define
arg min max
|1 - p(z)). l++1 = min Pc(x), ue+1 = max Pc(x) sEll.,e
(8)
pEPadd sE(fe,
=E(e..u)
The resulting composition p= :- pT о PT o o pi is optimal and the error is given by:
max
p=()
min pmpropy op EE(, erepa
max
|1 p(z) 1 tr+1-
(9)
(E((,u)
Furthermore the new error, lower and upper bounds can be computed through
- u+ = 2 le+1,
and =E(P..
max
|1 - pr(x)) = l111.
(10)
Proof. The proof of (10) an immediate consequence of Lemma C.2, since for each 1,..., is the optimal approximation in Podd to +
We now proceed with the proof of (9), which will be by induction. The proof for 1 is an immediate consequence of Lemma C.2 and we also have (e) = lz by (10). Now suppose the result true for all < - Thus
g(x) PT-1
opp (xx)
is the optimal solution of (9) for T 1. For t = 1, ...., T 1, note that the image of pt on [&, u.] exactly [Le+1, M441] by Lemma C.2. Hence, the image of on [e, u] is [Er, wr]. Furthermore, by Lemma C.2 we also have s()) Ey. Pick any such that + g and
f pT- P1.
for some ,I.....Prr Pr-1 podd Let the image of on [e, u] be [a, b]. We will prove that g S ur contradiction.
2 Suppose & > fr. Define e o-T and cf satisfies
Then, the image of the scaled function ef on [e, u) is [ea, ob)
max 11 - cf(x)) | = max {1 ca, ch - 1} b =E(E.u]
21

---

Published as a conference paper at ICLR 2026
Recall by inductive hypothesis, have max - g(x)l se((,u)
ur where the second
equality holds by (10). follows that
a
lt
>
b
b
20
+ b 1
17
+ b
max [E
1 g(x))
max |1 - cz(x)) zE[C,u]
which leads to a contradiction to our inductive hypothesis that is optimal. Hence, we must have fn
Consequently, using that ; S wy we will show for any Pr e Padd and for any f = Pr-1 ..... op,, that pr o cannot be better approximation than pr o g. In particular, we have
max 1 ÿT (f())) 2 min max 1 p((())) #E =Ef.uu)
PEP (E((.)
min max
- p(x)
PEP me(a
min max |1 p(x)
PEP me(a/b,1]
2 min max pEl =E/t//ur
|| - p(z)
min max PEl (f(.we))
|1 p(x)
min max | - r(g(())) PEP [f(M.)]
max 11- xE(r
= -1- = Lr+1s
where the second and third equality follow by changing variables y z/b that
min max 11 p(=) min max |1 p(by)) min max 1 p(v) pEP yEla/b.11
PET #E|a,b]
PEP vEla/b,11
and this last equality follows because the space IP is invariant under input rescaling: that is, for any b + 0, the map - be preserves the space span x² .............. ....., This concludes the proof.
D PROOF OF THEOREM 3.3
In this section we provide the proof of the convergence guarantee stated in Theorem 3.3.
Theorem 3.3. Let M be matrix normalized so that a (M) [4. 1]. Let X, ppp (M), where p" is the polynomial from Theorem 3.1 with 2q + 1. Then, we have
polar (M) Xr S |1 2- (0++)
(12)
Hence, for d = 3 and d the method converges quadratically and cubically, respectively.
Proof. Define
p"
arg min
|1 p(x))
pmpropT opl (t
P.EP-
22

---

Published as a conference paper at ICLR 2026
Then Algorithm returns = (M). Let h Pg be the /q/0] Padé-approximant to (1-z)-1/2 (Kenney & Laub, 1991, Section 3) and define p(x) zh(1 =²) Padd Define f = pp..... p as
the composition p with itself I times. Then, by Theorem 3.1, (Kenney & Laub, 1991, Theorem 3.1), and f()) 2 0 for æ 2 0 we have
|| sign(M) Xr|12
max aE(4,]]
p" (x)
max
f(z)
1 - =1=11(
max aE(4,]]
1 f())
|| - (d((1)))²
as required.
E PROOF OF EQUIVALENCE BETWEEN (5) AND (6)
In this section we provide a proof for the equivalence between (5) and (6). is sufficient to show that for any fixed polynomial p we have
[|polar(M) - p(M))l2 max |1 p(x)| = (M)
MER
=E(f.a)
For any fixed M, by the unitary invariance of the spectral norm we immediately have
[poola((M) - p(M)l22 = max 11 p(o,)) max o.Ea((M)
p(x))
zE(t,u)
Consequently, E1 S2..
Suppose that x² E [e, u] is chosen so that - [0=p((=)] )| max 1- 1-p(()) Without loss of
generality, assume 2 Letting M = x UV for any matrix U E Rmxx and V E RRax with
orthonormal columns, and noting polar (M) = UVT yields
E 2 polar p(M)
I, p(x=') Inll
1 p(x²))
max |1 - p(x)) 62 =E(
Consequently, E1 2 Hence, desired.
F REMEZ ALGORITHM
In this section, we show in detail how to solve (13). By Theorem 3.1, these solutions give the update rule for a single step of Polar Express. We give a closed form solution for = 3. We then describe how the Remez algorithm (Pachón & Trefethen, 2009; Parks & McClellan, 1972) can be used to approximate for arbitrary d. We then present Algorithm 2, simplified version of Remez for solving (13) with d = 5. Recall (13):
arg min max [1 - p(x) pEPaad zE(t.a]
We begin with the case when d' 3. We seek a polynomial of the form p(x) = az + bz The Equioscillation Theorem (Lemma C.I) stipulates that p must have equioscillating set of size 3. For p to achieve its maximum error at a point x, x must be local extremum of p(z) interval [e, u]. Thus, for to be eligible for membership in the equioscillating set, it must either be a true local extremun of p(x) 1 that happens to lie in [E, u), or else one the endpoints E, u. However, because is odd cubic, it has at most true local extremum on R Thus, to build
on the
23

---

Published as a conference paper at ICLR 2026
an equioseillating set of three points, must include p's unique positive local extremum and both endpoints. This local extremum of occurs Therefore, we seek a, b such that
=a
- E,
P
=1+E,
p(u) =1
(16)
3b
for some E. This a system of three equations three variables. The solution p(x) = ax + bx is most easily expressed as follows. Let PNs (a)
2 2 Then
3 and B + Eu(l + u)a3
p(x) (pns(ar), where
V 2² + tu +
(17)
One can verify that this polynomial satisfies the equioscillation condition of (16), with
36
and E = 3 - Therefore, it must necessarily be the optimal approximation from Pgdd Note that for xr (os()a)) the same polynomial derived in Chen & Chow (2014).
Unfortunately, for larger d, finding closed form expressions for optimal approximations from P becomes challenging, and we know of no closed form solution. However, we can approximate the optimal polynomial using the Remez algorithm. Let d = 2q + 1. Again recalling Lemma C.I, the optimal polynomial must satisfy the equioscillation property The Remez algorithm finds the equioscillation points A = {øo.... Zq+1 from Lemma C.I by iteratively refining a sequence of trial points ( {-o
of q + 2 points, in (16).
x 9+11 so that ( converges to A. From the sequence of trial points (A) the algorithm also finds sequence of polynomials p(*) so that pf*) converges to the optimal polynomial. The convergence is very fast, and usually 10 iterations sufficient converge to the optimal polynomial up double precision machine epsilon (Pachón & Trefethen, 2009). More commonly, the Remez algorithm is used to find optimal polynomial approximations to general continuous functions where d 100 or even d 1000. However, because the polynomial we build to approximate sign(x)) is a composition of polynomials, each of which has low degree, in our setting the degree is small, usually d = For algorithm simplifies significantly. We now describe this simplified algorithm.
the Remez
We first choose an initial set of trial points A(), which ideally should come close to satisfying the equioscillation property. From Lemma C.1, the unique optimal approximation p* E pgdd satisfies the equioseillation property four points in [ uj. Since the function we wish to approximate is constant, the equioscillation points must be extrema of p* on uf, Because px is a odd quintic, it have at most two local extrema the positive real line, and thus at most two local extrema on [e, uj. The other two cquioscillation points must therefore be the endpoints l and u. Since we know that l and must be equioscillation points 3
(k) f and 23 (k)
always set zig
for all k.
We initialize (1) and x22 to Ve + 3u and 1 l + 4 u, since we observe that as l + u these are approximately the other two equioseillation points.
We now show how to refine candidate set of trial points (A) to produce ((++1) well as an (k) (k)
approximately equioscillating polynomial Pk- For any fixed of trial points {L, .22 u). we can find d odd polynomial pk (x) a&d + by.² + -x.25 that satisfies
Ps (c) = Ee, pa (x(
(k)) = + Ek. px(x (K)) = Ex, (=( = + Fx (18)
for some Eg by solving a linear system in ake by, Ck and Eg. This can rewritten as follows:
1
0a
(k)
(x
ss
-1
by
x2
(k)
(x2
Chipa
1
Ck
=
(19)
Ex
(M) were the extrema of the enror function ex(x) = ps( on [e, u]. then they would be an equioscillating set for Ps. and Pa would be the solution. Therefore, to refine A we find the extrema of ex (x) = 1 - Pk (3)). These occur at e, u and the roots of (x). Setting (x(x))
24

---

Published as a conference paper at ICLR 2026
yields the quartic equation 5cg=** + 3b + ag O, whose two solutions are given explicitly by
the quadratic formula after the substitution = x² We set s1 ((++1) and x22 ((++1) to be the solutions
to this equation and let ((++1) = (L, a. (k+1) (k+1) .u. We repeat the procedure until [&| == I2
max [1 pk (a)) x max [1 Pk+1 (x)] =: Ek+1/- (E(f.u)
=E(f.u)
We note that the matrix appearing in (19) is Vandermonde matrix. Vandermonde matrices become notoriously ill- coldionitto as the degree grows large (Golub & Van Loan, 2013, Section 4.6). How- be small, there is ill-conditioning due large degrees.
ever, since in setting we choose Instead, we observe ill-conditioning when & = u. However, as 8/u + the optimal polynomial will 8 (15 10(x/u)² 3(z/u)") which can be verified by noting that as f/u - 1 all equioscillation points Fo. 21, 22, 2g must converge to u. For general d = 2q + 1, the polynomial will converge to (a/E))(11 - (z/ where h E Pg the [g/0) Padé approximant to (1 - z)!/² (Kenney & Laub, 1991). In fact, this polynomial is extremely close to the optimal polynomial for sufficiently large f. To see this, let p" be the optimal approximation from P2 and let p(x) = z/N (15 10(x/u)² + 3(x//)))) Then,
converge to the polynomial #/"
8
max (x*(2) p()) S max =E((,u]
p(x)] + max p² (x) xE((,u)
zE((,u]
< max |1 p(x)
#E [f.
< 2 (1 - t/u1²
where we invoked (Kenney & Laub, 1991, Theorem 3.1) and the fact that p* is the optimal approx-
imation to x + 1 from Pgdd. Hence, when f/u 2 1 d where Edouble R 1.1 x 10-16 16 is the double precision machine epsilon, then 1p* (1:) p(x)) dedull In other words, up double precision machine epsilon, p" is equal to p. Therefore, whenever É/u 2 1 - 1/3 double the algorithm simply returns the Padé approximant (that is, the scaled Newton-Schull polynomial).
The full algorithm is given in Algorithm 2. In our experiments, we never observed Algorithm 2 tak- ing more than five iterations to converge. This algorithm is implemented in full in Implementation 2.
Algorithm 2 Remez algorithm (degree 5 approximation for sign(x))
input: interval [e, u] for u > E> O.
output: Approximation p pgdd to p* arg min max [1 p(x)
pEpado sE((.u)
define Edouble = 1.11 x 10
if e/u 2
1/3 double then
Return plx) = 8 (15 (1(2////22 3(x/w)))
s/u
end if
21 (1)
1 ju,
222 te + 3 u.
En 00,
E. E-
-x
+0
while
|g
EEpa/l
Edouble do
1 k
l.
e5
aa
b
x
(k)
(x
k)
(k)
Ck
x22
(k)
(x2
12)3
Fx.
u
11+
30b
909
20
(k+1)
3b v
9913
-00
10cx
422
10cg
end while
Return p(x) axx b
+ c
25

---

Published as a conference paper at ICLR 2026
G FINITE PRECISION CONSIDERATIONS
As highlighted in Section 3.4, one must take care to implement Polar Express in finite preci- sion. In this section we outline modifications to our method to ensure stability in finite precision arithmetic.
The first issue arises when numerical round-off creates singular values that are slightly larger than our current upper bound u, Our optimal polynomials converge only when the singular values of X, are less than ue. In some cases have
Pr(u: + E) > Ur41 + E,
so over many iterations, singular value that slightly larger than u, large could grow to oo instead of converging to 1.
To fix this issue, we simply replace each polynomial x ~> (x(x) by x +> (o.111001). This safety factor corrects for round-off errors in previous iterations while only slightly changing the behavior of the polynomial on the interval [4,, uz]. though it does cause the singular values to converge to 0.999998 instead of to 1. To correct for this, the safety factor can be omitted in the final iteration. This fix is reflected in line 5 of Algorithm
The second issue was identified in Nakatsukasa & Higham (2012) and addressed in the context of polynomial iterations by Chen & Chow (2014). In general, iterative methods for polar(M aim to increase each singular value relative to the largest singular value; while ( < 0 (X)), after enough iterations, min (X,) RJ max (X) sJ 1. However, the convergence of each singular value max may not be monotonic. Over the domain [&,, wr]. our optimal polynomial pt oscillates repeatedly between 44+1 and Mt+1. so some singular values that are near may get mapped down to 4411- It so happens that this non-monotonicity. even at a single iteration can cause loss of precision. That is, problems occur if
max (z(
(r(o;) T
E Omin C Fmanx
mmax
where 0- mmin S 0z S max are singular values of X; (Nakatsukasa Higham, 2012). In the extreme case (r:)) < O, the ith singular vector will change sign, causing the method to converge to the polar factor of the wrong matrix. Unlike Newton-Schulz, unscaled Newton, or QDWH, our method is affected by this loss of precision.
To mitigate this issue, Chen & Chow (2014) propose modifying their update polynomials to enforce
a lower bound on the ratio pt(a,) This issue only occurs when & < u4; as & - uy, our optimal
polynomial approaches the Padé approximant and so P(()) 2 1 for all x E [0, uz]- We could fully solve the problem by using the Padé approximant instead of our optimal polynomial, but this would significantly slow down convergence. Instead we compromise. When le Z we/10, we find that
pa(w) 0.236. Therefore, whenever l, < ut/10 we select the update rule as though l =,/00 This change slows convergence, but only very slightly. (The choice of 10 is somewhat arbitrary. In Implementation 2, we different factor.) This fix reflected in line 4 of Algorithm 1.
The third change is copied from the original Muon implementation: normalize M by | Ir+ 10-2 instead of by M [r. As before, set U1 1. This fix is reflected in line 10 of Algorithm 1.
H ADDITIONAL EXPERIMENTAL RESULTS
In this section, we present additional experimental results.
H.1 CONVERGENCE OF POLAR EXPRESS AND ITS IMPACT ON MUON
Convergence in Frobenius Norm In Figure 8, we plot the convergence of Polar Express and three baselines as measured in the Frobenius norm. We also plot convergence in cosine similarity, which is defined with respect to the Frobenius inner product (A, B) Tr(A B). Formally, the (A.B) cosine similarity between A and B defined as TA BTF We use gradients of GPT-2 layers as test matrices. While Polar Express is designed to minimize the spectral norm error, convergence in the Frobenius norm is similar (compare with Figure 3).
26

---

Published as a conference paper at ICLR 2026
Before stabilizing
After stabilizing
0.236x
2.0
44-1
1.5
1.0
0.5
6+1
0.0
0.0
0.2
0.4
0.6
0.8
1.0
Figure 7: Effects of stabilizing the update rules with safety factor and cushioning, described in Appendix G. The blue curve is the optimal degree-5 polynomial for the interval [0.005, 1]. It is has numerical issues because it maps singular values 0.8 down to almost and maps +
U4+1 + 25c. The stabilized version is better because it ensures Px(2) z 0.236 and maps all < 1.01 to at most u1+1-
GPT-2 Gradient
GPT-2 Gradient
Relative Frobenius Error
Layer 1 attn.c_athe
Layer mlp.c
1.0
0.5
0.0
1.00
Cosine Similarity
0.75
0.50
0.25
2
4
6
8
10
2
6
8
10
Iteration
Iteration
=0=
Newton-Schulk (d 5)
Jordan
You
PolarExp, f= 10-3
Figure 8: Convergence of degree-5 polynomial methods measured in Frobenius norm and cosine similarity. Test matrices are gradients of two layers of mandomly-initilizze GPT-2 model on batch of language modeling data. Polar Express outperforms other methods.
27

---

Published as a conference paper at ICLR 2026
(In)sensitivity of Muon to Small Singular Values Figure 5 shows that using than five six iterations of Polar Express does not improve the performance of Muon. However, Figures and 8 show that five iterations is not enough for Polar Express or any other method to converge. In practice, Polas ar Express is taking steps in directions that are meaningfully different from the exact volar(M) (as computed by SVD), but still converging equally fast. One possible explanation for this observation is that Muon may not be sensitive to the convergence of small singular values of M. Intuitively, the singular vectors associated with these small singular values correspond to directions which have little effect on the output of the neural network; they may signify little more than noise in the stochastic gradients.
We now conduct an experiment to test this hypothesis. We compare three ways that a Muon like optimizer could handle the small singular values. Assume M has full rank, and partition the singular value decomposition of M into two parts
= UEV = UU1 U2]
Σ1
V
=-= = U + V.
(20)
Σg
where Σj contains the singular values larger than some threshold max and Σg contains those smaller than 7 max. where o max is the largest singular value of M. Recall that
polar (M) == UV = + U2V
(21)
obtained by mapping each singular value of M to 1. We define the truncated polar factor by mapping the larger singular values to 1 and the smaller singular values to O:
polar, = V
(22)
A third possibility is to map the small singular values to -1:
UV U.V - V2V
(23)
Note that - is in the opposite direction the Muon update. If the small singular values carry meaningful information about the loss landscape, then we expect this partly "uphill" step to hurt performance. Comparing the three update rules in Equations (21) to (23) can tell us how small singular values affect Muon.
We train GPT-2 Small using each of these three update rules with learning rate 0.05 and weight decay 0.1. We sweep three different options for the cutoff y that defines the 'small" singular values: 10-4, 10-3, and 10 2 The results are plotted in Figure 9. They show that the treatment of singular values smaller than 10- o
10-40mm mmax does not matter at all for the performance of Muon, and those smaller than 10 3 have very minor effect. Notably, even reversing the direction of the Muon step in the bottom singular subspace barely worsens performance, showing that the gradient information in this subspace not very informative. The bottom panel of Figure shows how five iterations of Polar Express (with é = 10-3) affect small singular values. Singular values greater than 10-3 are all mapped close to 1, while those smaller than 10-4 are all mapped close to O. Thus, while Polar Express does not fully converge after five iterations, it does converge in the ways that matter for Muon.
Convergence of Top Singular Values As discussed in the previous paragraph, we hypothesize that Muon may not be sensitive to the convergence of the small singular values of M when approx- imating polar(M M). Therefore, in Figure 10, we plot the convergence of Polar Express and the baselines when all singular values smaller than 10-3 3 are ignored. Specifically, if alg(M) denotes the output of an algorithm for approximating polar MM), then we compare
U,U alg(M) V.V.T
to po
where polarin-s1 (M) = U.V, = U,U,T (olar(M) V is the truncated polar factor defined above. The results show that Polar Express converges in just six iterations as measured in the relative Frobenius norm and just five iterations when measuring in cosine similarity. The other methods converge faster too, but Polar Express still outperforms them. These results may explain why the performance of Muon saturates at five or six iterations of Polar Express, as shown in Figure 5.
28

---

Published as a conference paper at ICLR 2026
5.000
4.884
Singular Value Map
4.750
yowax
Final Validation Loss
4.500
of > moa
4.250
mi (true polar)
4.000
3.750
3.669
3.589
3.587
3.612
3.588
Polar Express
0-
10
10
10
Cutoff (n)
Figure 9: Impact of small singular directions of momentum matrix optimization quality. We compare three variations of the Muon update rule. Exact Muon (green) processes the momentum M UVVT by mapping each singular value to 1: polar(M) UVT T'runcated Muon (orange) maps the larger singular values to and the smaller singular values to Reverse Muon (blue) maps the larger ones to I and the smaller ones to -1. Computations runs train GPT-2 Small on I billion tokens Web data with learning rate 0.05 and weight decay 0.1. When the cutoff that defines "large" and "small" singular values is = 10-3 all three methods perform well, showing that the small singular directions do not matter. Bottom panel shows the polynomial defined by composing five iterations of Polar Express. Five iterations just enough for singular values 2 10-3 to nearly converge.
performed in bfloat 32. All
GPT-2 Gradient
GPT-2 Gradient
Cosine Similarity Relative Probenius Eroo
Layer attn c_attn
Layer 3 mlp.c.fc
1.0
0.5
0.0
1.00
0.75
0.50
0.25
8
10
10
Iteration
Iteration
Newten-Schule (d 5)
Jordan
You
PolarExp,
Figure 10: Convergence of degree-5 polynomial methods, considering only singular values larger than max /103. Test matrices are gradients of two layers of randomly-initian inittli-iiii GPT-2 model on batch of language modeling data. Polar Express converges in just five six iterations and outperforms other methods.
29

---

Published as a conference paper at ICLR 2026
H.2 TRAINING GPT-2
Additional Metrics We report additional results from the experiment of Section 4.2. In addition to showing validation loss vs. learning rate and training step, we also report training loss vs. learning rate and training time. The results are shown in Figures 1 lla and 11b. The upper rows of each subfigure are identical to Figure and Figure 4, and are repeated here for ease of comparison.
Weight Decay As described in Section 4.3, we reran our GPT-2 training runs with weight decay of 0.1. This change had little effect on the results, shown in Figure 12.
Number of Training Tokens We also reran some of our GPT-2 training runs using 10 billion tokens of training data instead of 1 billion. As described in Section 4.3, 10 billion tokens roughly matches the Chinchilla scaling rule for GPT-2-Large and exceeds it for GPT-2-Small. Results are shown in Figure 13. Note that the top row of Figure 13a is identical to Figure 6. Polar Express still outperforms the baselines across all conditions, but the gap shrinks as the training loss con- verges.
H.3 IMAGE CLASSIFICATION
We conducted experiments on the CIFAR-IO and CIFAR-100 image classification benchmarks (Krizhevsky, 2009) using ResNet-20 and ResNet-110 110 architectures with batch normalization (He et al., 2016). We used range learning rates in the range 10 6 to 1 with a constant learning-rate 10-6
schedule, a batch size of 128, and 50 epochs of training data. We used three different random seeds for each hyperparameter setting to assess stability and variability. As a baseline, we also included AdamW and SGD with momentum (Kingma & Ba, 2015). Results are given in Figures 14 and 15. For these experiments we see that all the Muon variants performed well, matching or exceeding the training loss and validation accuracy of AdamW and sgd while also being more stable with respect to the choice of learning rate. However, we do not see a marked difference between the varieties Muon. Indeed, even Newton-Schut (degree 5) performs equally well in this context, despite being significantly less accurate than PolarExpress, Jordan or You.
Next we train a Vision Transformer (patch size 4, embedding dimension 512, depth heads, MLP dimension 512, dropout 0.1) on CIFAR-10 for 200 epochs with batch size 512 using a constant learning rate schedule. Results are shown in Figure 16. Muon with Polar Express achieved the best training and validation loss (closely followed by Jordan's 's and You's methods), However, improved loss did not entirely translate to better accuracy: both Muon and Newton-Schulz and Adam performed well in terms of validation accuracy. Overall, these experiments do not show consistent advantage for Polar Express. Further work may be beneficial to fully realize the potential benefits of Muon and to further tune Polar Express for these settings.
I INITIALIZATION FOR MATRICES WITH LARGE SPECTRAL GAPS
In Section 3, we constructed a sequence of polynomials that adapted the range of the singular values [e. ul. Assuming nothing else about the input, these polynomials are optimal since they pro- vide good approximation
across the entire interval. However, many applications, the spec- trum has large gaps; that is, there are several large outlying singular values that are well-separated from the rest. For these matrices, it is not necessary for the polynomial be accurate on the entire interval [e, u).. only on the range of the small singular values plus few other isolated points. In this section, we take advantage of this structure to accelerate our method by preprocessing the matrix to eliminate the largest singular values.
The first step is to find small intervals containing each of these large singular values. To find lower bounds, we use subspace iteration, which is generalization of the power method that approximates multiple singular values simultaneously. Fix k. the number of singular values we wish to eliminate. Letting σi 2 2 on denote the singular values of M, subspace iteration produces estimates
30

---

Published as a conference paper at ICLR 2026
4.4
reas You Fass--Pdaddip tasn--koria
Flail Valorite 1,00
3.44
Loss
4.
3.42
3.40
3.50
3.56
Learting Hate
Epurch
4.4
maas Voa
tnans- PaarEep
Jondan
4.2
3.44
4.0
Floil 1100
3.42
3.40
3.38
10-2
4000
4000
Leariing Itte
(a) Garee-Large (774M params). Best final validation losses You (r 0.02): 3.399.
muon-Jordan (Ir = 0.02): 3.398 and muon-Polarkas (Ir 0.02): 3.340.
4.6
4.
Flaal Valitate Lana
4.2
Valtciation Loon.
422
4.2
102
Leuning Race
Epoch
4.6
4.2
Find Look
4.0
6.0
3.8
200
43C
100
100)
1000
1300
Tremnin
Time
(b) GPal2-Small (124M params). Best final validation losses were adamw (Ir 0.001): 4.197. muon-Jordan (Ir 0.01): 3.639, muu (Ir 0.01): 3.629 and muon-FolarExp (Ir 0.005): 3.588.
Figure 11: Training GPT-2 on billion tokens of FineWeb data Arroc--Ouelleti et al., 2023) with- out weight decay. The label muon-<me thod> denotes Muon with iterations of <oothod> to compute polar [M). Top left: final validation loss vs. learning rate. Bottom left: final training loss vs. learning rate. Top right: validation loss vs. number iterations best learning rate. Bottom right: training loss vs. time for best learning rate.
31

---

Published as a conference paper at ICLR 2026
3.46
Flui Vildation L665
Validation
3.9
3.8
3.4
3.50
0.4
0.2
0.4
0.0
0.8
Epoch
3.40
tmann-kmias tease. PatarEp
3.40
Flad 1,00
3.42
3.40
3.38
3.50
314
2000
4000
2000
0000
Louring Tatn
Time
(a) Lar-e-Large (774M params), Best final validation losses were muon-You You (Ir = 0.02): 3.390,
muon Jorddnn (Ir = 0.02): 3.401 and muon- (Ir 0.02): 3.344.
3.72
Floal Vil0a330 La0
3.68
3,06
Viltcdation
19 TTT
3.02
3.8
3.58
0.11
0.2
0.4
13.0
Egoch
3.72
3.70
3.68
Flni 1000
3.000
3.64
3.00
Learsing
Time (4)
(b) GPT-2-Small (124M params). Best final validation losses were muon-Jordan (Ir = 0.01): 3.638, muon-You (Ir 0.005): 3.641 and muon- Poonnroolarpppp (Ir 0.005): 3.587.
Figure 12: Training GPT-2 I billion tokens of Fibe data (Aroca-Ouellette et al., 2023) with weight decay 0.1. The label muon-<methooe denotes Muon with 5 iterations of <method to compute polar(M). Top left: final validation loss vs. learning rate. Bottom left: final training loss vs. learning rate. Top right: validation loss vs. number of iterations for best learning rate. Bottom right: training loss vs. time for best learning rate.
32

---

Published as a conference paper at ICLR 2026
3.8
3,00
mass--koitee musss-Piadtaa rease Neas
J4
3.04
Fiea Valdatio Lod
3.02
13
3.00
0.2
2.08
2.00
2.04
7.00
10
0.2
0.8
1.5
Loaring Rate
Epoch
3.5
3,00
Juise
3.4
Voa
per coree com. Periee
Floil Lo55
3.00
2.08
3.00
2.08
2.00
10
10
10000 30000
20003
43000
50000
02000
Learinng Rate
(5)
(a) GPT-2-Large (774M params) with weight decay 0.1. Best final validation losses muon-Jordan (Ir 0.002): 2.921, m (Ir= 0.002); 2.919 and muon-Poloadeee PolarExp 0.002): 2.913.
3,50
3.45
Flai Vallatinn 1000
3.40
3,35
3,30
3.25
3.20
0.11
0.2
0.6
04
0.8
3.0
Ejurc
3.50
maats-Pdarfaa adazne
3.5
3.40
Plui Lo00
emt 200
3.50
3.20
20000)
21002
30000
Learning Hate
Time
(b) GPT-2-Small (124M params) without weight decay. Best final validation losses were adamw (Ir= 0.0005): 3.370, muen-Jerdan (Ir = 0.005): 3.233, muon-You (Ir 0.005): 3.234 and muon-Polarkip (Ir 0.005): 3.231.
Figure 13: Training GPT-2 on 10 billion tokens of FineWeb data (Aroca-Ouellette et al., 2023). The label muon-<met hod> denotes Muon with iterations of <method> to compute polar( M )
Top left: final validation loss vs. learning rate. Bottom left: final training loss vs. learning rate. Top right: validation loss vs. number of iterations for best learning rate. Bottom right: training loss vs. time for best learning rate.
33

---

Published as a conference paper at ICLR 2026
alucc
mnse You
ulune
moen-Nortoo mann-Fdaarog xpl-m
munu-Yoo You
mon-Jiido
0.90
0.85
Valkinko so0re
Thinins Joss
0.80
0.03
10"
0.00
Loarning rate
Lourzing rute
Figure 14: CIFARIO with a RESNET20. Shaded regions show range over three random seeds. The best validation accuracy for each method was sgd-m (Ir 0.1): 0.855 Adamw (Ir = 0.01): 0.878 mucn-You (Ir 0.001): 0.887, muon-Newton (Ir 0.001): 0.890, muon-Jordan (Ir 0.001): 0.891, muon Farxxp (Ir 0.001): 0.893.
1pki
moon monnn Poder -
splim
0.70
0.65
Validaten score
L0P
Toining loos
0.55
0.50
0.45
10-7
10-0
10-4
10
10
10
10P
10-5
10
10-5
10-2
10
10"
Learaing
Figure 15: CIFAR100 with RESNET110 100. Shaded regions show range over three random seeds. The best validation accuracy for each method was sgd-m (Ir 0.1): 0.602, Adamw (Ir 0.01): 0.643, muon-Jorde (Ir 0.001): 0.660, muon Newton (Ir 0.001): 0.663. muon-PolarExe (Ir = 0.001): 0.663, muon-You (Ir 0.001): 0.665.
ailusty
Monn-Nertun
Yon
slune
monn-Norton
munn You
Pop
0.9
0.8
Val22412 x00ee
0.7
Training Lose
0.6
10-1
0.5
0.4
0.3
10
10 Learinng rade
10
10
10
10
10
100
Louraing rato
Figure 16: CIFAR10 with a VIT. Shaded regions show range over three random seeds. The best
validation accuracy for each method sgd-m (Ir 10 10-1): 0.809, muon-PolarEwe (Ir 10-5); 0.860, Adamw (Ir = 10-3): 0.861, muon mwwn- -Jrdan (Ir = 10-5): 0.861, muon mon- (Ir 10= 0.865, Nuwn-N (Ir 10- 0.874.
34

---

Published as a conference paper at ICLR 2026
01 2 oy satisfying 0; 2 or for all i E I,...,,.""
To find upper bounds each ds, can
use the fact that |M
f = Eド= j=1 0f as follows:
of = M i
Σ
o3
M
Σ
of
S
M
P
j
(24)
==1
j=1
jp
That is, for each i E [r]].
0₂ E
M E
Setting: i k 1, the above also provides an upper bound for the tail of the spectrum, 0h4-1-.- On-
The second step is to find an odd polynomial that well-aaproxsmmtes the constant function on each of these intervals and the tail simultaneouslyy. For simplicity, we treat only the k = 1 case here. Assume that M is normalized to || = 1 and let 2 = Ø_ be the lower bound produced by subspace iteration (which reduces to the power method in this case). Then (24) gives 01 [z, 1] and 02 ........ On S /1 22. Assume that these intervals do not overlap. that is, v1 22 S 2 + N 1/v2. Then we construct the unique odd cubic polynomial p(z) = az + bx that satisfies p(v1 22) = 1 and p2z) = 1 by setting
a = 2² + v1 2²) 2
V1 2-
b
(25)
2v1 (22 11)
2V1 (2
Because p(0) = 0 and has at most one local extremum on P20, these conditions immediately guarantee that p is concave-increan increasing on [0, vI=2 so it must lie above the line x + x/ x/v1== 2 Furthermore, p is decreasing on [0i, 1). so maps 0r E [z, 1] to (p(1), 1]. By minimizing p(1) over all valid z (that is, over the interval 2 E [1/v2, 1)), one can further show that p(1) > 1/v2,, so 01 cannot be decreased very much by applying p. Thus, the largest singular value of p(M) is still at most 1, while the smaller singular values have increased by a potentially large factor of 1/v11 2?. When there is large outlying singular value, z is close to 1 and this initialization scheme makes much more progress than standard iteration of Pol a arExpeess would have.
In Figure 17, we demonstrate the benefit of using the given by (25) on a synthetic matrix whose spectrum follows power law decay. That is, 0, (M) = j-b so this matrix has large outlying singular value σ >> a2. Applying (25) costs almost as much as performing an iteration of a degree- 5 polynomial method, so for fair comparison, we count it both Newton-Schull and Polar Express, performing the extra spectrum spectrum-aware aware initialization step described in this section leads to significant speedups in convergence.
additional iteration in this plot. For
J FAST POLYNOMIAL ITERATION FOR RECTANGULAR MATRICES
In this section, we describe simple method for applying an iterative polynomial method rectan- gular matrix. For matrices with large aspect ratio, this method yields significant computational sav- ings. We emphasize that this method is applicable to any computation of the form (PT o.. o..) (X), where each pæ is an odd polynomial. Thus, it can be used to apply Cewton-Schull or Jordan's 's poly- nomials in addition to our own.
As preliminary, we first describe the baseline approach. Let X E R= with m 2 n, where a == m/n 2 1 called the aspect ratio. Any odd polynomial p degree d == 2q + can be represented as p(x) = x2?),? where h is polynomial of degree g. Thus, p(X) = X xh X). Furthermore, can be written in a factored form called Horner's rule to reduce the number of multiplications. For instance, if fvy) = a + by + cy²² + dyy", Horner's rule gives h(y) = + y (6 + y + dy)). For matrix, ((Y)= = al + Y (bl + (cl + d))). Thus for Y E Ruxn computing (YY costs about
Let Qo E m be a random matrix with orthonormal columns and define Q--1, Re+1 qr M :Q.), where qr is the QR decomposition. Subspace iteration outputs the singular values Ø1,.... of MQr, 2 2g. By the Cauchy interlacing theorem, or S Tb.
35

---

Published as a conference paper at ICLR 2026
1.0-
0.8
Spectall Erta
0.6
Method:
Newton-Sehulz
0.4
[il = 5)
PolatExp
Initialier::
Standard
Spectron-Awar
5
10
15
20
25
Figure 17: Benefits of the spectrum-award initialization scheme of Appendix Using this scheme improves convergence of both Newton-Schulz and Polar Express on a synthetic 32 x 32 matrix with ( = jj Note that we count the spectrum-aware initialization as an additional iteration.
X) costs 2mn2 + ( 1) n2 ( + 2a) n operations. This process could be repeated for each iteration Pl...., Pr. Notice that if we instead computed ((X X X
(deg(h) - 11-n""
n2 operations, and computing p(x) X
, the result would be the same but the cost would be higher.
A major drawback of this naive approach is that has a strong dependence on a, since two rectan- gular matrix multiplications must be performed in each of the T iterations. When m > n, these two multiplications dominate the cost. In Algorithm 3, we introduce simple trick that dramatically reduces this cost, using just two rectangular matrix multiplications to compute all T iterations.
Algorithm Fast Polynomial Iteration for Rectangular Matrices
input: X E Rmx" with m > 1.5n, odd polynomials pı (I) wh1 (I²))....,Pi PP(()) = (1) output: The matrix (pr o...o p.) (X).
Y x
Let Q1 =
for t 1,2.... do
R 1,1,, YQ
> 2n
Q, Q.-10e (Re
Horee's mile: deg(h.) h
end for return XQr
b mm
To why this works, define 90(x:) 1,
(e(x)) (pr .... ....oppil p1) Pt oz)((())) Pe (xq-1(x))
(26)
x
x
rqq. zq.-1 (x)) h (104)- (x))) qz) hg (x² (x)))
(27)
and rz(z) x² 91-1 1 (x)2. is clear by induction that R2 = (X). Q: qe (X). and XQr (Pr o... (X). As promised, this algorithm uses no rectangular multiplications in the for- T If each is degree d, then the total cost is 2 + 2a) n². When > 1.5 7-1- this is smaller than the naive method. We can use this criterion to select either Algorithm 3 or the baseline method runtime.! 1!
Algorithm 3 can introduce numerical errors, especially when working in a low precision format like bfloat16 We identify two sources of numerical trouble and propose remedies for each. The first is due to the ill-conditioning of X. Let X = UEV be the SVD. For large T, (pr ... (X) XQr polar(X) UVT Thus, Qr V When X has very small
"Notice that Qr + 1-2 This shows that the Polar Express polynomials also give method of computing the inverse square of PSD matrix.
36

---

Published as a conference paper at ICLR 2026
1024
2048
10a
4096
Restart Interval
Runtime (sec)
3
10-1
..*.* 6
10-2
100
10
Aspect Ratio (a m/n)
Figure 18: Effects of using Algorithm& 3 on runtime on a GPU. We run T 6 iterations of degree-5 polynomial method on matrices with various dimensions n and aspect ratios a. Restart interval is Algorithm 3, restart interval = is equivalent to the baseline (that is, not using Algorithm 3), and restart interval = 3 is an intermediate method that calls Algorithm once to the first three iterations and again to do the last three iterations for greater stability. When > 1, increasing the restart interval significantly reduces the runtime.
singular values and the floating point precision is very low, instantiating QT may be unstable. To mitigate this issue, we use a restarting strategy. Notice that the issue arises only for large T, for which (pr o... Pe)(e) 1. Limiting ourselves to T iterations improves the conditioning of Qr because (pT .o. < Thus, to compute T > 3 iterations, we begin with Xo and apply Algorithm 3 with the first three polynomials, producing X3. When then apply Algorithm again with the next three polynomials to X3. prodduiin X6. and so on. As X, approaches convergence, its conditioning improves and may longer need to restart at all. Note that restarting Algorithm after every iteration is exactly the same as the baseline method.
Second, while the matrix Y is positive definite in exact arithmetic,, numerical round-off can intro- duce spurious negative eigenvalues that cause the method to diverge to infinity. To combat this issue, we instead set
X X + 10 during the first application of Algorithm 3. (We also normalize by || + 10 3 instead of |X (e.)) In subsequent restarts of Algorithm 3. we set Y = X X as before. This akin slightly increasing each of the singular values of X, but does not change the polar factor of X. Thus, while the output will be slightly different in the early iterations, the algorithm still converges to the correct answer.
Figure 18 shows that using Algorithm 3 can significantly improve runtime on the GPU when the as- pect ratio large enough. As expected, using Algorithm 3 for many iterations significantly reduces the dependence of the runtime the aspect ratio. Running six iterations of a degree-5 polynomial method when a = 4 (as with the linear transformations in each MLP block of a transformer) we obtain almost 2x speedup, and when iterations, the trend is the same but the runtime savings are somewhat smaller.
32, obtain 5x speedup. If we restart every three
J.1 APPLICATION TO MUON
these problems can mitigated, the speed afforded by Algorithm suggests an improvement in the way Muon is applied to transformers. In sum, the idea is to replace one large matrix with a small aspect ratio by many smaller matrices with large aspect ratios and apply Algorithm 3 to all of them in parallel. Each multi-head attention layer contains four square weight matrices Wo. WK. Wy and Wo E ddd The orthogonalization step of Muon is either applied separately to these four matrices or else to [Wa Wx Wv] and Wo. since typical implementations of multi-head attention store the weights in this concatenated form. However, weight matrices to be a concatenation of many smaller linear transformations, each corresponding single attention head. If H is the number of heads, each of these smaller matrices has size
believe is natural consider each of these four
37

---

Published as a conference paper at ICLR 2026
d x it' that is, they have aspect ratio H. The gradient matrices of We WK Wv] and Wo can be reshaped into 3-tensors which each slice one of these smaller matrices. Since typical transformers like GPT-3 can have many as 96 heads, this variation of Muon has the potential to reduce the runtime.
We this idea train GPT-Small model on FineWeb1BB B. We compare four conditions:
1. The baseline approach used in the rest of this paper (not splitting (Wo WW not using Algorithm 3)
Wv] and
2. Splitting up the gradient matrices of [WQ Wx Muon to each piece, as described above
Wv] and Wo by head and applying
3. Using Algorithm 3, restarted after three iterations, on all rectangular weight matrices
4. Splitting by head and using Algorithm
We used Polar Express with weight decay of 0.1 for all conditions and swept learning rates 0.003, 0.005, 0.01. Otherwise, all hyperparameters were the same as in Section 4.2.
Our results showed that these changes had negligible effect in this setting. They did not affect the optimization quality. Compared to the baseline, splitting by heads actually reduced the final loss slightly from 3.59 to 3.55; using Algorithm 3 increased the loss very slightly, from 3.59 3.60 when not splitting by head, and from 3.55 to 3.56 when we did split. However, the runtimes of all 12 runs were nearly identical, showing that this scale, the FLOP savings of Algorithm is not beneficial. The embedding size of GPT-Small is just 768. These techniques may be more impactful when using a larger model. may also have more impact outside of deep learning, where Polar Express would be run for more than the iterations used in our experiments. We leave exploration of these settings to future work.
38


Reviews:
[{"review_id": "ZIUFGUg9F3", "reviewer": "Reviewer_caWQ", "summary": "This paper proposes a new algorithm for computing the matrix polar decomposition for the Muon optimizer used in training large neural networks. Unlike classical methods designed for high numerical accuracy, Polar Express is optimized for GPU efficiency and relatively low-precision computation (bfloat16). The authors derive the algorithm by solving a minimax optimization problem that guarantees optimal convergence in the worst case and demonstrate substantial empirical improvements over existing methods.", "questions": "see weakness section.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper introduces a theoretically grounded algorithm.\n2. The method is practical, using GPU-efficient GEMMs and is numerically stable.\n3. Experiments on GPT-2 with the Muon optimizer show consistent improvements in validation loss over existing methods.", "weaknesses": "1. It would be interesting to see how well the method works on models other than language models. \n2. Lack of detailed runtime or throughput benchmarks. The experiments show the method works better than baselines under the same number of iterations, but it is unclear whether it is faster.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HWd8vocCkq", "reviewer": "Reviewer_Xn5s", "summary": "The paper proposes a better method for computing polar factor decomposition, practical for GPU implementation. Optimality under certain conditions is proven, and improved empirical performance for implementation into the Muon algorithm for synthetic problems and GPT model training.", "questions": "Does the approximation quality depend on the dimension of the problem?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "An important and timely problem is considered with strong potential for practical impact.\n\nA novel method is proposed and analyzed mathematically.\n\nReasonable empirical experiments illustrate and support the claims.", "weaknesses": "The set of practical ML problems can be diversified a bit, e.g., by considering a vision domain.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3QgLVF5gAi", "reviewer": "Reviewer_kcTx", "summary": "This paper seeks to improve the Muon optimizer which is being used in several recent results as an alternative to Adam/AdamW for training small-sized language models. In Muon, one subroutine is calculating the polar decomposition. This is the main task handled in this work. The standard solution involves running Newton-Schulz which works well but convergence can be slow. More recent heuristics to speed it up can work but may not converge to the correct solution. Polar Express proposes to fix by adapting the polynomial at each iteration (a minimax problem is solved). This gives good initial speed of the heuristics and guarantees on convergence. On the practical side, GPT-2 models are trained to a better validation loss. The paper is well written. The technical insights are very nice and the empirical results are convincing.", "questions": "1. You may consider slightly expanding the scope of experiments. The paper is still a valuable contribution, but even if the idea is not equally effective in other cases, it will be good to point out where the gains are limited.", "limitations": "", "rating": 10, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1.  The paper replaces recent heuristic strategies which can be fast but plateau out with a provably optimal algorithm. The initial progress is rapid and has similar convergence guarantees like classical methods. I found this result very interesting and new. \n\n2. I was familiar with the Nakatsukasa/Freund paper, but the way in which the minimax paradigm is adapted here is surprising in a good way, giving Thm 3.1 and 3.3. I feel that the idea is quite novel and can be useful beyond the Muon use case. \n\n3. The experimental results back up the main claims. Fig. 3 is sufficient validation that the idea in this paper deserves consideration. \n\n4. What is also nice is that this is not just a algorithmic results paper. The authors discuss numerical instability by adjusting the polynomials, showing the adjustments needed for practical gains on modern GPUs. Same with precomputing coefficients etc. Excellent balance of technical findings with practicality.", "weaknesses": "1. Relatively minor. It will be good to show how much slower is Zolo-PD using newer QR implementations. \n2. Also minor. Experiments are only on GPT-2 with Muon. If Muon is effective, why not expand the scope of experiments to other architectures and even finetuning?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nn7ke876m2", "reviewer": "Reviewer_WacA", "summary": "This paper proposes a new method for orthogonalizing a matrix, with applications in the Muon optimizer. The method, PolarExpress, greedily computes polynomial approximations of the sign function in a residual manner using the Remez algorithm, and then applies these approximations during training to approximately orthogonalize the momentum buffer in Muon. The authors prove results about the convergence of PolarExpress and its optimality with respect to the L_\\infty norm to the sign function, and empirically show that PolarExpress improves upon existing orthogonalization methods used in current Muon implementations.", "questions": "See above.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I will preface this review by saying that I am not an expert in numerical methods. However, the presented algorithm seems novel and interesting, and the there is theory showing that PolarExpress is guaranteed to converge at least as fast as the canonical \"NewtonShultz5\"  method. This paper also has some other nice components, including\n- Theoretical justification for using a greedy residual polynomial construction\n- Empirical results on toy matrices showing that the method converges faster than existing baselines (Figure 3)\n- Results showing the effect of varying the lower bound guess $\\ell$ of the smallest singular value in the gradient matrix\n- Empirical GPT training experiments showing that PolarExpress achieves better initial convergence over prior Muon implementations.", "weaknesses": "- There seems to be a missing \"key\" baseline, which is running Muon with the actual `polar' function and computing the SVD online. Although this is obviously impractical, it gives a sense of how far PolarExpress and these other methods are from \"true\" Muon. Can you run an experiment showing how much benefit there is from getting the optimal solution in real-world GPT training?\n- Likewise, does minimizing the L_\\infty norm of the sign function maximize the directional similarity to the actual \"true\" Muon update? Can you generate a plot of the cosine distance of PolarExpress's update (and other baselines) to the SVD update to get a sense of how close these approximations are?\n- The empirical experiments only go out to 1B tokens. 125M parameters @ 1B tokens is way below Chinchilla optimal, and my own experiments in the past have suggested that Muon's initial performance gap diminishes in the \"overtrained regime.\" If you have time, can you run an experiment past Chinchilla (125M parameters @ 10B tokens should be easily doable in a day on the 4 H100s used in the paper) to see how much of a gap there is?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)