We consider the problem of estimating the spectral norm of a matrix using only matrix-vector products. We propose a new Counterbalance estimator that provides upper bounds on the norm and derive probabilistic guarantees on its underestimation. Compared to standard approaches such as the power method, the proposed estimator produces significantly tighter upper bounds in both synthetic and real-world settings. Our method is especially effective for matrices with fast-decaying spectra, such as those arising in deep learning and inverse problems.
Figures & tables
δ
θcb
g(θcb)
maximizing ρ
0.1
1.40
0.094291
1.2391
0.05
1.76
0.045134
1.1631
0.01
3.00
0.008700
1.0577
0.001
6.45
0.000861
1.0126
Table 1: Numerical maximization of the matrix-agnostic coupled bound g(θ) at the parameter values used in the experiments.
δ
θcb
θv
0.1
1.4
1.72
0.05
1.76
2.17
0.01
3.00
3.70
0.001
6.45
7.98
Table 2: Values of parameter θ for Counterbalance ( 7 ) and Vanilla ( 3 ) estimators at different δ .
Figure 1: Empirical density comparison for different methods
Figure 2: Empirical density comparison with Lanczos-type estimator.
Figure 3: Finite-budget comparison of different estimators as the number of matrix-vector products increases. Solid lines show the mean estimates; shaded areas indicate empirical variance.
δreal
Hilbert
Rank 2
Dominant 0.1
Roget’s Thesaurus
Dominant 0.5
Fréchet derivative
ResNet Jacobian
Counterbalance
0.034
0.031
0.045
0.0001
0.0001
0
0
Vanilla
0.011
0.019
0.016
0.013
0.0
0
0
Dixon
0.019
0.029
0.031
0.022
0.0
0
0
Lanczos-type
0.027
0.030
0.034
0.044
0.04
0
0
Table 3: Values of underestimation probability δreal=P(T(θ,X)≤∥A∥2) measured for different test statistics with 3 matvecs and δ=0.05 .
MAE
Hilbert
Rank 2
Dominant 0.1
Roget’s Thesaurus
Dominant 0.5
Fréchet derivative
ResNet Jacobian
Counterbalance
1.01
1.06
0.97
1.28
1.99
8.65
25.82
Vanilla
2.04
1.98
1.97
3.35
3.77
13.46
36.78
Dixon
1.65
1.60
1.6
2.97
3.26
13.44
38.60
Lanczos-type
1.44
1.25
1.31
2.32
1.78
6.95
19.73
Table 4: MAE for different test statistics δ=0.05 .
Figure 4: Counterbalance and randomized power method comparison
Appendix figures & tables2 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 5: Singular values for matrices from Section 6 .
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.
The Davis-Kahan theorem is a fundamental tool in spectral analysis, providing quantitative control over the distance between the eigenspaces of a symmetric matrix and its perturbation. However, when the matrix dimension is large, computing leading eigenvectors is computationally expensive, limiting the practical use of spectral methods in modern large-scale applications. This paper addresses this problem by proposing an independent Bernoulli sampling scheme and proves that the leading left singular vectors of the subsampled matrix faithfully approximate the target subspace of a low-rank symmetric matrix. Our main result is a subsampled Davis-Kahan bound that gives an explicit error bound depending directly on the sampling probability. The bound reveals the trade-off: the computational cost scales linearly with the sampling probability, while the statistical error scales as the inverse square root of the sampling probability. Our result thus extends the Davis-Kahan theorem to the subsampled setting, enabling scalable spectral analysis of large-scale symmetric matrices.
Huan Qing
School of Economics and Finance, Chongqing University of Technology, Chongqing, 400054, China
Motivated by an application in machine learning optimization, this paper focuses on the challenges of sampling a matrix uniformly from the unit spectral norm ball. It is proven that all singular values of sampled matrices converge to 1 almost surely as the matrix dimensions increase. This result provides the theoretical justification for a proposed simple sampling method applicable for large dimension sizes matching matrices found in modern large language models. Experimental results demonstrate both the convergence of the singular values, as well as the exact and proposed approximate sampling methods.