Newton-Schulz Iteration

Purpose: Compute a matrix inverse (or pseudo-inverse) iteratively using only matrix multiplication — no division, no pivoting, and with quadratic convergence:

$X_{k+1}=X_k(2I-A{X}_k).$

The Scalar Analogue

For a number, Newton’s method applied to $f(x)=1/x-a$ gives

$x_{k+1}=x_k(2-a{x}_k),$

which converges quadratically to $1/a$ from any start with $|1-a{x}_0|<1$. The matrix version is the same formula with $1\to I$.

Convergence

Define the residual $E_k=I-A{X}_k$. Then

$E_{k+1}=I-A{X}_k(2I-A{X}_k)=(I-A{X}_k{)}^2=E_k^2.$

So $\Vert E_{k+1}\Vert \le \Vert E_k{\Vert }^2$ — the number of correct digits doubles each iteration. Convergence requires $\Vert I-A{X}_0\Vert <1$; the standard safe start is

$X_0=\frac{A^T}{\Vert A{\Vert }_1\Vert A{\Vert }_{\infty }}.$

Complexity

• Per iteration: two matrix multiplications, $O(n^3)$ (or $O(n^{\omega })$)
• Iterations: $O(\log \log (1/\epsilon ))$ for accuracy $\epsilon$, plus $O(\log \kappa )$ to get into the convergence basin
• Total: worse than $O(n^3)$ Gauss-Jordan in the sequential sense

Why use it at all

Not for a single dense inverse on one CPU. It wins where matrix multiplication is much cheaper than elimination:

• GPUs and distributed systems. Matrix multiplication is embarrassingly parallel; Gaussian elimination is inherently sequential in its pivot order. Newton-Schulz turns inversion into a handful of GEMM calls.
• No division. Useful in settings where division is expensive or undefined — fixed-point hardware, homomorphic encryption, certain algebraic rings.
• Refinement. Given an approximate inverse from any source, one Newton-Schulz step squares the error. This is the standard way to “polish” a low-precision inverse into a high-precision one.

The same iteration for power series

The identical Newton scheme is how formal power series operations are computed in $O(n\log n)$:

Operation	Newton iteration
Series inverse	$B_{k+1}=B_k(2-A{B}_k)$
Square root	$B_{k+1}=\frac{1}{2}(B_k+A/B_k)$
Exponential	$B_{k+1}=B_k(1-\log B_k+A)$
Reciprocal of an integer	$x_{k+1}=x_k(2-a{x}_k)$ — used inside big-integer division

Each doubles the number of correct coefficients, so the total cost is dominated by the last step: $T(n)=T(n/2)+O(n\log n)=O(n\log n)$. Recognising Newton’s method in all these places is genuinely useful — it turns “how do I invert a power series?” into a solved problem.

Variants / Use Cases

• Moore-Penrose pseudo-inverse — the same iteration converges to $A^{+}$ for rectangular $A$
• Matrix sign function and polar decomposition — closely related iterations
• Formal Power Series — where this iteration matters most in competitive programming
• Hensel lifting — the $p$-adic version of the same doubling scheme
• Newton’s Method — the scalar root-finding original