Newton–Schulz Iterations

• Newton–Schulz iterations are fixed-point matrix recurrences that replace factorizations with GEMM operations to approximate inverses, square roots, and polar factors.
• They are employed in modern optimization for efficient orthogonalization and retractions on the Stiefel manifold, reducing reliance on SVD or QR routines.
• Researchers tune finite-step variants and polynomial designs to balance convergence accuracy and computational efficiency in applications like Muon optimizers and group synchronization.

Searching arXiv for papers on Newton–Schulz iterations, orthogonalization, Muon, and Stiefel applications. Searching arXiv for papers on Newton–Schulz iterations, orthogonalization, Muon, and Stiefel applications. Newton–Schulz iterations are fixed-point Newton-type matrix recurrences that replace factorizations and solves by matrix–matrix multiplications. In classical numerical linear algebra they are used for matrix inversion and related matrix functions; in contemporary optimization they are used to approximate inverse square roots, matrix sign functions, and polar factors, thereby enabling efficient orthogonalization, retractions on the Stiefel manifold, and semi-orthogonal update rules in Muon-type optimizers. Their contemporary importance derives from two linked properties: rapid local convergence under spectral scaling, and an implementation profile dominated by GEMMs rather than SVD or QR routines (Stotsky, 2022, Stotsky, 2020, Grishina et al., 12 Jun 2025, Shulgin et al., 22 Oct 2025).

1. Canonical formulations and normalization conditions

The classical second-order Newton–Schulz iteration for the inverse of an invertible matrix $A$ is

$X_{k+1} = X_k(2I - AX_k),$

with convergence under a contraction condition such as $\rho(I-AX_0)<1$ or $\|I-AX_0\|<1$ (Stotsky, 2022, Stotsky, 2020, Maity, 29 Sep 2025). In this form the residual obeys the quadratic recursion

$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$

and higher-order “hyperpower” variants generalize this to

$$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$

which gives order-$n$ local convergence when the initial residual is spectrally contractive (Stotsky, 2022).

For symmetric positive definite matrices, Newton–Schulz is also used to approximate the inverse square root. A standard form is

$X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$

or an equivalent implementation such as

$X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$

with scaling chosen so that the spectrum is mapped into a neighborhood of the identity, typically $(0,1]$ after normalization (Grishina et al., 12 Jun 2025, Shulgin et al., 22 Oct 2025, Tang et al., 25 Jun 2026). Closely related coupled recurrences track both a square root and an inverse square root,

$X_{k+1} = X_k(2I - AX_k),$0

again under appropriate scaling (Maity, 29 Sep 2025).

The orthogonalization use-case follows from polar decomposition. For $X_{k+1} = X_k(2I - AX_k),$1 with full column rank,

$X_{k+1} = X_k(2I - AX_k),$2

is the polar factor, and Newton–Schulz can approximate $X_{k+1} = X_k(2I - AX_k),$3 or act directly on $X_{k+1} = X_k(2I - AX_k),$4 through the “direct NS” polar iteration

$X_{k+1} = X_k(2I - AX_k),$5

with $X_{k+1} = X_k(2I - AX_k),$6 chosen so that $X_{k+1} = X_k(2I - AX_k),$7 (Peng et al., 7 Apr 2026, Shulgin et al., 22 Oct 2025). Rectangular row-oriented and column-oriented forms differ only in whether the smaller Gram is $X_{k+1} = X_k(2I - AX_k),$8 or $X_{k+1} = X_k(2I - AX_k),$9. This direct orthogonalization viewpoint is the form most visible in modern optimization systems.

Historically, the classical second-order Newton–Schulz iteration dates back at least to Schulz (1933), and later work systematized higher-order forms, factorized hyperpower schemes, and integrations with Richardson iteration for parameter estimation and generalized inverse computation (Stotsky, 2020, Stotsky, 2022).

2. Singular-value action, polar geometry, and manifold structure

Newton–Schulz orthogonalization is most naturally understood through the singular value decomposition. If

$\rho(I-AX_0)<1$0

then an odd-polynomial orthogonalization step

$\rho(I-AX_0)<1$1

acts diagonally on singular values: $\rho(I-AX_0)<1$2 Thus the iteration preserves left and right singular vectors while transforming singular values toward unity (Grishina et al., 12 Jun 2025, Huang, 29 May 2026). In Muon-style iterations, this is realized through low-degree matrix polynomials in a Gram matrix rather than by explicit computation of a polar factor.

The geometric target is the polar factor. For a thin SVD $\rho(I-AX_0)<1$3, the polar factor is

$\rho(I-AX_0)<1$4

In the Muon analysis it is a partial isometry with $\rho(I-AX_0)<1$5, and it satisfies

$\rho(I-AX_0)<1$6

which identifies it as the maximally aligned direction in the nuclear–operator norm geometry used for convergence analysis (Kim et al., 27 Jan 2026). This is why orthogonalization by polar approximation is not merely a numerical convenience: it is the geometry-aware descent direction underlying the optimizer.

For polynomial Newton–Schulz updates of the form

$\rho(I-AX_0)<1$7

the singular subspaces are preserved. If $\rho(I-AX_0)<1$8, then

$\rho(I-AX_0)<1$9

so the $\|I-AX_0\|<1$0 factors are unchanged, the range of $\|I-AX_0\|<1$1 is preserved, and $\|I-AX_0\|<1$2 is invariant under the iteration (Kim et al., 27 Jan 2026). This invariance is central to finite-step analyses of practical orthogonalization routines.

On the Stiefel manifold,

$\|I-AX_0\|<1$3

the tangent and normal spaces at $\|I-AX_0\|<1$4 are

$\|I-AX_0\|<1$5

The order-$\|I-AX_0\|<1$6 Newton–Schulz correction

$\|I-AX_0\|<1$7

lies in the normal space, and the paper on retraction-free second-order methods interprets Newton–Schulz as motion along the normal bundle of a layered manifold surrounding the Stiefel constraint set (Xiong et al., 4 May 2026). This geometric interpretation distinguishes Newton–Schulz from generic fixed-point orthogonalizers: it reduces infeasibility by an explicitly normal-space displacement while leaving tangent-space objective reduction to a separate mechanism.

3. Convergence orders, residual recursions, and inexactness models

The classical convergence mechanism is residual squaring. For inverse iteration one has $\|I-AX_0\|<1$8, and for higher-order hyperpower variants $\|I-AX_0\|<1$9 (Stotsky, 2022, Stotsky, 2020). For polar orthogonalization, the same principle is expressed through contraction of the Gram error or orthogonality residual once the spectrum has been normalized into the admissible region.

A modern finite-step analysis appears in the study of Muon with Newton–Schulz. There the momentum matrix $F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$0 is first scaled by

$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$1

which guarantees $F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$2 and places the spectrum of $F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$3 in $F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$4 (Kim et al., 27 Jan 2026). The degree-$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$5 Newton–Schulz polynomial is

$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$6

and the induced orthogonality residual contracts according to

$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$7

Hence,

$F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$8

and the polar approximation error $F_k = I - X_kA,\qquad F_{k+1}=F_k^2,$9 yields a multiplicative convergence factor

$$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$0

The paper proves that $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$1 doubly exponentially in $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$2, so Muon with a finite number of Newton–Schulz steps matches the SVD-polar idealization up to a rapidly vanishing constant factor (Kim et al., 27 Jan 2026). This is a formal justification of the empirical observation that “a few NS steps” can already behave like exact polar orthogonalization in iteration count.

The same paper adopts a nonconvex stationarity measure in nuclear norm,

$$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$3

and proves that Muon with Newton–Schulz converges to a stationary point at the same rate as the SVD-polar idealization, up to the constant $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$4 (Kim et al., 27 Jan 2026). It also proves that Muon removes the typical square-root-of-rank loss compared to SGD with momentum under the same nuclear-norm stationarity metric.

A complementary perspective comes from the analysis of the inexact Muon update. There, the exact linear minimization oracle solution $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$5 is replaced by an approximate $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$6 with additive error

$$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$7

The deterministic and stochastic bounds show that lower oracle precision requires a smaller step size and a larger momentum parameter; in particular, larger inexactness degrades the convergence bound through factors such as $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$8, while the optimal schedules satisfy $$G_k = \Bigl(I + F_{k-1} + F_{k-1}^2 + \cdots + F_{k-1}^{\,n-1}\Bigr)G_{k-1},\qquad F_k = F_{k-1}^n,$$9 and $n$0 in the stochastic setting (Shulgin et al., 22 Oct 2025). This elevates the number of Newton–Schulz steps from a low-level implementation choice to a hyperparameter that couples directly to learning-rate and momentum schedules.

Outside optimization, NS-RGS for orthogonal group synchronization establishes that blockwise Newton–Schulz retraction achieves linear convergence to the target solution up to near-optimal statistical noise levels, provided spectral initialization places the iterate in the appropriate local region and the inexact retraction error stays below a prescribed bound (Peng et al., 7 Apr 2026). On the Stiefel manifold, the retraction-free SOL and SOL-sym methods obtain local quadratic convergence in the exact case, and superlinear convergence of order $n$1 under Eisenstat–Walker forcing for inexact linear solves, with the order-$n$2 Newton–Schulz normal component providing quadratic reduction of the infeasibility $n$3 (Xiong et al., 4 May 2026).

4. Role in optimization and synchronization

In Muon, Newton–Schulz is used to orthogonalize the momentum matrix rather than the raw gradient. The update analyzed in the convergence paper is

$n$4

where $n$5 is obtained by scaling $n$6, applying $n$7 Newton–Schulz steps, and using the result as an approximation to $n$8 (Kim et al., 27 Jan 2026). The theoretical significance is that the polar direction is optimally aligned in operator–nuclear geometry, and the practical significance is that Newton–Schulz makes this geometry usable without an SVD at every layer and step.

The group-synchronization setting places Newton–Schulz in a Riemannian optimization loop over $n$9. The objective is

$X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$0

and the NS-RGS method updates each block by a Riemannian gradient step followed by blockwise Newton–Schulz sign/polar retraction (Peng et al., 7 Apr 2026). The convergence analysis uses spectral initialization and refined leave-one-out arguments to show linear contraction to the target up to statistical error. Empirically, the method attains accuracy comparable to the generalized power method while achieving nearly a $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$1 speedup, and on Lucy 3D global alignment the reported runtime is $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$2s for NS-RGS versus $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$3s for GPM, with relative error approximately $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$4 for RTR, GPM, and NS-RGS (Peng et al., 7 Apr 2026).

On the Stiefel manifold, Newton–Schulz appears in a different role: not as a retraction applied after a tangent step, but as the normal component of a retraction-free second-order landing scheme. The update is

$X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$5

where $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$6 solves a modified Newton equation designed to account for the normal-induced perturbation of the gradient (Xiong et al., 4 May 2026). The resulting SOL and SOL-sym methods are proved to enjoy local quadratic convergence, or superlinear convergence for the inexact variant, and numerical experiments on orthogonal Procrustes, PCA, and real-data ICA show better performance than the compared methods.

These applications reveal a common pattern. Newton–Schulz is not only a generic orthogonalizer; it is a mechanism for embedding matrix-function structure into an optimization algorithm without factorization-heavy projections. In Muon, the structure is a geometry-aware update rule; in group synchronization, it is blockwise Riemannian retraction; in Stiefel optimization, it is an explicit normal-space feasibility correction.

5. Polynomial design, factorization strategies, and finite-step variants

A large part of the modern literature concerns not the existence of Newton–Schulz iterations, but the design of better finite-step polynomials. The Chebyshev-optimized Newton–Schulz method CANS replaces fixed coefficients by minimax odd polynomials on a prescribed singular-value interval $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$7. For degree $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$8, the paper derives a closed form via Chebyshev alternance; for higher degrees it uses a Remez algorithm to compute optimal coefficients (Grishina et al., 12 Jun 2025). This preserves the all-matmul pipeline while improving convergence per GEMM, and the paper reports faster convergence than classical Newton–Schulz on random $X_{k+1} = \tfrac12 X_k(3I - AX_k^2),$9 matrices, better test loss than Muon’s original polynomial at the same matmul budget in NanoGPT training, and the lowest per-epoch time among the compared Stiefel retractions in Wide ResNet experiments (Grishina et al., 12 Jun 2025).

Turbo-Muon modifies the first stage rather than the polynomial coefficients. Its “Almost Orthogonal Layer” preconditioner rescales columns by

$X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$0

so that the initial Gram matrix is better conditioned before Newton–Schulz begins (Boissin et al., 4 Dec 2025). Because the first Gram computation is reused, the overhead is negligible relative to the baseline NS routine. The paper reports that this allows the removal of one iteration out of the usual five without degrading approximation quality, yields up to a $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$1 speedup in the Newton–Schulz approximation, and improves end-to-end training runtime by $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$2–$X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$3 in realistic training scenarios (Boissin et al., 4 Dec 2025).

Another line of work changes the granularity of the operator. HiMuon partitions each momentum-gradient matrix into $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$4 tiles, applies the same finite Newton–Schulz map independently to each tile, and reassembles the output. For fixed finite $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$5, this is explicitly not a convergent approximation to the full-matrix update; it is a local matrix-function map that preserves spectral interactions within tiles and discards them across tile boundaries (Tang et al., 25 Jun 2026). The leading work changes from

$X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$6

for the full-matrix map to

$X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$7

for the tiled map, enabling tile-size-dependent kernels, cross-layer batching, memory-bounded chunking, and runtime tile-size schedules (Tang et al., 25 Jun 2026).

A separate finite-step design question is addressed by the relaxed cubic schedule proposed to study how much orthogonalization Muon actually needs. Instead of seeking a more accurate polar solver, the paper derives a five-step cubic schedule that moves protected singular values into a relaxed target band $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$8 under bfloat16 constraints. Each cubic step uses two dominant multiplications rather than the three required by a quintic step, so cubic5 uses ten dominant matrix multiplications versus fifteen for five quintic iterations (Huang, 29 May 2026). The central empirical conclusion is that training quality is not governed monotonically by polar-decomposition accuracy: truncated Polar Express, Muon-Jordan quintic, cubic5, and an explicit FP32 SVD polar factor can reach nearly indistinguishable final loss on GPT-2 Small, and cubic5 matches Muon-Jordan within about $X_{k+1} = \tfrac12 X_k(3I - X_kAX_k),$9 validation loss on hybrid MoE/Mamba models from one billion to four billion parameters (Huang, 29 May 2026).

The earlier matrix-inversion literature developed analogous factorization ideas in a different idiom. Unified factorization schemes express a high-order Newton–Schulz series with order $(0,1]$0 through structured groupings such as

$(0,1]$1

reducing the number of matrix–matrix multiplications per cycle and defining an efficiency index

$(0,1]$2

for the factorized scheme (Stotsky, 2020). This older factorization viewpoint and the newer polynomial-design viewpoint serve the same end: improved contraction per unit of matrix-multiplication budget.

6. Computational profile, hardware alignment, and limitations

The practical appeal of Newton–Schulz is that its dominant kernels are GEMMs. NS-RGS emphasizes that the method uses only dense matrix multiplications and additions, which align well with GPUs, TPUs, and tensor cores, while SVD and QR have more sequential dependencies and are typically memory- and latency-bound (Peng et al., 7 Apr 2026). In Muon-style applications, a direct polar Newton–Schulz step requires roughly one Gram computation and one multiplication by a small dense polynomial in that Gram, so total cost grows linearly with the number of inner steps (Shulgin et al., 22 Oct 2025).

The Muon convergence paper makes the wall-clock comparison explicit for a layer of shape $(0,1]$3. Exact thin SVD has

$(0,1]$4

floating-point complexity, while $(0,1]$5 degree-$(0,1]$6 Newton–Schulz steps implemented by Horner’s rule have

$(0,1]$7

complexity (Kim et al., 27 Jan 2026). Since GEMM utilization on GPUs is much higher than SVD throughput, the paper argues that a few low-degree Newton–Schulz steps achieve near-SVD behavior at a fraction of the wall-clock cost. Similar arguments underlie the results of Turbo-Muon, CANS, and HiMuon, each of which improves either the contraction per step or the execution efficiency per step (Boissin et al., 4 Dec 2025, Grishina et al., 12 Jun 2025, Tang et al., 25 Jun 2026).

Several limitations recur across the literature. First, convergence depends on scaling or on staying within a local basin. NS-RGS states the classical sufficient condition

$(0,1]$8

for quadratic convergence of the sign iteration, and notes that if $(0,1]$9 is far from orthogonal one may need scaling or preconditioning to enter the basin (Peng et al., 7 Apr 2026). The inexact Muon analysis likewise treats orthogonalization error as an explicit optimization parameter because finite-step Newton–Schulz can be materially suboptimal when the number of steps is too small (Shulgin et al., 22 Oct 2025).

Second, more exact polar approximation is not always better for downstream training. The cubic5 study explicitly finds that training quality is not governed monotonically by polar-decomposition accuracy, and that an explicit FP32 SVD polar factor does not outperform strong approximate Newton–Schulz variants in the tested GPT-2 Small regime (Huang, 29 May 2026). This suggests that the relevant object in optimizer design is not only approximation error to the exact polar factor, but the induced spectral shaping of the update.

Third, acceleration mechanisms introduce their own biases and failure modes. Turbo-Muon identifies an AOL-induced “polar bias” that persists as the iteration count grows, even though at practical budgets the decrease in approximation error dominates that bias (Boissin et al., 4 Dec 2025). CANS notes instability in Remez-based optimization for high degrees and reports that degrees beyond $X_{k+1} = X_k(2I - AX_k),$00 did not improve performance empirically (Grishina et al., 12 Jun 2025). HiMuon emphasizes that tiling defines a distinct local operator rather than a convergent approximation to full-matrix Newton–Schulz unless $X_{k+1} = X_k(2I - AX_k),$01, so very small tiles can under-whiten and reduce capacity (Tang et al., 25 Jun 2026).

Finally, Newton–Schulz is not the only route to semi-orthogonal updates. AuON proposes a linear-time alternative that avoids constructing semi-orthogonal matrices altogether, while Hybrid-AuON applies a single Newton–Schulz step before an alternative normalization scheme (Maity, 29 Sep 2025). The existence of such alternatives does not diminish the importance of Newton–Schulz; rather, it clarifies the trade-off space. Newton–Schulz remains the canonical GEMM-only mechanism for approximating polar structure, but modern work increasingly treats its polynomial degree, number of steps, preconditioning, and coupling range as tunable design variables rather than fixed classical formulas.