The Polar Express Algorithm

10 October 2026 · Download PDF

In this note we look at how the sign function can be approximated by compositions of polynomials, as in the Polar Express algorithm [1]. This simple idea has a number of applications in modern machine learning and optimization, as it can be used to approximate matrix functions using only matrix-matrix products, which are very efficient to compute on modern hardware such as GPUs. Applied to the singular values of a matrix, the approximation of the sign function gives the polar factor $𝑈 𝑉^{⊤}$ . This is, e.g., the basis of the update of the Muon optimizer [2], which would otherwise require an SVD, an expensive computation that is not very GPU friendly. Muon is explained as steepest descent in the spectral norm in Which Way Is Down? It is furthermore possible to approximate the ReLU function in a similar way, which can be used to project matrices onto the cone of positive semidefinite matrices.

1. The Singular Value Transformation

The singular value transformation is defined for any odd function $𝑓$ , as

𝑓 (𝑀) ≔ 𝑈 𝑓 (Σ) 𝑉^{⊤}, 𝑀 = 𝑈 Σ 𝑉^{⊤} .

(1)

When $𝑀$ is square and symmetric, the singular value transformation of an odd $𝑓$ coincides with the matrix function $𝑓 (𝑀) = 𝑄 𝑓 (Λ) 𝑄^{⊤}$ , which applies $𝑓$ to the eigenvalues. Writing $𝑀 = 𝑄 Λ 𝑄^{⊤} = (𝑄 sign (Λ)) | Λ | 𝑄^{⊤}$ shows that $𝑈 = 𝑄 sign (Λ)$ , $Σ = | Λ |$ and $𝑉 = 𝑄$ , and since $𝑓$ is odd,

𝑓 (𝑀) = 𝑄 𝑓 (Λ) 𝑄^{⊤} = 𝑄 sign (Λ) 𝑓 (| Λ |) 𝑄^{⊤} = 𝑈 𝑓 (Σ) 𝑉^{⊤},

(2)

where the sign function is defined as the odd function

sign (𝑥) : = {\begin{cases} 1, & 𝑥 > 0, \\ 0, & 𝑥 = 0, \\ - 1, & 𝑥 < 0, \end{cases}

(3)

2. Approximating the Polar Normalization

Polar normalization corresponds to the singular value transformation when $𝑓$ itself is the sign function, i.e.

polar (𝑀) = 𝑈 diag (sign (𝜎_{1}), \dots, sign (𝜎_{𝑘})) 𝑉^{⊤} = 𝑈 𝑉^{⊤},

(4)

where the last equality takes $𝑈$ and $𝑉$ from the compact SVD, which keeps only the singular vectors of the nonzero singular values. Zero singular values stay zero, as $sign (0) = 0$ .

Computed directly, the polar factor requires an SVD, which is computationally expensive and not very GPU friendly (it is primarily memory bound rather than compute bound). The Polar Express algorithm instead approximates the matrix sign function using just matrix-matrix products [1]. The key idea is rather simple: Approximate the function $𝑓$ using a composition of polynomials, e.g.,

𝑓 (𝑥) \approx 𝑝_{𝑘} (𝑝_{𝑘 - 1} (\dots 𝑝_{2} (𝑝_{1} (𝑥)) \dots)) .

(5)

Since the above is a composition of polynomials, it can be evaluated using just additions and multiplications.

Example 1: Newton–Schulz

Example 1 (Newton–Schulz) The idea of approximating functions using repeated polynomial application is not new. An example is the Newton–Schulz algorithm which is based on repeatedly applying the same polynomial, i.e. $𝑝_{1} = 𝑝_{2} = \dots = 𝑝_{𝑘}$ , in order to approximate the sign function. For order 3 and 5 the polynomials look as follows

𝑝 (𝑥) = \frac{3}{2} 𝑥 - \frac{1}{2} 𝑥^{3}, 𝑝 (𝑥) = (15 𝑥 - 10 𝑥^{3} + 3 𝑥^{5}) / 8 .

(6)

While the Newton–Schulz method converges for all $𝑥 \in [- 1, 1]$ , the convergence is only fast for values close to $\pm 1$ . Conversely, for values that are close to zero the convergence is slow. The slow convergence is exactly what the Polar Express method explained shortly aims to fix. The graph shows the cubic and quintic polynomials composed $𝑘$ times, together with the sign function (dotted). The shaded interval $[- 1, 1]$ is where they approximate it.

𝑘

A note here is that for non-square matrices we cannot simply do $𝑀^{𝑛}$ . Instead, what is done is that $𝑀^{2 𝑛 + 1} = {𝑀 (𝑀^{⊤} 𝑀)}^{𝑛} = {(𝑀 𝑀^{⊤})}^{𝑛} 𝑀$ .

3. The Polar Express Algorithm

The basic question that the Polar Express algorithm aims to solve is as follows: Given a set of $𝑇$ polynomials ${𝑝_{𝑖} (𝑥) = 𝑎_{1}^{𝑖} 𝑥 + 𝑎_{3}^{𝑖} 𝑥^{3} + 𝑎_{5}^{𝑖} 𝑥^{5}}_{𝑖 = 1}^{𝑇}$ , can we find an optimal set of coefficients that minimizes the maximum error over a given interval $[ℓ, 𝑢]$ ? Define their composition as $𝑝 = 𝑝_{𝑇} ⚬ 𝑝_{𝑇 - 1} ⚬ \dots ⚬ 𝑝_{1}$ . We then solve

\begin{aligned} (𝑝_{1}^{⋆}, \dots, 𝑝_{𝑇}^{⋆}) & = \underset{𝑝_{1}, \dots, 𝑝_{𝑇}}{arg min} max_{\begin{matrix} 𝑀 \in ℝ^{𝑚 \times 𝑛} \\ 𝜎 (𝑀) \subset [ℓ, 𝑢] \end{matrix}} ‖ polar (𝑀) - 𝑝 (𝑀) ‖_{2} \\ subject to & 𝑝_{𝑡} \in ℝ_{5}^{odd} [𝑥], 𝑡 = 1, \dots, 𝑇 . \end{aligned}

(7)

Here $ℝ_{𝑑}^{odd} [𝑥]$ denotes the space of real odd polynomials of degree at most $𝑑$ . Since $𝑝$ is odd, $𝑝 (𝑀) = 𝑈 𝑝 (Σ) 𝑉^{⊤}$ , so the error is ${max}_{𝑖} | 1 - 𝑝 (𝜎_{𝑖}) |$ , and the problem is the scalar one of approximating $1$ on $[ℓ, 𝑢]$ , or by oddness the sign function on $[- 𝑢, - ℓ] \cup [ℓ, 𝑢]$ .

In [1] they show that the above can be solved by a greedy approach, with a caveat that the error bound only holds if the smallest singular value is at least $ℓ$ , which we do not know in advance. Smaller singular values still converge, but more slowly. A comment here is that the main focus in the paper is on bfloat16. Here the spacing between representable numbers at 1 is $2^{- 7} = 0.0078125$ , so singular values much smaller than this are at the level of rounding errors. They therefore suggest setting $ℓ = 10^{- 3}$ , and $𝑢 = 1$ after scaling $𝑀$ by its Frobenius norm so that $𝜎_{1} \leq 1$ . In addition, they provide a few modifications in order to make the method numerically robust. These additions are skipped here, but can be found in the original paper.

Interactive figure 1: Polar Express against Newton–Schulz

Interactive figure 1 (Polar Express against Newton–Schulz) The graph compares the first $𝑘$ Polar Express polynomials for $[ℓ, 𝑢] = [10^{- 3}, 1]$ , composed, with the Newton–Schulz quintic composed $𝑘$ times, drawn as in Example 1. Newton–Schulz approaches $\pm 1$ slowly from below, while Polar Express first lifts the small values and then oscillates around $\pm 1$ , with every step shrinking the oscillation. The oscillation is the price for the smallest worst-case error. By the equioscillation theorem, the polynomial with the smallest maximum error on an interval reaches that error several times with alternating sign [3], so each optimal step has to oscillate around $\pm 1$ . Newton–Schulz instead matches the sign function closely at $\pm 1$ , where $𝑝^{'} (1) = 0$ , which makes it accurate close to $\pm 1$ but slow close to zero. After a few steps the oscillations become too dense to be told apart, and soon after that too small to be seen. The readout gives both compositions at $ℓ$ , which is too close to zero to be seen in the graph.

𝑘

Values at $𝜎 = ℓ$ , Newton–Schulz and Polar Express:

Interactive figure 2: Polar Express on a matrix

Interactive figure 2 (Polar Express on a matrix) The same polynomials applied to a $2 \times 2$ matrix $𝑀$ . A matrix maps the unit circle (dotted) to an ellipse whose semi-axes are its singular values. The dashed ellipse belongs to $𝑀$ , scaled so that $𝜎_{1} \leq 1$ ,¹¹ The scaling is a requirement: the polynomials are designed for singular values in $[ℓ, 1]$ , the positive half of $[- 1, 1]$ in the previous figure, and larger singular values need not converge. Dividing $𝑀$ by its Frobenius norm guarantees $𝜎_{1} \leq 1$ , as the Frobenius norm is at least the largest singular value. and its columns (the two points) can be dragged. The solid ellipse and arrows show the matrix after $𝑘$ steps. Each step applies the next polynomial to the singular values, so the ellipse turns into the unit circle, which is the ellipse of $𝑈 𝑉^{⊤}$ .

𝑘

After $𝑘$ steps:

4. Projection onto the semidefinite cone

The positive semidefinite cone is the set of all real symmetric $𝑛 \times 𝑛$ positive semidefinite matrices. That is

𝕊_{+}^{𝑛} : = {𝑋 \in ℝ^{𝑛 \times 𝑛} | 𝑋 = 𝑋^{⊤}, 𝑋 ⪰ 0} .

(8)

Semidefinite programming (SDP) problems are optimization problems for which matrix variables are constrained to be positive semidefinite, a special case of the conic programs in Conic Optimization. As such, algorithms that aim to solve SDPs need to ensure that the solution satisfies the positive semidefinite constraint. In interior point methods this is done by only taking steps for which the next iterate remains inside the cone, while first-order methods most often rely on projecting back onto the set after taking a step, i.e., for a symmetric $𝑋 = 𝑈 Λ 𝑈^{⊤}$ ,

Π_{𝕊_{+}^{𝑛}} (𝑋) = \underset{𝑌 \in 𝕊_{+}^{𝑛}}{arg min} ‖ 𝑌 - 𝑋 ‖_{𝐹} = 𝑈 diag (ReLU (𝜆_{1}), \dots, ReLU (𝜆_{𝑛})) 𝑈^{⊤} .

(9)

While the idea of projecting back onto the semidefinite cone using just matrix-matrix products is not new (a similar fixed-point approach was suggested in [4]), the main takeaway is that the ideas of Polar Express can be used to find a faster converging method [5]. An initial idea would be to directly find polynomials that approximate the ReLU function, i.e., solving

\begin{aligned} inf_{𝑓_{1}, \dots, 𝑓_{𝑇}} max_{𝑥 \in [- 1, 1]} | 𝑓_{𝑇} ⚬ 𝑓_{𝑇 - 1} ⚬ \dots ⚬ 𝑓_{1} (𝑥) - ReLU (𝑥) | \\ subject to & 𝑓_{𝑡} \in ℝ_{𝑑_{𝑡}} [𝑥], 𝑡 = 1, \dots, 𝑇, \end{aligned}

(10)

where $𝑇$ is the total number of polynomials in the composition, $𝑑_{𝑡}$ is the degree of the $𝑡$ ’th polynomial, and the degrees are fixed in advance as part of a “matrix-matrix” multiplication budget. Unfortunately, directly solving the above is challenging, and instead a two-stage approach is proposed. The idea is to realize that we can write the ReLU function as

ReLU (𝑥) = \frac{1}{2} 𝑥 (1 + sign (𝑥)) .

(11)

The first step is to compute optimal coefficients for approximating the sign function using the Polar Express approach. The second step is then to refine the coefficients of the polynomials found in the first step by minimizing the error to the ReLU function on the entire interval $[- 1, 1]$ , i.e., minimizing the following loss function

𝐿 (𝑓_{𝑇}, \dots, 𝑓_{1}) : = max_{𝑥 \in [- 1, 1]} | \frac{1}{2} 𝑥 (1 + 𝑓_{𝑇} ⚬ 𝑓_{𝑇 - 1} ⚬ \dots ⚬ 𝑓_{1} (𝑥)) - ReLU (𝑥) | .

(12)

Interactive figure 3: From the sign function to the ReLU

Interactive figure 3 (From the sign function to the ReLU) The graph shows $\frac{1}{2} 𝑥 (1 + 𝑝 (𝑥))$ , with $𝑝$ the Polar Express composition from the previous figure, together with the ReLU (dotted). The shaded interval $[- 1, 1]$ is where it approximates the ReLU. The oscillation is inherited from $𝑝$ , and the factor $𝑥$ damps the error of $𝑝$ where it is largest, close to zero.

𝑘

Interactive figure 4: Projecting a matrix onto the semidefinite cone

Interactive figure 4 (Projecting a matrix onto the semidefinite cone) The same idea for a symmetric $2 \times 2$ matrix $𝐴$ with eigenvalues $𝜆_{1}, 𝜆_{2}$ and an eigenvector at angle $𝜑$ . For each unit vector $𝑢$ the graph draws a point at distance $| 𝑢^{⊤} 𝐴 𝑢 |$ in the direction of $𝑢$ , solid where $𝑢^{⊤} 𝐴 𝑢 > 0$ and dashed where it is negative, so the lobes point along the eigenvectors and have lengths $| 𝜆_{1} |$ and $| 𝜆_{2} |$ . $𝐴$ is positive semidefinite exactly when nothing is dashed.²² If $𝐴$ is already positive semidefinite, the curve first moves away from $𝐴$ and then returns. The first polynomials are tuned for singular values near $ℓ = 10^{- 3}$ and temporarily push values close to $1$ far off. The faint lobes are $𝐴$ , the dotted curve its exact projection, and the bold curve the approximation $\frac{1}{2} 𝐴 (𝐼 + 𝑝 (𝐴))$ after $𝑘$ steps. As $𝑘$ grows, the dashed lobes vanish and the solid ones grow slightly, since the negative eigenvalue no longer subtracts from them.

𝜆_{1}

0.8

𝜆_{2}

-0.5

𝜑

𝑘

After $𝑘$ steps:True projection:

Cite this post

bibtex


@misc{paltorp2026polarexpress,
  author = {Paltorp, Mikkel},
  title = {{The Polar Express Algorithm}},
  year = {2026},
  month = oct,
  url = {https://mipals.github.io/pdf/2026/01-03-polar-express/polar-express.pdf},
  note = {Blog post}
}

Bibliography

[1] N. Amsel, D. Persson, C. Musco, and R. M. Gower, “The Polar Express: Optimal Matrix Sign Methods and Their Application to the Muon Algorithm.” [Online]. Available: https://arxiv.org/abs/2505.16932
[2] K. Jordan et al., “Muon: An optimizer for hidden layers in neural networks.” [Online]. Available: https://kellerjordan.github.io/posts/muon/
[3] L. N. Trefethen, Approximation Theory and Approximation Practice, Extended Edition. Society for Industrial, Applied Mathematics, 2019. doi: 10.1137/1.9781611975949.
[4] J. B. Francisco and D. S. Gonçalves, “A fixed-point method for approximate projection onto the positive semidefinite cone,” Linear Algebra and its Applications, vol. 523, pp. 59–78, 2017, doi: 10.1016/j.laa.2017.02.014.
[5] S. Kang, H. Han, A. Groudiev, and H. Yang, “Factorization-free Orthogonal Projection onto the Positive Semidefinite Cone with Composite Polynomial Filtering.” [Online]. Available: https://arxiv.org/abs/2507.09165