Orthogonal Polynomial Approximation for Matrix Log Normalization in Global Covariance Pooling

Md Rifat Ur Rahman; Md Raihan Khan; Md Sakib Hossain Shovon; Pietro Liò; Mohammad Ali Moni

cs.CVAug 19, 2026

Orthogonal Polynomial Approximation for Matrix Log Normalization in Global Covariance Pooling

Authors: Md Rifat Ur Rahman, Md Raihan Khan, Md Sakib Hossain Shovon, Pietro Liò, Mohammad Ali Moni

Organizations: NeuronTreeAI · Bangladesh University of Engineering and Technology · Khulna University of Engineering and Technology · KAIST ·

Abstract

Global Covariance Pooling (GCP) improves deep networks by capturing second-order feature statistics, and is especially effective for fine-grained recognition. Because covariance matrices live on the Symmetric Positive Definite (SPD) manifold, a normalization step is required before the Euclidean classifier. The faithful choice is the matrix logarithm (MLN-COV), which maps the SPD manifold to its tangent space; in practice it was abandoned in favour of the matrix square root because its eigendecomposition-based gradient is numerically unstable. We show that this instability is an artifact of computing the logarithm spectrally, not of the logarithm itself. Approximating the logarithm with finite polynomials in the covariance matrix removes the eigendecomposition from both passes: every operation becomes a General Matrix Multiplication (GEMM), the gradient stays bounded on the spectral support of the pre-normalized covariance, and the unstable 1/(lambda_i-lambda_j) term never appears. The key ingredient is a mean-eigenvalue pre-normalization that centres the spectrum near 1, away from the singularity of log, with a scalar post-compensation that returns the singular part of log(A) in closed form. Our recommended normalizer is a degree-8 Chebyshev expansion evaluated by a three-term matrix recurrence, with a matching reverse recurrence for the backward pass; Legendre, Laguerre, Taylor and Pade expansions are studied as controls that isolate the roles of the basis and of the target function. On three fine-grained benchmarks and ImageNet-1k the decomposition-free logarithm is both faster and more accurate than the spectral logarithm and than the square-root approximations it replaces, and at matched basis and degree the log target beats the square-root target, confirming that the gain comes from the faithful Riemannian map rather than from a better polynomial family.

Figures & tables

Panel of Figure 1 — Figure 1: Top: the GCP meta-layer as it exists today; every normalizer it can host targets either the logarithm computed spectrally (MLN-COV [ 11 ] , needing EIG/SVD) or the square root (MPN-COV [ 16 ] , iSQRT-COV [ 15 ] , MPN-MTA/MPN-MPA [ 26 ] ). Bottom: the meta-layer proposed here, where the target is the faithful logarithm and the computation is decomposition-free. MLN-Chebyshev (highlighted) is the recommended setting; MLN-MTA/MLN-MPA are the exact log counterparts of MPN-MTA/MPN-MPA, and MLN-Legendre/MLN-Laguerre vary the orthogonal weight. The contrast between the rows is the contribution: the same box, with the target changed from $A^{1/2}$ to $\log(A)$ and the implementation from EIG/SVD to GEMM.

Polynomial	Coefficients	Recurrence Relation	Expression (first 3 terms)
Chebyshev	$c_{k}$ via quadrature on $[a,b]$ , Eq. ( 4 )	$T_{k+1}=2M(\tilde{A})\,T_{k}-T_{k-1}$	$c_{0}I+c_{1}T_{1}(\tilde{A})+c_{2}T_{2}(\tilde{A})+\cdots$
Legendre	$c_{k}$ via quadrature on $[a,b]$ , Eq. ( 6 )	$(k+1)P_{k+1}=(2k+1)$	$-I+\tfrac{3}{2}(2\tilde{A}-I)-\tfrac{5}{6}(6\tilde{A}^{2}-6\tilde{A}+I)+\cdots$
Legendre	$c_{k}$ via quadrature on $[a,b]$ , Eq. ( 6 )	$M(\tilde{A})P_{k}-kP_{k-1}$
Taylor	$c_{k}=\frac{(-1)^{k+1}}{k}$	N/A	$(\tilde{A}-I)-\tfrac{1}{2}(\tilde{A}-I)^{2}+\tfrac{1}{3}(\tilde{A}-I)^{3}+\cdots$
Padé	Matching Taylor up to $m+n$	N/A	$\tfrac{2(\tilde{A}-I)}{\tilde{A}+I},\;\;\tfrac{6(\tilde{A}-I)+3(\tilde{A}-I)^{2}}{6I+6(\tilde{A}-I)+(\tilde{A}-I)^{2}}$
Laguerre	$c_{0}=-\gamma,\;\;c_{k}=-\tfrac{1}{k}\;\;(k\geq 1)$	$(k+1)L_{k+1}^{(\nu)}=(2k+\nu+1$	$-\gamma I+(\tilde{A}-I)-\tfrac{1}{2}L_{2}(\tilde{A})+\cdots$

Table 1: Polynomial and rational families used to approximate

\log(\tilde{A})

, with the mean-eigenvalue normalization

\tilde{A}=dA/\operatorname{tr}(A)

(Eq. ( 2 )). The scalar

\log(\operatorname{tr}(A)/d)\,I

term is restored via Eq. ( 3 );

M(\tilde{A})=\tfrac{2}{b-a}\tilde{A}-\tfrac{a+b}{b-a}I

is the affine map onto

[-1,1]

\gamma

is the Euler–Mascheroni constant and

\nu

is the Laguerre weight parameter (we use

\nu=0

). The leading terms shown for Legendre use the legacy interval

[a,b]=[0,1]

; the implementation uses

[a,b]=[0.05,3.5]

and quadrature. The Padé column gives the

[1/1]

and

[2/2]

approximants of

\log(1+x)

evaluated at

x=\tilde{A}-I

Figure 2: Why the spectral logarithm is the worst case for gradient stability, and why the polynomial route is not. (a) The diagonal term of the backward pass: $1/\lambda$ for the exact logarithm, unbounded as $\lambda\!\to\!0$ , against $|p_{K}^{\prime}(\lambda)|$ for a degree-8 expansion on $[a,b]$ , bounded by $\approx 12$ . The two coincide over the bulk and separate exactly where the exact operator becomes unstable. Trace normalization places the mean eigenvalue at $1/d$ , deep inside the unstable region; Eq. ( 2 ) places it at $1$ . (b) The off-diagonal term $K_{ij}$ of Eq. ( 8 ), which diverges as the eigenvalue gap closes and which the polynomial formulation never forms. Both panels are closed-form properties of $\log$ and of the degree-8 expansions on $[a,b]=[0.05,3.5]$ ; no training data is involved.

Figure 3: Eigenvalues of $\tilde{A}$ pooled over $3{,}000$ mini-batches, for both backbones. (a) Log-scaled histogram. Eq. ( 2 ) pins the mean at $1$ ; trace normalization would place it at $1/d\approx 0.004$ , marked on the left, deep inside the region where the spectral-logarithm gradient diverges (Fig. 2 a). The shaded band is $[a,b]=[0.05,\,3.5]$ , holding $99.2\%$ of the ResNet-50 eigenvalues and $99.0\%$ of the EfficientNetV2-M ones. The two backbones give almost the same distribution, which is what makes one fixed interval workable across architectures. (b) The same spectrum on a linear axis against the orthogonal weights. The Laguerre weight $e^{-x}$ decays across exactly the region where the spectral mass sits, which is why MLN-Laguerre is weakest in Tables 3 – 5 ; Chebyshev and Legendre are instead fitted on $[a,b]$ , which brackets the whole support.

Normalization	Forward	Backward	Time (ms)
MLN-COV [ 11 ]	EIG	EIG, Eq. ( 8 )	69.5
MPN-COV [ 16 ]	EIG	EIG	67.0
iSQRT-COV [ 15 ]	$3N{=}15$ GEMM	coupled, $\sim\!3N$ GEMM	51.1
MPN-MTA [ 26 ]	$K{=}8$ GEMM	$O(K^{2})$ GEMM + Lyapunov	43.1
MPN-MPA [ 26 ]	$8$ GEMM + solve	$O(K^{2})$ GEMM + Lyapunov	47.9
MLN-MTA (ours)	$K{=}8$ GEMM	$O(K^{2})$ GEMM	44.7

Table 2: Normalization-step cost (forward+backward) on a

256{\times}256

covariance, batch size 32, P100, averaged over many iterations. All polynomial rows use

K{=}8

; Padé is

[4/4]

; iSQRT-COV uses

N{=}5

Newton–Schulz iterations. GEMM counts are per sample.

Normalization	ResNet-50			EfficientNetV2-M
Normalization	CUB	Aircraft	Cars	CUB	Aircraft	Cars
MLN-COV [ 11 ]	$70.4\pm 0.4$	$80.1\pm 0.3$	$86.7\pm 0.2$	$78.5\pm 0.3$	$87.9\pm 0.2$	$91.3\pm 0.2$
MPN-COV [ 16 ]	$71.4\pm 0.3$	$81.2\pm 0.3$	$87.5\pm 0.2$	$79.1\pm 0.3$	$88.6\pm 0.2$	$92.0\pm 0.2$
iSQRT-COV [ 15 ]	$72.8\pm 0.2$	$82.0\pm 0.2$	$88.3\pm 0.2$	$79.6\pm 0.2$	$89.1\pm 0.2$	$92.5\pm 0.2$
MPN-MTA [ 26 ]	$73.1\pm 0.3$	$82.3\pm 0.2$	$88.6\pm 0.2$	$79.9\pm 0.2$	$89.4\pm 0.2$	$92.8\pm 0.2$
MPN-MPA [ 26 ]	$73.2\pm 0.2$	$82.5\pm 0.2$	$88.7\pm 0.2$	$80.0\pm 0.2$	$89.5\pm 0.2$	$92.9\pm 0.2$

Table 3: Top-1 accuracy (%) on three FGVC benchmarks, mean

\pm

std over three seeds, all at

K{=}8

. Spectral MLN-COV is the weakest; making the same logarithm decomposition-free recovers and surpasses it. The italicized rows are matched-basis square-root controls, isolating the role of the Riemannian map from that of the polynomial basis.

Table 7

Normalization	$K{=}6$		$K{=}8$		$K{=}10$
Normalization	Acc	Time	Acc	Time	Acc	Time
MLN-MTA	82.13	34.6	83.61	44.7	84.52	59.6
MLN-MPA	82.28	38.4	83.73	48.3	84.66	63.4
MLN-Legendre	82.47	31.1	83.92	41.1	84.88	56.3
MLN-Laguerre	81.41	30.9	82.91	40.7	83.83	55.5
MLN-Chebyshev	82.56	26.8	84.01	36.4	84.97	50.2

Table 6: Degree ablation

K\!\in\!\{6,8,10\}

on ImageNet-1k with EfficientNetV2-M: Top-1 (%) and FP+BP time (ms). Padé is

[\tfrac{K}{2}/\tfrac{K}{2}]

. The

K{=}8

columns are the rows of Table 5 .

Appendix figures & tables2 assets

Supplementary material from the paper’s appendix.

Appendix

Family	$\alpha_{k}$	$\beta_{k}$	$\gamma_{k}$	$(\tau,\mu)$	$\kappa$
Chebyshev	$2$	$1$	$0$	$\bigl(\tfrac{2}{b-a},-\tfrac{a+b}{b-a}\bigr)$	$0$
Legendre	$\tfrac{2k+1}{k+1}$	$\tfrac{k}{k+1}$	$0$	$\bigl(\tfrac{2}{b-a},-\tfrac{a+b}{b-a}\bigr)$	$0$
Laguerre	$\tfrac{1}{k+1}$	$\tfrac{k+\nu}{k+1}$	$\tfrac{2k+\nu+1}{k+1}$	$(-1,\,0)$	$\nu+1$

Appendix

Table 7: Scalars of the unified recurrence ( 21 ).

Normalization	$d{=}128$			$d{=}256$			$d{=}512$
Normalization	$32$	$64$	$128$	$32$	$64$	$128$	$32$	$64$	$128$
MLN-COV	8.6	17.3	34.6	69.5	138.1	276.5	552.1	1105.8	2198.3
MPN-COV	8.2	16.6	33.2	67.0	132.8	266.1	531.2	1064.5	2115.4
iSQRT-COV	6.3	12.5	25.4	51.1	101.4	202.8	405.6	811.4	1621.0
MPN-MTA	5.2	10.4	21.3	43.1	85.3	171.0	341.1	684.2	1354.5
MPN-MPA	5.8	11.6	23.6	47.9	94.7	190.1	379.8	760.3	1507.2

Appendix

Table 8: Normalization-step cost (forward+backward, ms) across covariance dimension and batch size. All polynomial rows use

K{=}8

; Padé is

[4/4]

; iSQRT-COV uses

N{=}5

Newton–Schulz iterations.

Cost-free Spectral Estimation for Adaptive Newton--Schulz in Matrix Optimizers
Sep 27, 2026Kristi Topollai, Anna ChoromanskaMatrix Optimization Newton-Schulz Iteration
New York University
Low-Rank Dependence Decomposition via Accelerated Symmetric Non-negative Matrix Factorization
Jul 27, 2026Lavinia Ghita, Dhruv Desai, Jake Goldberg +1Non-Negative Matrix Factorization Low-Rank Matrix Decomposition
NVIDIA
Inverse-Free Sparse Variational Gaussian Processes
Apr 1, 2026Stefano Cortinovis, Laurence Aitchison, Stefanos Eleftheriadis +1Sparse Gaussian Processes Gaussian Processes
University of Oxford · University of Bristol

Sep 27, 2026J/K move · Enter open · S save

Cost-free Spectral Estimation for Adaptive Newton--Schulz in Matrix Optimizers

Kristi Topollai, Anna Choromanska

New York University

Matrix Optimization Newton-Schulz Iteration Muon Optimizer Spectral Methods

Open paper AskPDF arXiv

Matrix optimizers such as Muon transform each momentum matrix through an approximate orthogonalization, typically implemented by a small number of Newton-Schulz matrix multiplications. The quality and cost of this approximation depend strongly on the singular-value spectrum of its input, yet existing implementations use the same fixed polynomial routine for every layer and throughout training. We show that this uniform treatment is unnecessary: the computations in the Newton-Schulz method already reveal enough information to make the method adaptive. The Gram matrices formed inside Newton-Schulz iterations yield spectral moments through inexpensive scalar reductions, requiring no additional matrix multiplications. From these moments, we recover an estimate of the empirical singular-value distribution and use it to select a polynomial routine specialized to the current matrix. This turns Newton--Schulz orthogonalization into a spectrum-adaptive procedure that responds to differences across both layers and training time. On saved momentum matrices, spectral estimation substantially reduces orthogonalization error at a fixed iteration budget or reaches the same accuracy with fewer iterations, and in GPT pretraining up to 1B parameters it lowers the validation loss of two matrix optimizers. Our results suggest that matrix-function operations inside optimizers need not be designed for a conservative worst-case spectrum: they can cheaply measure the spectrum they are already processing and specialize computations accordingly.

Abstract

Explore similar work

Cost-free Spectral Estimation for Adaptive Newton--Schulz in Matrix Optimizers

Low-Rank Dependence Decomposition via Accelerated Symmetric Non-negative Matrix Factorization

Inverse-Free Sparse Variational Gaussian Processes