Chebyshev-Optimized Newton-Schulz (CANS) Method

• CANS is a polynomial-iteration framework that combines the traditional Newton–Schulz method with Chebyshev-optimality to adapt to spectral properties for matrix function approximation.
• It significantly improves convergence rates and reduces computational time in applications like deep learning optimizers and large-scale numerical linear algebra by leveraging efficient matrix multiplications.
• Using techniques such as the Remez algorithm for higher-degree polynomials, CANS offers exponential error reduction while maintaining scalability on modern parallel architectures.

The Chebyshev-Optimized Newton-Schulz (CANS) method is a polynomial-iteration framework for matrix function approximation, specifically tailored for tasks that involve matrix inversion, inverse square roots, or matrix orthogonalization. It synergistically combines the structure of Newton–Schulz (or Hotelling’s) iteration with the minimax property of Chebyshev polynomials, resulting in a spectrally-aware, highly parallelizable, and computationally efficient approach. CANS has demonstrated measurable improvements in convergence and wall-clock performance for large-scale numerical linear algebra, deep learning optimizers (such as Muon), and constrained optimization on matrix manifolds, among other applications (Grishina et al., 12 Jun 2025, Bergamaschi et al., 2020).

1. Classical Newton–Schulz Iteration and Its Limitations

The Newton–Schulz iteration is a canonical matrix-only method for computing the matrix inverse or its inverse square root, appealing for its exclusive use of matrix multiplications, which are highly efficient on modern parallel computing architectures:

$X_{k+1} = \frac{1}{2} X_k (3I - A X_k^2)$

for the inverse square root problem, converging quadratically provided $\|I - A X_0^2\| < 1$. The method is also applicable to orthogonalization, where, given a rectangular $X$, one iterates: $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$ However, the scalar coefficients (e.g., $\frac{3}{2},-\frac{1}{2}$) are fixed as Taylor approximants about the expansion point $t=1$, and do not adapt to the spectral distribution of the matrix at hand. This lack of spectral adaptivity imposes limitations on convergence speed and error uniformity, especially in cases where the eigenvalues of the input matrix are poorly clustered or widely spread (Grishina et al., 12 Jun 2025, Bergamaschi et al., 2020).

2. Chebyshev-Optimality: Deriving Spectrum-Aware Polynomial Updates

CANS introduces Chebyshev-optimality to the construction of matrix polynomial iterations. For the cubic case, the goal is to find the best odd cubic polynomial $p(x) = \alpha_1 x + \alpha_3 x^3$ minimizing the maximum uniform deviation $\epsilon$ from unity on the spectral interval $[a, b]$: $\varepsilon = \max_{x\in[a,b]} |p(x) - 1|.$ Applying Chebyshev’s alternance theorem, there exist extremal points $\|I - A X_0^2\| < 1$0 with alternating error, yielding explicit closed-form optimal coefficients: $\|I - A X_0^2\| < 1$1 Thus, the cubic CANS update is: $\|I - A X_0^2\| < 1$2 This analytic solution is only available in the cubic case; for higher degrees, the construction is performed numerically.

3. Remez Algorithm and Higher-Degree CANS Polynomials

For higher-order odd-degree polynomial approximations ($\|I - A X_0^2\| < 1$3), CANS employs a discrete Remez algorithm on the spectrum interval $\|I - A X_0^2\| < 1$4:

• Initialize with a set of alternation points $\|I - A X_0^2\| < 1$5.
• Solve the linear system: $\|I - A X_0^2\| < 1$6, $\|I - A X_0^2\| < 1$7, to obtain the polynomial coefficients and uniform error.
• Identify the new extremal points, update the alternation set, and iterate until convergence in $\|I - A X_0^2\| < 1$8.

This algorithm computes the unique best uniform odd polynomial $\|I - A X_0^2\| < 1$9 that approximates unity over the specified spectral interval, ensuring robust minimax optimality regardless of the underlying spectrum (Grishina et al., 12 Jun 2025).

4. General Formulation and Iteration Structure

Defining the residual matrix at iteration $X$0 as $X$1, the degree-$X$2 CANS iteration applies the minimax polynomial $X$3:

$X$4

Each term in the polynomial expansion corresponds to a matrix–matrix multiplication, which is efficient on GPU hardware for moderate $X$5, as the algorithm avoids explicit factorizations or square roots, and does not require explicit knowledge of the full spectrum, only its extremal points.

5. Theoretical Error Analysis and Convergence Behavior

Letting $X$6 describe the effective spectral condition number, the minimax polynomial's maximum error takes the form:

$X$7

so that after one iteration: $X$8 and, after $X$9 iterations, the residual decays as $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$0. Thus, increasing the degree $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$1 rapidly reduces the one-step error, yielding an exponential convergence rate improvement compared to the classical fixed-coefficient iteration for fixed spectrum (Grishina et al., 12 Jun 2025, Bergamaschi et al., 2020).

In the context of the Conjugate Gradient (CG) method, using a CANS preconditioner of degree $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$2 transforms the effective condition number as: $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$3 where $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$4, and yields a similar rate reduction in CG residual norm.

6. Computational Complexity and Parallelization

Each degree-$$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$5 CANS step entails:

• One matrix–matrix multiplication for $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$6.
• One multiplication to form $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$7.
• $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$8 multiplications to evaluate powers of $$X_{k+1} = \frac{3}{2} X_k - \frac{1}{2} X_k (X_k^T X_k) X_k.$$9 applied to $\frac{3}{2},-\frac{1}{2}$0.

Thus, the per-iteration complexity is $\frac{3}{2},-\frac{1}{2}$1 large matrix products. On modern GPUs, these operations approach peak compute throughput due to architectural optimization for such kernels; increasing $\frac{3}{2},-\frac{1}{2}$2 trades per-step work for fewer total iterations. In practice, typical choices are $\frac{3}{2},-\frac{1}{2}$3, balancing convergence rate and resource consumption. For parallel sparse linear systems, the method is implemented in a fully matrix-free manner with block-row data distribution across MPI ranks and no global reductions except for the required scalar products in CG (Bergamaschi et al., 2020).

7. Applications in Large-Scale and Machine Learning Contexts

CANS is incorporated in several computational contexts:

A. Deep Learning Orthogonalization (Muon Optimizer):

Muon requires per-step approximate orthogonalization of gradient matrices $\frac{3}{2},-\frac{1}{2}$4, wherein CANS is used with spectrum estimates from $\frac{3}{2},-\frac{1}{2}$5. Empirical studies with NanoGPT document that composite CANS polynomials (of varying degree per iteration) attain singular-value deviation $\frac{3}{2},-\frac{1}{2}$6 and match or slightly improve convergence rates versus alternative polynomials, all with identical matrix-multiplication costs (Grishina et al., 12 Jun 2025).

B. Riemannian Optimization on the Stiefel Manifold:

Retraction of perturbed points $\frac{3}{2},-\frac{1}{2}$7 to the manifold $\frac{3}{2},-\frac{1}{2}$8 is expedited by using a low-degree CANS iteration as an efficient substitute for explicit polar or SVD-based retraction. For instance, in Wide-ResNet experiments on CIFAR-10, CANS retraction matches the accuracy of NLA-based schemes (Cayley or QR) but improves per-epoch training time by $\frac{3}{2},-\frac{1}{2}$9–$t=1$0% on a V100 GPU (Grishina et al., 12 Jun 2025).

C. Preconditioning for Large-Scale Conjugate Gradient:

In systems up to billions of unknowns, CANS preconditioners of degree $t=1$1 reduce CG iteration counts by factors of 10–100 with minor per-iteration overhead. Weighted weak-scaling experiments report time-to-solution speedups between $t=1$2 and $t=1$3 relative to diagonal preconditioning, and up to $t=1$4% reduction in synchronization overhead on over 2,000 MPI ranks (Bergamaschi et al., 2020).

Use Case	Performance Metric	Outcome
Muon/NanoGPT	Singular value deviation	$t=1$5
Stiefel (WRN/CIFAR-10)	Epoch time/speedup	$t=1$6–$t=1$7% faster retraction
CG Preconditioning	CG iteration reduction	Factors of 10–100; up to 2.4$t=1$8 speedup

8. High-Level Algorithmic Realization

For the general degree-$t=1$9 method, CANS computes one update as follows (see procedural details in (Grishina et al., 12 Jun 2025) and (Bergamaschi et al., 2020)):

1. Form $p(x) = \alpha_1 x + \alpha_3 x^3$0.
2. Form $p(x) = \alpha_1 x + \alpha_3 x^3$1.
3. Compute residual $p(x) = \alpha_1 x + \alpha_3 x^3$2.
4. Apply Horner's rule to evaluate $p(x) = \alpha_1 x + \alpha_3 x^3$3 on $p(x) = \alpha_1 x + \alpha_3 x^3$4:

$p(x) = \alpha_1 x + \alpha_3 x^3$7

1. Return $p(x) = \alpha_1 x + \alpha_3 x^3$5.

This composition is compatible with both dense (GPU) and large-scale sparse (distributed, matrix-free) environments. Polynomial coefficients may be precomputed offline using the Remez algorithm for the relevant spectral interval.

9. Relationships, Limitations, and Extensions

CANS generalizes and strictly subsumes both fixed-coefficient Newton–Schulz iterations and Chebyshev polynomial methods, with a proved exact equivalence for parameters $p(x) = \alpha_1 x + \alpha_3 x^3$6. A "de-clustering" modification, implemented by inflating the Chebyshev interval slightly, mitigates extremal eigenvalue condensation and yields measurable acceleration for Krylov methods (Bergamaschi et al., 2020). The method retains its efficacy across orthogonalization, matrix function evaluation, and preconditioning, and requires only approximate knowledge of the spectral extremities of the argument matrix.

Limitations stem from the up-front cost of computing spectral bounds (potentially 5–10 power method or DACG steps) and increased matrix-multiplication cost for higher-degree polynomials. Nevertheless, the algorithm is highly resilient to problem size growth and is particularly advantageous on architectures where matrix multiplications are inexpensive relative to reductions or factorizations.

A plausible implication is that CANS constitutes a unifying computational paradigm for iterative polynomial matrix functions in large-scale numerical optimization and machine learning, especially where matrix-multiplication efficiency is paramount and spectral information is approximately available.