---
title: Newton–Schulz Orthogonalization
url: https://www.emergentmind.com/topics/newton-schulz-orthogonalization
type: topic
---

# Newton–Schulz Orthogonalization

Newton–Schulz Orthogonalization refers to a class of iterative, matrix-polynomial algorithms designed to efficiently orthogonalize matrices, particularly in scenarios where explicit QR factorization or singular value decomposition (SVD) is computationally prohibitive. By leveraging only matrix–matrix multiplications, the Newton–Schulz (NS) scheme is highly suitable for large-scale problems, GPU acceleration, and applications requiring repeated or controlled approximate orthogonalization. Key domains of application include eigensolvers, modern stochastic and geometry-aware optimizers (such as Muon), and Riemannian optimization on matrix manifolds.

## 1. Mathematical Formulation and Derivation

The Newton–Schulz orthogonalization process is rooted in the polar decomposition of a nonsingular matrix $A \in \mathbb{R}^{n \times n}$, which admits the unique factorization
$$
A = H P
$$
where $H$ is orthogonal ($H^T H = I$), and $P = (A^T A)^{1/2}$ is symmetric positive definite. In terms of the SVD ($A = U \Sigma V^T$), $H = U V^T$ and $P = V \Sigma V^T$.

To extract the orthogonal factor $H$ without explicit SVD, NS leverages Newton’s method on the matrix equation $F(X) = X^T X - I$:
$$
X_{k+1} = \frac12 (X_k + X_k^{-T})
$$
This iteration converges quadratically to $H$ for nonsingular $A$. By avoiding explicit inversion, a Schulz-type approximation delivers the update:
$$
X_{k+1} = X_k \cdot \frac12 (3I - X_k^T X_k)
$$
or equivalently,
$$
X_{k+1} = X_k (3I - X_k^T X_k)/2
$$
The quadratic convergence is preserved when the spectrum of $A$ (or $X_k^T X_k$) is inside $(0, \sqrt{3})$, and convergence is rapid once the singular values of $A$ are well-scaled [2509.00495].

## 2. Algorithmic Implementation and GPU Efficiency

Newton–Schulz orthogonalization is particularly attractive for large matrices and GPU contexts due to its reliance solely on matrix–matrix products, with no data-dependent branching or explicit orthonormalization. The practical mixed-precision eigensolver in [2509.00495] demonstrates MATLAB pseudocode making direct use of the above iteration, requiring only 2–3 matrix multiplications for double precision accuracy. Similarly, in Muon-style optimizer kernels, each NS step is implemented as a series of batched fused GEMMs, enabling hardware-level acceleration [2512.04632][2506.10935][2509.00495].

Polynomial variants, such as fixed-degree (cubic/quintic) and Chebyshev-optimized schemes, can be constructed to further tailor accuracy/cost tradeoffs, with higher-degree polynomials produced using a Remez algorithm for best approximation on the target spectral interval [2506.10935]. This enhances convergence per matrix-multiply, at the cost of increased per-iteration operations.

## 3. Theoretical Convergence, Error Bounds, and Conditioning

Each iteration contracts the orthogonality residual $R_k = X_k^T X_k - I$ quadratically:
$$
R_{k+1} = \frac14 (3I - R_k) R_k^2
$$
Hence, for any subordinate matrix norm,
$$
\|R_{k+1}\| \leq \frac{3+\|R_k\|}{4} \|R_k\|^2 = O(\|R_k\|^2)
$$
In practice, with appropriate initial scaling (e.g., spectral or Frobenius norm), 2–5 steps suffice to reduce errors to machine epsilon [2509.00495][2512.04632]. For applications involving ill-conditioned matrices (e.g., low-rank, anisotropic gradient moments), convergence constants depend on the condition number $\kappa$, and residual errors decay geometrically as $(1-1/\kappa)^{2^i}$ after $i$ NS steps [2505.24749]. This dependence motivates either specialization to well-conditioned subspaces (as in SUMO) or improvement via polynomial tuning [2506.10935].

## 4. Newton–Schulz Orthogonalization in Optimization Algorithms

Newton–Schulz orthogonalization has been extensively adopted in modern optimizers designed for highly anisotropic or constrained landscapes. The Muon optimizer, in particular, performs spectral norm–inducing step normalization by applying $q$ NS steps of degree-$\kappa$ polynomial to the momentum matrix, constructing
$$
O_t = X_{t,q}, \quad X_{t,j+1} = p_\kappa(X_{t,j} X_{t,j}^T) X_{t,j}
$$
where $p_\kappa(\lambda)$ is a carefully constructed Taylor or Chebyshev polynomial [2601.19156][2510.19933]. For typical $q = 2$–$5$, the orthogonality error decays doubly-exponentially, so that the convergence rate gap to the ideal SVD-polar update becomes negligible after few iterations [2601.19156]. SUMO, in contrast, uses exact SVD in a low-rank subspace to circumvent the condition-number dependency encountered by standard NS [2505.24749].

A critical practical insight is the fundamental coupling between NS approximation error and optimizer hyperparameters: smaller NS iteration count (i.e., looser orthogonalization) demands smaller learning rate and higher momentum to maintain optimization stability [2510.19933].

## 5. Polynomial Acceleration and Chebyshev Optimized Variants

The generic NS update can be recast in a polynomial framework: at each step, the SVD decomposition is used to show that
$$
X_{k+1} = U p_d(S_k) V^T,\quad p_d(x) = \sum_{j=0}^{d-1} \alpha_{2j+1} x^{2j+1}
$$
Chebyshev-optimized Newton–Schulz (CANS) replaces the fixed-degree polynomial of classical NS by constructing best uniform odd polynomial approximations (using the alternation theorem and Remez algorithm) on the estimated spectral interval $[a, b]$ of the current iterate [2506.10935]. In practice, this polynomial tailoring yields a reduction in iteration count for a given level of residual, directly translating to computational savings on accelerators. However, construction becomes numerically unstable for polynomial degrees beyond 5, so practical implementations favor low-degree CANS.

## 6. Preconditioning and Algorithmic Enhancements

Preconditioning is key to further improving Newton–Schulz orthogonalization performance, especially for high aspect-ratio or badly conditioned matrices. The "Almost-Orthogonal Layer" (AOL) diagonal preconditioner rescales input columns to shrink the spectral spread before NS iteration, reducing necessary step count and total number of expensive matrix multiplications [2512.04632]. Empirically, this enables the removal of one iteration in the Muon pipeline—translating to a ∼20% reduction in orthogonalization kernel cost and up to 10% end-to-end optimizer speedup on modern hardware, with no loss in output accuracy.

## 7. Applications, Limitations, and Comparative Analysis

**Applications** of Newton–Schulz orthogonalization encompass:
- Mixed-precision eigensolvers: Achieving double-precision spectral decompositions by orthogonalizing single-precision approximate eigenvectors [2509.00495].
- Large-scale stochastic and geometry-aware optimization: Core to Muon, Turbo-Muon, and similar step-normalizing optimizers [2601.19156][2512.04632][2510.19933].
- Efficient retractions in Riemannian optimization on the Stiefel manifold, providing controlled accuracy at fixed computational cost [2506.10935].

**Comparative performance analyses** show that Newton–Schulz orthogonalization outperforms classical Householder QR in both speed and accuracy for large $n$, owing to optimized BLAS matmul throughput and quadratic error decay. While exact SVD or QR orthogonalization is more stable in highly ill-conditioned or subspace-adaptive scenarios, matrix-polynomial/NS variants often dominate in overall wall-clock efficiency for moderately conditioned, full-rank contexts [2509.00495][2506.10935][2505.24749].

**Limitations** include slower convergence in the presence of extreme condition numbers (e.g., during LLM training with rapidly evolving gradients) unless augmented by preconditioning or restricted to adaptively chosen subspaces (as in SUMO).

**Summary Table: Algorithmic Tradeoffs (processes and outcomes present verbatim in the data)**

| Method          | Matmuls per Step | Orthogonalization Error    |
|-----------------|------------------|---------------------------|
| Householder QR  | ≈(10/3) $n^3$    | $\lesssim 10^{-13}n$      |
| Classical NS    | 2–3 $n^3$        | $\lesssim 10^{-14}$       |
| Chebyshev NS    | $d$ $n^3$        | $O(\varepsilon^k)$        |
| SUMO (SVD)      | O($nr^2$)        | Exact in subspace         |

This table summarizes the computational and accuracy properties for the main orthogonalization approaches, as numerically reported in [2509.00495][2506.10935][2505.24749].

## References

- "A Mixed Precision Eigensolver Based on the Jacobi Algorithm" [2509.00495]
- "Beyond the Ideal: Analyzing the Inexact Muon Update" [2510.19933]
- "Accelerating Newton-Schulz Iteration for Orthogonalization via Chebyshev-type Polynomials" [2506.10935]
- "Convergence of Muon with Newton-Schulz" [2601.19156]
- "SUMO: Subspace-Aware Moment-Orthogonalization for Accelerating Memory-Efficient LLM Training" [2505.24749]
- "Turbo-Muon: Accelerating Orthogonality-Based Optimization with Pre-Conditioning" [2512.04632]

Source: https://www.emergentmind.com/topics/newton-schulz-orthogonalization