Newton–Schulz Iterations: Matrix Inversion & Beyond

• Newton–Schulz iterations are a family of fixed-point schemes that compute matrix inverses and polar factors with quadratic or higher-order convergence.
• High-order extensions employ optimized polynomial factorizations to minimize matrix multiplications while achieving efficiency indices around 1.41–1.49.
• Integrating these iterations with Richardson’s method enables robust preconditioning and fast orthogonalization in applications like failure detection and machine learning.

The Newton–Schulz iteration is a family of high-order fixed-point schemes for constructing matrix inverses, orthogonal polar factors, and polynomial approximations of linear operators. Originating in classical work on iterative inversion and rapidly expanding applications in computational linear algebra, optimization, and machine learning, the family encompasses both the original second-order Newton–Schulz iteration and its numerous high-order and composite generalizations. Central to its efficacy is its amenability to efficient implementation using only matrix-matrix multiplications, making it especially relevant in contexts such as GPU-accelerated computation, preconditioning of ill-conditioned or rank-deficient systems, and fast orthogonalization in large-scale optimization algorithms.

1. General Formulation and Basic Properties

The classical (second-order) Newton–Schulz iteration for inverting a nonsingular matrix $A$ is given by

$$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$

where $X_0$ is an initial preconditioner. Defining the residual $F_k = I - X_k A$, the iteration can be recast as

$X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$

which reveals the quadratic convergence: the norm of the residual is squared at each step. This iteration is a special case of the family of hyperpower or Durand–Schulz methods, whose $n$-th order generalization employs the truncated Neumann series: $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$ and $F_{k+1} = F_k^n$, ensuring local convergence of order $n$ (i.e., $\|F_{k+1}\| \leq C \|F_k\|^n$) (Stotsky, 2022, Stotsky, 2020).

The convergence of these schemes is ensured under the spectral condition $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$0 (where $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$1 denotes spectral radius), typically guaranteed by choosing $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$2 for appropriate $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$3. For symmetric positive definite (SPD) $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$4, taking $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$5 with eigenvalues in the reciprocal range of $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$6 suffices.

2. High-Order Extensions and Unified Polynomial Factorization

While a naïve implementation of the $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$7-th order Newton–Schulz scheme would entail computing each power $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$8 for $$X_{k+1} = 2X_k - X_k A X_k, \qquad X_0 \approx A^{-1}$$9, sophisticated factorization of the underlying residual polynomial significantly reduces computational cost. Given

$X_0$0

one can decompose $X_0$1 into products of lower-degree factors, enabling a reduction in the number of matrix multiplications per iteration.

• If $X_0$2, the factorization

$X_0$3

requires only $X_0$4 multiplications (rather than the naïve $X_0$5).

• The efficiency index (EI) $X_0$6 quantifies the per-multiply convergence acceleration.

Typical factorization examples (with 6 products per cycle): | $X_0$7 | Polynomial Factorization $X_0$8 | $X_0$9 | |-----|---------------------------------------------------------|----------| | 8 | $F_k = I - X_k A$0 | $F_k = I - X_k A$1 | | 9 | $F_k = I - X_k A$2 | $F_k = I - X_k A$3 | | 10 | $F_k = I - X_k A$4 | $F_k = I - X_k A$5 | | 11* | Bespoke (prime) factorization | $F_k = I - X_k A$6 |

Optimal choice of $F_k = I - X_k A$7 and factorization yields a balance between per-step cost and asymptotic convergence (Stotsky, 2022, Stotsky, 2020).

3. Composite and Parallelized Newton–Schulz Schemes

Composite expansions can further enhance computational throughput, especially on parallel hardware. By writing the high-order polynomial as a product of several low-order terms, one schedules independent blocks across compute units with different capabilities. Given exponents $F_k = I - X_k A$8, define

$F_k = I - X_k A$9

and perform the update $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$0, ensuring rapid decay of the residual $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$1 controlled by the composite contraction $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$2. Recursive Horner-type schemes and nested factorizations enable scaling to high orders while maintaining minimal multiply count. For instance, a 45th-order Newton–Schulz iteration can be realized in only 10 matrix-matrix multiplies (Stotsky, 2020).

4. Integration with Richardson Iteration

Richardson’s method solves $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$3 via the fixed-point scheme

$X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$4

with slow practical convergence. By embedding a Newton–Schulz approximate inverse $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$5,

$X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$6

or, equivalently,

$X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$7

and updating $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$8 via a high-order Newton–Schulz step at each outer Richardson iteration, one achieves both high-order convergence of the inverse and accelerated reduction of the solution residual. Skeleton pseudocode for the combined scheme:

1. Initialize $X_{k+1} = (I + F_k) X_k, \qquad F_{k+1} = F_k^2$9 (preconditioner), $n$0 (e.g., zero).
2. For $n$1 do:
   • Compute residual $n$2.
   • Newton–Schulz update: $n$3 (with efficient factorization).
   • Set $n$4.
   • Richardson update: $n$5.
3. Repeat until convergence.

This nested structure enables dramatic acceleration, especially when robust preconditioning is required (e.g., ill-conditioned or rank-deficient $n$6) (Stotsky, 2022, Stotsky, 2020).

5. Practical Guidelines and Computational Trade-Offs

Key considerations for real-world use:

• For SPD problems of moderate size ($n$7), 8th–9th order factorizations ($n$8–$n$9) provide a favorable balance of speed and stability.
• In very large or sparse systems, repeated low-order ($$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$0) Newton–Schulz steps per Richardson iteration can increase robustness and reduce implementation complexity.
• For matrices with condition number $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$1, a single 4th–6th order NS update per iteration generally suffices; more extreme ill-conditioning may warrant higher $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$2 or tighter inner-outer iteration nesting.
• The spectral radius condition $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$3 is ensured in practice by setting $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$4 with $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$5. A typical choice is $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$6.
• Parallel compute resources can be exploited by distributing factorizations across units according to their computational capabilities (Stotsky, 2022, Stotsky, 2020).

6. Applications in Robust Failure Detection and Machine Learning

Failure Detection in Electrical Networks

A central application is in robust parameter estimation for failure detection in electrical networks, where streaming least-squares problems are afflicted by nearly rank-deficient design matrices (e.g., due to missing or faulty data channels): $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$7 Regularizing with small $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$8: $$X_{k+1} = \left(I + F_k + F_k^2 + \cdots + F_k^{n-1}\right) X_k$$9, the systems frequently become severely ill-conditioned. Adaptive preconditioning with

$F_{k+1} = F_k^n$0

ensures $F_{k+1} = F_k^n$1. The combined Newton–Schulz–Richardson iteration yields stable solutions even as $F_{k+1} = F_k^n$2 ($F_{k+1} = F_k^n$3), enabling accurate tracking of voltage sags/swells where direct (e.g., LU-based) solvers fail. On real DOE/EPRI data, this approach achieves real-time performance with fewer than 10 Richardson steps per window shift and recovers first-harmonic amplitudes with high fidelity (Stotsky, 2022).

Orthogonalization in Machine Learning Optimizers

In large-scale optimization (e.g., the Muon optimizer), Newton–Schulz iterations provide fast, inexact polar decomposition updates as cheap alternatives to SVD-based projection. The key recurrence

$F_{k+1} = F_k^n$4

is run for a small, fixed number of steps (typically 5), balancing approximation error and computational budget (Shulgin et al., 22 Oct 2025).

Theory predicts that lowering the number of Newton–Schulz steps (incurring greater inexactness $F_{k+1} = F_k^n$5) degrades convergence rate by explicit factors (e.g., a deterministic rate inflated by $F_{k+1} = F_k^n$6). Optimal learning rate and momentum must therefore be co-tuned with the number of Newton–Schulz steps: fewer iterations require a smaller step size and increased momentum, as confirmed in empirical training of deep models such as NanoGPT (Shulgin et al., 22 Oct 2025).

Accelerated variants, such as Chebyshev-optimized Newton–Schulz (CANS), further reduce matmul cost and improve contraction rates by minimizing the worst-case approximation error over the relevant singular value interval, and their use in Muon or Stiefel-manifold optimization achieves faster convergence and lower wall-clock time than classical Newton–Schulz or QR/Cayley retractions (Grishina et al., 12 Jun 2025).

7. Conclusions and Outlook

The Newton–Schulz family, bolstered by unified polynomial factorizations and composite-block constructions, constitutes a highly flexible and tunable class of iterative solvers for inverse, orthogonal, and system-solution computations. Their seamless integration into outer iterative schemes (notably Richardson iteration) imparts both finite-precision safety and exceptional adaptability to scale, conditioning, and computational infrastructure. In both numerical linear algebra and modern machine learning, these schemes have demonstrated robust real-world utility—outperforming direct solvers in speed, robustness, and flexibility across challenging application domains (Stotsky, 2022, Stotsky, 2020, Shulgin et al., 22 Oct 2025, Grishina et al., 12 Jun 2025). Future research directions involve optimal composite design for particular hardware or problem structure, improved contraction analysis for blockwise or layerwise errors, and principled adaptive selection of the order and composition at runtime in large-scale distributed systems.