Stability of the Lanczos algorithm on matrices with regular spectral distributions

TL;DR

This work analyzes the Lanczos algorithm’s stability when run on problems whose eigenvector empirical spectral distribution is close to a regular reference measure with suitable orthogonal polynomials. By extending Knizhnerman’s modified Chebyshev moment framework and employing a backwards-stability construction to a nearby measure $\mu_*$, the authors show backwards stability and, under a square-root density condition at the spectrum edges, forward stability for many large random matrix models. The analysis connects moment perturbations to Jacobi-recurrence coefficients via a Riemann–Hilbert approach and yields explicit, polynomial-in-$k$ bounds under reasonable regularity assumptions. Practically, this justifies using random matrices to test Lanczos-based methods and explains why the algorithm can behave nearly deterministically in finite-precision arithmetic when $k$ is small relative to $N$.

Abstract

We study the stability of the Lanczos algorithm run on problems whose eigenvector empirical spectral distribution is near to a reference measure with well-behaved orthogonal polynomials. We give a backwards stability result which can be upgraded to a forward stability result when the reference measure has a density supported on a single interval with square root behavior at the endpoints. Our analysis implies the Lanczos algorithm run on many large random matrix models is in fact forward stable, and hence nearly deterministic, even when computations are carried out in finite precision arithmetic. Since the Lanczos algorithm is not forward stable in general, this provides yet another example of the fact that random matrices are far from "any old matrix", and care must be taken when using them to test numerical algorithms.

Figures (5)

• Figure 1: Here $(\mathbf{A},\mathbf{b})$ corresponds to a $2000\times 2000$ random matrix, drawn from the Gaussian orthogonal ensemble (see \ref{['sec:wigner']}), and independent vector. In the large $N$ limit, the VESD of matrices drawn from this ensemble converge to the semicircle distribution on $[-1,1]$ (density $\propto \sqrt{1-x^2}$). Therefore the Lanczos coefficients $\alpha_i$ and $\beta_i$ from the "exact" computation (with reorthogonalization in quadruple precision arithmetic) respectively converge to $1/2$ and $0$; i.e. the Lanczos algorithm exhibits deterministic behavior. In our particular experiment we observe fluctuations on the order of $10^{-2}$ around the limiting values due to finite $N$ effects. Remarkably, the coefficients $\overline{\alpha}_n$ and $\overline{\beta}_n$ output by the Lanczos algorithm run in single precision floating point arithmetic without reorthogonalization are within the unit roundoff ($\approx 10^{-7}$) of $\alpha_n$ and $\beta_n$, at least while $n$ is sufficiently small; i.e. the algorithm is forward stable.
• Figure 2: We use the same $2000\times 2000$ random matrix $\mathbf{A}$ and fixed vector $\mathbf{b}$ from \ref{['fig:motivating_experiment']}. After running the Lanczos algorithm on $(\mathbf{A},\mathbf{b})$ in single precision finite precision arithmetic without reorganization, we use \ref{['eqn:new_measure']} and \ref{['rem:backwards_finite']} to construct a slightly perturbed $\mathbf{b}_*$. Lanczos run on $(\mathbf{A},\mathbf{b}^*)$ "exactly" (with reorthogonization in quadruple precision floating point arithmetic) produces a nearly identical output as the original single precision finite precision computation. For reference, $\|\mathbf{b} - \mathbf{b}_*\|\approx 3.6\cdot 10^{-5}$ is only a few orders of magnitude above the machine precision in which the original computation was carried out.
• Figure 3: Maximum value of orthogonal polynomial $p_n(\cdot\,;\ref{['eqn:VESD']})$ over $[-1,1]$ and the maximum value in the $N \to \infty$ limit (), where \ref{['eqn:VESD']} is drawn from the same random matrix model as in \ref{['fig:motivating_experiment', 'fig:GOE_backwards']}. For each $n,N$, the violin plot gives the distribution of $\|\ref{['eqn:op']}(\cdot\,;\ref{['eqn:VESD']})\|_{[-1,1]}$, with the 5%, 50%, and 95% quantiles marked explicitly. Note that for $k$ growing sufficiently slow with $N$, the maximum value of $p_n$ has polynomial growth for all $n\leq k$.
• Figure 4: Output of Lanczos run on $(\mathbf{A}_N,\mathbf{b}_N)$ in single precision arithmetic, where $\mathbf{A}_N$ is a GOE matrix of size $N$ and $\mathbf{b}_N$ is an independent vector; see \ref{['fig:motivating_experiment']} for more details.
• Figure 5: Error of Lanczos used to solve the system $\mathbf{A}\mathbf{x} = \mathbf{b}$ in single precision floating point arithmetic. Here $\mathbf{b}_N$ is proportional to the all ones vector and $\mathbf{A}_N = N^{-1} \mathbf{X}\mathbf{X}^\mathsf{T}$ where the entries of $\mathbf{X}$ are either iid standard normal random variables or iid Rademacher random variables ($\pm 1$ with equal probability). For each $n,N$, the violin plot gives the distribution of error and the 5%, 50%, and 95% quantiles are marked. Notice the convergence to the "deterministic" behavior () as $N$ increases, at least until the maximal accuracy is reached.

Theorems & Definitions (49)

• Definition 1.1
• Remark 2.1
• Definition 2.2
• Proposition 2.3: informal; see paige_70paige_80
• Remark 2.4
• Theorem 3.1: Backwards stability
• Corollary 3.2
• Lemma 3.3
• Corollary 3.5
• Corollary 3.6