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.

Sep 27, 2026

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

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.

Kristi Topollai, Anna Choromanska

New York University

Matrix Optimization Newton-Schulz Iteration Muon Optimizer Spectral Methods

Ask PDF

Jul 27, 2026

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

Symmetric non-negative matrix factorization (SymNMF) recovers latent group structure from a dependence matrix, but its dense, quadratic-memory objective has confined prior work to moderate sizes. We present a large-scale GPU study of seven algorithm families (over 30 configurations) on absolute Pearson correlation and tail pairwise dependence matrices from Extreme Value Theory, two proxies for empirical risk-factor estimation on large portfolios. A trace-identity reformulation eliminates all

n \times n

intermediates, so a single GPU reaches

n \approx 10^5

and multi-node distribution scales to

n = 10^6

and beyond. Under a two-phase protocol, eleven methods converge at moderate scale; six remain efficient enough at

n = 10^5

(five AdaGrad-family plus ADMM), and five AdaGrad-family methods still converge at

n = 10^6

: AdaGrad, RMSprop, and three we introduce (Piecewise AdaGrad, Row-Stochastic SVRG, Block-SVRG AdaptGrow). At

n = 10^6

the fastest solver tracks the matrix spectrum: Block-SVRG AdaptGrow wins on the flat, ill-conditioned tail-dependence spectrum, where its lower per-iteration cost decides a long factorization, and full-batch AdaGrad wins on the dominant-low-rank correlation spectrum, where the run is short. We also benchmark spherical K-means as a hard-label baseline: cheaper when angular cluster structure is present, yet provably degenerate once the matrix collapses toward a single common factor, where the soft factorization remains necessary.

Lavinia Ghita, Dhruv Desai, Jake Goldberg +1

NVIDIA

Non-Negative Matrix Factorization Low-Rank Matrix Decomposition Matrix Optimization Adaptive Gradient Methods Stochastic Gradient Descent

Ask PDF

Apr 1, 2026

Inverse-Free Sparse Variational Gaussian Processes

Gaussian processes (GPs) offer appealing properties but are costly to train at scale. Sparse variational GP (SVGP) approximations reduce cost yet still rely on Cholesky decompositions of kernel matrices, ill-suited to low-precision, massively parallel hardware. While one can construct valid variational bounds that rely only on matrix multiplications (matmuls) via an auxiliary matrix parameter, optimising them with off-the-shelf first-order methods is challenging. We make the inverse-free approach practical by proposing a better-conditioned bound and deriving a matmul-only natural-gradient update for the auxiliary parameter, markedly improving stability and convergence. We further provide simple heuristics, such as step-size schedules and stopping criteria, that make the overall optimisation routine fit seamlessly into existing workflows. Across regression and classification benchmarks, we demonstrate that our method 1) serves as a drop-in replacement in SVGP-based models (e.g., deep GPs), 2) recovers similar performance to traditional methods, and 3) can be faster than baselines when well tuned.

Stefano Cortinovis, Laurence Aitchison, Stefanos Eleftheriadis +1

University of Oxford · University of Bristol

Sparse Gaussian Processes Gaussian Processes Variational Inference Gaussian Process Regression Natural Gradient Methods

Ask PDF

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