% Source: https://kexue.fm/archives/11736
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage{booktabs,longtable,array}
\usepackage{xcolor}
\usepackage{listings}

\geometry{a4paper, margin=1in}
\definecolor{lightgray}{gray}{0.95}
\definecolor{skyblue}{rgb}{0.0,0.7,1.0}
\lstset{
    backgroundcolor=\color{lightgray},
    basicstyle=\ttfamily\footnotesize,
    numbers=left,
    numberstyle=\tiny,
    frame=single,
    extendedchars=true,
    literate={*}{}{0\columnwidth}
}

\title{How to Estimate Matrix Spectral Norms More Scientifically?}
\author{Su Jianlin}
\date{2026-05-04}

\begin{document}
\maketitle

\section{Spectral Norm} \label{spectral-norm}
Spectral norm is a core concept in matrix analysis and plays an important role in deep learning---from Lipschitz constraints required in the WGAN era to training stability in modern LLMs, and to the emerging Muon optimizer, all are closely related to the spectral norm of matrix parameters. Therefore, how to efficiently and accurately estimate spectral norms has become increasingly important and deserves in-depth study.

It is well known that power iteration is the standard method for estimating spectral norms, but there is still considerable room for improvement. This article briefly reviews several ideas for estimating spectral norms, including accelerating the convergence of power iteration and estimating strict upper bounds of spectral norms.

The spectral norm is defined as
\begin{equation}
\Vert\boldsymbol{W}\Vert_2 = \max_{\Vert\boldsymbol{x}\Vert_2=1} \Vert \boldsymbol{W}\boldsymbol{x}\Vert_2
\end{equation}
where $\boldsymbol{W}\in\mathbb{R}^{n\times m}, \boldsymbol{x}\in\mathbb{R}^m$, and without loss of generality, assume $n\geq m$. By definition, the spectral norm represents the ``expansion rate'' of a linear layer---after an input vector passes through the linear layer, its norm expands by at most a factor of $\Vert\boldsymbol{W}\Vert_2$. The spectral norm is also equal to the largest singular value of the matrix (see the proof in \href{https://kexue.fm/archives/10407}{``Low-Rank Approximation Path (II): SVD''}). We won't elaborate on these basic concepts further.

My first encounter with spectral norms was during the rise of GANs, when \href{https://kexue.fm/archives/6280}{WGAN} emerged, emphasizing the necessity for discriminators to satisfy Lipschitz constraints. The Lipschitz constant of a linear layer $\boldsymbol{W}\boldsymbol{x}$ is precisely the spectral norm of matrix $\boldsymbol{W}$, leading to techniques such as spectral normalization and spectral regularization. Interested readers may refer to \href{https://kexue.fm/archives/6051}{``Lipschitz Constraints in Deep Learning: Generalization and Generative Models''}.

In recent years, the increasingly popular \href{https://kexue.fm/tag/muon/}{Muon} optimizer is also deeply connected with spectral norms, often viewed as steepest descent under the spectral norm. Meanwhile, in \href{https://kexue.fm/archives/11729}{``Beyond MuP: 4. Maintaining Parameter Stability''}, we proposed constraining the spectral norm of parameters to ensure training stability. These developments demand accurate computation of spectral norms.

\section{Power Iteration} \label{power-iteration}
Computing the spectral norm via SVD is the most direct approach, but clearly too expensive. We usually use power iteration:
\begin{equation}
\boldsymbol{v}^{(t)} = \frac{\boldsymbol{W}^{\top}\boldsymbol{W}\boldsymbol{v}^{(t-1)}}{\Vert\boldsymbol{W}^{\top}\boldsymbol{W}\boldsymbol{v}^{(t-1)}\Vert_2}
\end{equation}
which converges to the right principal singular vector $\boldsymbol{v}_1$ at rate $(\sigma_2/\sigma_1)^{2t}$. After $T$ iterations, we have
\begin{equation}
\sigma_1 \approx \Vert\boldsymbol{W}\boldsymbol{v}^{(T)}\Vert_2
\end{equation}
The proof is available in \href{https://kexue.fm/archives/10648}{``From Spectral Norm Gradient to Modern Weight Decay Considerations''}. In practice, if we only care about the spectral norm and not the singular vectors, the convergence speed is often more optimistic than $(\sigma_2/\sigma_1)^{2t}$. Additionally, the result of power iteration is easily affected by the initial value $\boldsymbol{v}^{(0)}$, causing fluctuations. We can consider computing $k$ paths in parallel and then taking the maximum.

If we compute $k$ paths simultaneously and replace the L2 normalization at each step with orthogonalization of these $k$ vectors, i.e.,
\begin{equation}
\boldsymbol{V}^{(t)} = \mathrm{QR}(\boldsymbol{W}^{\top}\boldsymbol{W}\boldsymbol{V}^{(t-1)})
\end{equation}
then the result converges to the first $k$ right singular vectors, enabling computation of the first $k$ singular values. The principle can be found in \href{https://kexue.fm/archives/11710}{``Muon Implementation Based on Streaming Power Iteration: 4. Principle''}. However, we won't expand on this extension here and will instead focus on estimating the spectral norm---i.e., the largest singular value.

\section{Taking Gradients} \label{taking-gradients}
First, here is a reference implementation in Jax for power iteration:
\begin{lstlisting}[language=Python]
import jax, jax.lax as lax
import jax.numpy as jnp

@jax.jit(static_argnums=(1,))
def l2_normalize(x, axis=-2):
    return x / jnp.linalg.vector_norm(x, axis=axis, keepdims=True)

@jax.jit(static_argnums=(1, 2, 3))
def spec_norm_v1(w, T=10, k=1, key=42):
    v_shape = w.shape[:-2] + w.shape[-1:] + (k,)
    v_init = jax.random.normal(jax.random.PRNGKey(key), v_shape)
    v_step = lambda i, v: l2_normalize(w.mT @ (w @ v))
    v = lax.fori_loop(0, T, v_step, v_init)
    return jnp.linalg.vector_norm(w @ v, axis=-2).max(axis=-1)
\end{lstlisting}
It is clear that the complexity of power iteration is $\mathcal{O}(mnTk)$, which is already quite ideal. If we only need monitoring, this implementation is sufficient. However, for purposes such as spectral normalization or spectral regularization, we need gradients of the spectral norm. With the above implementation, backpropagation must go through the power iteration loop, resulting in high computational complexity.

In fact, in \href{https://kexue.fm/archives/10648}{``From Spectral Norm Gradient to Modern Weight Decay Considerations''}, we derived the gradient of the spectral norm $\nabla_{\boldsymbol{W}}\sigma_1 = \boldsymbol{u}_1 \boldsymbol{v}_1^{\top}$, so we can directly define this gradient to reduce the complexity of backpropagation. Alternatively, we can use $\sigma_1 = \boldsymbol{u}_1^{\top}\boldsymbol{W}\boldsymbol{v}_1$ and output during the forward pass:
\begin{equation}
\sigma_1 = \color{skyblue}{\mathrm{sg}}[\boldsymbol{u}_1^{\top}] \boldsymbol{W} \color{skyblue}{\mathrm{sg}}[\boldsymbol{v}_1]
\end{equation}
where $\color{skyblue}{\mathrm{sg}}[\cdot]$ is the stop gradient operator. Thus, during automatic differentiation, only $\boldsymbol{W}$ is differentiated, yielding $\boldsymbol{u}_1 \boldsymbol{v}_1^{\top}$, avoiding backpropagation into the internal power iterations of $\boldsymbol{u}_1,\boldsymbol{v}_1$.

\section{No Waste} \label{no-waste}
Through $T$ steps of power iteration, we obtain the sequence $\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)},\cdots,\boldsymbol{v}^{(T)}$, but ultimately we only use $\boldsymbol{v}^{(T)}$ to estimate $\sigma_1$, which seems quite ``wasteful''. Therefore, some have attempted to utilize all of them to achieve better estimates.

Unless unusual, $\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)},\cdots,\boldsymbol{v}^{(T)}$ are different but increasingly close to $\boldsymbol{v}_1$. We can think of them as ``$\boldsymbol{v}_1$ plus some noise'', and combining them may yield more accurate results. Specifically, we consider constructing a better approximation for $\boldsymbol{v}_1$ via linear combinations of $\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)},\cdots,\boldsymbol{v}^{(T)}$. The subspace spanned by these vectors is called the ``Krylov subspace.''

For simplicity, perform one orthogonalization to obtain an equivalent orthonormal basis $\boldsymbol{Q} = [\boldsymbol{q}_1,\boldsymbol{q}_2,\cdots,\boldsymbol{q}_T]\in\mathbb{R}^{m\times T}$; then find coefficients $\boldsymbol{x}=[x_1,x_2,\cdots,x_T]^{\top}\in\mathbb{R}^T$ such that $\sum_{i=1}^T x_i \boldsymbol{q}_i = \boldsymbol{Q}\boldsymbol{x}$ satisfies
\begin{equation}
\boldsymbol{W}^{\top}\boldsymbol{W}\boldsymbol{Q}\boldsymbol{x} \approx \sigma_1^2\boldsymbol{Q}\boldsymbol{x}
\end{equation}
Multiplying both sides by $\boldsymbol{Q}^{\top}$ gives $\boldsymbol{Q}^{\top}\boldsymbol{W}^{\top}\boldsymbol{W}\boldsymbol{Q}\boldsymbol{x} \approx \sigma_1^2\boldsymbol{x}$, indicating that $\sigma_1^2, \boldsymbol{x}$ are the eigenvalue and eigenvector of $(\boldsymbol{W}\boldsymbol{Q})^{\top}\boldsymbol{W}\boldsymbol{Q}$. This means we only need to perform eigenvalue decomposition, take the largest eigenvalue, and take its square root for improved estimation. Since this is only a $T\times T$ matrix, when $T$ is small, eigenvalue decomposition is not costly, so calling existing functions suffices.

This approach is essentially a simplified version of the ``Lanczos algorithm'' used in modern SVD solvers, with actual SVD employing more refined processing. Additionally, this acceleration idea shares similarities with randomized SVD, which typically uses Gaussian random matrices for projection, while here we use orthonormal bases in the Krylov subspace for projection.

\section{Acceleration Trick} \label{acceleration-trick}
The reference code for power iteration incorporating Krylov subspace acceleration is:
\begin{lstlisting}[language=Python]
@jax.jit(static_argnums=(1, 2))
def spec_norm_v2(w, T=10, key=42):
    v_shape = w.shape[:-2] + w.shape[-1:] + (1,)
    v_init = jax.random.normal(jax.random.PRNGKey(key), v_shape)
    v_step = lambda v, _: (l2_normalize(w.mT @ (w @ v)),) * 2
    v = lax.scan(v_step, v_init, length=T)[1].swapaxes(0, -1)[0]
    v = jnp.linalg.qr(v)[0]
    return jnp.linalg.eigvalsh((u := w @ v).mT @ u)[..., -1]**0.5
\end{lstlisting}
Practical testing shows that subspace acceleration is highly effective: when $T=5$, its accuracy exceeds that of the original power iteration at $T=10$, and it is faster. At $T=10$, about $2/3$ of the function's runtime is due to power iteration, $1/3$ due to QR decomposition, with eigenvalue decomposition accounting for a very small portion.

Readers familiar with the ``streaming power iteration'' series might consider using Cholesky QR introduced there to accelerate QR decomposition, but unfortunately, as power iteration results grow increasingly collinear, the condition number of the matrix to be QR-decomposed deteriorates continuously, which is exactly where Cholesky QR ``fails'', with a very high failure rate, making it basically unusable.

A more effective acceleration technique is to use only the last few steps of the power iteration, e.g., $\boldsymbol{v}^{(T-2)},\boldsymbol{v}^{(T-1)},\boldsymbol{v}^{(T)}$, to construct the Krylov subspace acceleration. This already captures most of the benefits while reducing computational costs of QR and eigenvalue decomposition, thus achieving speedup.

\section{Estimating Upper Bound} \label{estimating-upper-bound}
Power iteration and its accelerated variant can strictly speaking only obtain a lower bound of the spectral norm, though this is often sufficient. However, in some scenarios, we need an upper bound---for example, if we want to enforce spectral normalization ensuring that a matrix's norm is strictly less than 1, alternative methods are required.

A basic approach here is computing the \href{https://en.wikipedia.org/wiki/Schatten_norm}{Schatten norm}, which serves as an upper bound for the spectral norm. Let the SVD of matrix $\boldsymbol{W}$ be $\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}$, with singular values $\sigma_1\geq\sigma_2\geq\cdots\geq\sigma_m$. The Schatten-$p$ norm is defined as:
\begin{equation}
\Vert\boldsymbol{W}\Vert_{S,p} = \sqrt[p]{\sum_{i=1}^m \sigma_i^p} \quad\geq\quad \sigma_1
\end{equation}
The larger $p$, the closer it is to the true value. When $p$ is even, computation of the Schatten-$p$ norm is relatively feasible, primarily using the following identities:
\begin{gather}
\Vert\boldsymbol{W}\Vert_{S,2k} = \sqrt[2k]{\sum_{i=1}^m \sigma_i^{2k}} = \sqrt[2k]{\mathrm{tr}((\boldsymbol{W}^{\top}\boldsymbol{W})^k)} \label{eq:S-2k} \\
\Vert\boldsymbol{W}\Vert_{S,4k} = \sqrt[4k]{\sum_{i=1}^m \sigma_i^{4k}} = \sqrt[2k]{\Vert(\boldsymbol{W}^{\top}\boldsymbol{W})^k\Vert_{S,2}}\label{eq:S-4k}
\end{gather}
Note that the Schatten-2 norm is the Frobenius norm, computable via element-wise squares, sum, and a square root. Thus, these two identities show that once $(\boldsymbol{W}^{\top}\boldsymbol{W})^k$ is computed, $\Vert\boldsymbol{W}\Vert_{S,2k}$ and $\Vert\boldsymbol{W}\Vert_{S,4k}$ can be derived at low cost.

The first step $\boldsymbol{W}^{\top}\boldsymbol{W}$ has complexity $\mathcal{O}(nm^2)$, and repeated squaring has cost $\mathcal{O}(m^3)$ per step, so $\mathcal{O}(nm^2 + Tm^3)$ suffices to compute $(\boldsymbol{W}^{\top}\boldsymbol{W})^{2^T}$. In theory, this complexity far exceeds power iteration, and even when $n,m$ are large, the first step $\boldsymbol{W}^{\top}\boldsymbol{W}$ alone far exceeds power iteration. Therefore, if we only aim to estimate the spectral norm without requiring an upper bound, power iteration should be preferred.

However, matrix multiplication is typically highly parallelizable, so when $n,m$ are small or $n\gg m$, Schatten-$p$ norms do not necessarily become a computational bottleneck, making them worth considering. Alternatively, when a strict upper bound is required, this might be the most direct approach.

\section{Preserving Numeric Stability} \label{preserving-numeric-stability}
To compute $\Vert\boldsymbol{W}\Vert_{S,2^{T+2}}$, a naive implementation starts from $\boldsymbol{M} = \boldsymbol{W}^{\top}\boldsymbol{W}$, applies $\boldsymbol{M}\leftarrow \boldsymbol{M}^2$ repeatedly $T$ times to obtain $(\boldsymbol{W}^{\top}\boldsymbol{W})^{2^T}$, and then applies formula \eqref{eq:S-4k}. However, due to super-exponential operations, values quickly explode to NaN or collapse to zero.

In general, $\Vert\boldsymbol{W}\Vert_{S,2^{T+2}}$ itself does not numerically explode; the issue lies in computing $(\boldsymbol{W}^{\top}\boldsymbol{W})^{2^T}$ explicitly. The solution is to renormalize at each squaring step: $\boldsymbol{M}\leftarrow \boldsymbol{M}^2/\mathrm{tr}(\boldsymbol{M}^2)$, preventing numerical explosion and ensuring scaling remains tight, avoiding collapse to zero. Meanwhile, accumulate normalization factors in the log domain for final computation of $\Vert\boldsymbol{W}\Vert_{S,2^{T+2}}$.

The computation flow is summarized as follows:
\[
\boxed{
\begin{array}{@{}l@{}}
\text{Compute } \Vert\boldsymbol{W}\Vert_{S,2^{T+2}} \text{ as an upper bound for spectral norm} \\
\hline
\begin{array}{ll}
1: & \text{Initialize } \log S = \log\mathrm{tr}(\boldsymbol{W}^{\top}\boldsymbol{W}), \, \boldsymbol{M}=\frac{\boldsymbol{W}^{\top}\boldsymbol{W}}{ \mathrm{tr}(\boldsymbol{W}^{\top}\boldsymbol{W})} \\
2: & \textbf{For } t=1,2,\cdots,T \textbf{ do } \\
3: & \qquad \log S\leftarrow 2\log S + \log \mathrm{tr}(\boldsymbol{M}^2) \\
4: & \qquad \boldsymbol{M} \leftarrow \frac{\boldsymbol{M}^2}{ \mathrm{tr}(\boldsymbol{M}^2)} \\
5: & \text{Output } \exp\left(\frac{\log S + \log\Vert\boldsymbol{M}\Vert_F}{2^{T+1}}\right)
\end{array} \\
\hline
\end{array}
}
\]

A simple reference implementation:
\begin{lstlisting}[language=Python]
@jax.jit
def tr(w):
    return w.trace(axis1=-1, axis2=-2)[..., None, None]

@jax.jit(static_argnums=(1,))
def spec_norm_v3(w, T=5):
    m = (m := w.mT @ w) / (s := tr(m))
    ms_step = lambda i, ms: ((m := ms[0] @ ms[0]) / (s := tr(m)), 2 * ms[1] + jnp.log(s))
    m, logs = lax.fori_loop(0, T, ms_step, (m, jnp.log(s)))
    logf = 0.5 * jnp.log((m**2).sum(axis=[-1, -2], keepdims=True))
    return jnp.exp((logs + logf) / 2**(T + 1))
\end{lstlisting}

\section{Higher-Order Moments} \label{higher-order-moments}
To compute $\Vert\boldsymbol{W}\Vert_{S,2^{T+2}}$, we implicitly compute $\boldsymbol{W}^{\top}\boldsymbol{W}, \cdots, (\boldsymbol{W}^{\top}\boldsymbol{W})^{2^{T-1}}, (\boldsymbol{W}^{\top}\boldsymbol{W})^{2^T}$, meaning we also obtain $\Vert\boldsymbol{W}\Vert_{S,2}, \cdots, \Vert\boldsymbol{W}\Vert_{S,2^{T+1}}, \Vert\boldsymbol{W}\Vert_{S,2^{T+2}}$, but only use the last---again seeming somewhat ``wasteful.''

Can we, like with Krylov subspace methods in power iteration, leverage all these results to improve estimation accuracy? Indeed! \href{https://jiha-kim.github.io/posts/fast-spectral-norm-bounds-from-scaled-gram-moments/}{``Fast Tight Spectral-Norm Bounds''} provides a non-linear programming approach to refine estimation. Let us first consider a simple case: suppose we have computed the 2nd and 4th moments of singular values:
\begin{equation}
\sum_{i=1}^m \sigma_i^2 = S_2, \qquad \sum_{i=1}^m \sigma_i^4 = S_4
\end{equation}
Then from previous sections, $S_4^{1/4}$ is closer to the spectral norm than $S_2^{1/2}$, so we return $S_4^{1/4}$ as the upper bound and discard $S_2$. But is $S_2$ truly useless? Consider the case $m=2$: two equations, two unknowns---$\sigma_1,\sigma_2$ can theoretically be solved! When $m>2$, exact solutions aren't possible, but $\sigma_1$'s range can still be narrowed. Specifically, we have:
\begin{equation}
\frac{S_2 - \sigma_1^2}{m-1} = \frac{1}{m-1}\sum_{i=2}^m \sigma_i^2 \leq \sqrt{\frac{1}{m-1}\sum_{i=2}^m \sigma_i^4} = \sqrt{\frac{S_4 - \sigma_1^4}{m-1}}
\end{equation}
i.e., $(S_2 - \sigma_1^2)^2\leq (m-1)(S_4 - \sigma_1^4)$, which is essentially a univariate quadratic inequality, yielding:
\begin{equation}
\sigma_1 \leq \sqrt{\frac{S_2+\sqrt{(m-1)(mS_4-S_2^2)}}{m}}\label{eq:s2-s4}
\end{equation}
This gives an upper bound for $\sigma_1$ using $S_2$ and $S_4$, better than $S_4^{1/4}$. \href{https://jiha-kim.github.io/posts/fast-spectral-norm-bounds-from-scaled-gram-moments/}{``Fast Tight Spectral-Norm Bounds''} provides results using $S_2,S_4,S_6,S_8$, but it's more complex overall; interested readers can refer to the original paper. In theory, using more moments improves estimation, but in practice, solving increasingly complex non-linear systems reduces practical value.

\section{Limitations} \label{limitations}
Theoretically, result \eqref{eq:s2-s4} can be generalized to arbitrary $S_{2k}$ and $S_{4k}$ for more accurate upper bounds:
\begin{equation}
\sigma_1 \leq \sqrt[2k]{\frac{S_{2k}+\sqrt{(m-1)(mS_{4k}-S_{2k}^2)}}{m}}
\end{equation}
However, when $k$ is large, this has limited practical value, as it requires explicitly computing $S_{2k}$ and $S_{4k}$, which explode or collapse when $k$ is large---unrealistic in practice. A possible improvement is to separate out $S_{4k}^{1/4k}$:
\begin{equation}
\sqrt[2k]{\frac{S_{2k}+\sqrt{(m-1)(mS_{4k}-S_{2k}^2)}}{m}} = S_{4k}^{1/4k}\cdot\sqrt[2k]{\frac{S_{2k}/S_{4k}^{1/2}+\sqrt{(m-1)(m-S_{2k}^2/S_{4k})}}{m}}
\end{equation}
Let $S_{2k}/S_{4k}^{1/2}=e^{\epsilon}$, and approximate in the log domain:
\begin{equation}
\sqrt[2k]{\frac{S_{2k}/S_{4k}^{1/2}+\sqrt{(m-1)(m-S_{2k}^2/S_{4k})}}{m}} \approx 1 - \frac{\epsilon^2}{4k(m-1)}
\end{equation}
However, our goal is a spectral norm upper bound; ensuring the upper bound property during approximation seems complex. On the other hand, if we already compute large $k$, then $S_{4k}^{1/4k}$ is already highly accurate, so further improving precision via optimization offers diminishing returns.

\section{Conclusion} \label{conclusion}
This article briefly summarizes several ideas for estimating spectral norms. In practical applications, if only approximate monitoring of the spectral norm is needed, power iteration and its subspace-accelerated variants are usually sufficient. If a strict upper bound is required, computing Schatten norms can be considered. The two paths each extend a ``no waste'' philosophy---Krylov subspace techniques utilize iteration history, while nonlinear programming utilizes higher-order moment information---providing valuable algorithm design insights.

\begin{center}
\color{orange}
\emph{\textbf{Please cite the original URL when reposting:} \href{https://kexue.fm/archives/11736}{https://kexue.fm/archives/11736}}\\
\emph{\textbf{For detailed reposting guidelines:} \href{https://kexue.fm/archives/6508\#article-repost}{``Scientific Space FAQ''}}
\end{center}
\end{document}
