Lanczos Iteration and Shift-Invert for Extreme Eigenvalues

Computing all eigenvalues of a symmetric matrix costs O(n^3) and produces n numbers, most of which are irrelevant. Applications usually want a handful: the lowest vibrational mode, the largest principal component, the eigenvalue closest to a stability boundary. For those, the matrix is often available only as a product, never as an array, and the only operation possible is multiplication by a vector.

The Lanczos iteration fits that constraint exactly. It builds a small tridiagonal matrix whose eigenvalues approximate the extreme eigenvalues of the original, using nothing but matrix-vector products and a three-term recurrence, and it does so with storage that does not grow with the number of iterations.

From Arnoldi to a three-term recurrence

The Arnoldi process builds an orthonormal basis v_1, ..., v_k of the Krylov subspace by orthogonalizing each new product A v_j against all previous basis vectors. The coefficients form an upper Hessenberg matrix H_k of size k x k, and the eigenvalues of H_k, called Ritz values, approximate eigenvalues of A.

For a symmetric matrix, H_k is symmetric as well as Hessenberg, which forces it to be tridiagonal. The orthogonalization against all previous vectors collapses to an orthogonalization against the last two, because symmetry makes the inner products with older vectors vanish in exact arithmetic. What remains is a recurrence with three terms:

beta_j v_{j+1} = A v_j - alpha_j v_j - beta_{j-1} v_{j-1}

with alpha_j = v_j^T A v_j and beta_j the norm of the residual after the projection. Each step costs one matrix-vector product, two inner products, and a few vector updates, and it stores three vectors of length n plus the two diagonal and subdiagonal sequences.

#include <math.h>
#include <stdlib.h>

static double dot(int n, const double *a, const double *b)
{
    double s = 0.0;
    for (int i = 0; i < n; i++) s += a[i] * b[i];
    return s;
}

/* Lanczos with full reorthogonalization. V holds k+1 basis vectors of
   length n in column-major order. alpha and beta receive the tridiagonal. */
void lanczos(int n, int k, void (*matvec)(const double *, double *),
             double *V, double *alpha, double *beta)
{
    double *w = malloc((size_t)n * sizeof *w);
    double *h = malloc((size_t)(k + 1) * sizeof *h);
    if (!w || !h) { free(w); free(h); return; }
    for (int i = 0; i < n; i++) V[i] = 0.0;
    V[0] = 1.0;                                   /* start vector: e_1 */
    beta[0] = 0.0;

    for (int j = 0; j < k; j++) {
        double *vj = V + (size_t)j * n;
        matvec(vj, w);
        if (j > 0)
            for (int i = 0; i < n; i++)
                w[i] -= beta[j] * V[(size_t)(j - 1) * n + i];
        alpha[j] = dot(n, w, vj);
        for (int i = 0; i < n; i++) w[i] -= alpha[j] * vj;

        /* Reorthogonalize against the whole basis; O(j*n) work per step. */
        for (int m = 0; m <= j; m++) {
            h[m] = dot(n, w, V + (size_t)m * n);
            for (int i = 0; i < n; i++) w[i] -= h[m] * V[(size_t)m * n + i];
        }
        /* Fold the correction into the diagonal so the tridiagonal stays
           consistent with the basis that was actually built. */
        alpha[j] += h[j];

        double nw = sqrt(dot(n, w, w));
        beta[j + 1] = nw;
        if (nw <= 1e-14) {                        /* invariant subspace found */
            for (int i = 0; i < n; i++) V[(size_t)(j + 1) * n + i] = 0.0;
            break;
        }
        for (int i = 0; i < n; i++)
            V[(size_t)(j + 1) * n + i] = w[i] / nw;
    }
    free(w); free(h);
}

The reorthogonalization block is the honest version of the algorithm. Without it, the recurrence is cheaper but the basis loses orthogonality as the Ritz values converge, and the loss has a specific and misleading symptom: duplicate copies of already-converged eigenvalues appear in the tridiagonal spectrum. Those copies are not eigenvalues of the matrix; they are an artifact of the arithmetic.

What the subspace can resolve

The Ritz values converge to eigenvalues of A, but not uniformly. The extreme eigenvalues converge first, and the rate depends on the gap between them and the rest of the spectrum.

The mechanism is polynomial approximation, as it is for the linear solvers described in Krylov Subspace Iteration and Why GMRES Restarts. After k steps, the residual for a Ritz pair is bounded by the best uniform approximation of zero on the spectrum by a degree-k polynomial that takes the value one at the Ritz value. Chebyshev polynomials are optimal on an interval, and their convergence rate is governed by the ratio of the interval’s width to its distance from the point being approximated.

The consequence is that an eigenvalue separated from the rest by a large gap converges quickly, in a number of iterations that depends weakly on n. An eigenvalue in the middle of a dense cluster converges slowly, because no low-degree polynomial can be small on the cluster while being large at one point inside it. Computing interior eigenvalues by plain Lanczos is therefore impractical, no matter how many iterations are run.

The same analysis explains why the method is used for the largest and smallest eigenvalues, and why it is combined with a transformation when interior ones are needed.

Shift-invert changes which eigenvalues are extreme

Multiplying by (A - sigma I)^{-1} maps an eigenvalue lambda to 1/(lambda - sigma). Eigenvalues close to sigma map to values of large magnitude; eigenvalues far away map to values near zero. Applying Lanczos to the inverted operator therefore converges to the eigenvalues of A nearest sigma, because those are the extreme ones of the transformed problem.

The transformation is exact in exact arithmetic and preserves symmetry when A is symmetric and sigma is real, since the inverse of a symmetric matrix is symmetric. The cost is a factorization of A - sigma I, which is the same sparse direct problem described in Multifrontal Factorization and the Dense Front Matrix, followed by one triangular solve per iteration. For a single eigenvalue, the factorization dominates; for many eigenvalues near the same shift, it is amortized.

There are three practical complications. The first is that A - sigma I may be indefinite or singular. If sigma is an exact eigenvalue, the matrix is singular and the factorization fails; a small perturbation of the shift, or a factorization that handles the singularity, is needed. If the matrix is indefinite, the factorization must use a symmetric indefinite kernel with 2 x 2 pivots, not the Cholesky kernel.

The second is that the inverted operator is not available as a matrix, only as a solve. The product is exact up to the accuracy of the triangular solves, which is a backward stable computation but still introduces error at the level of the condition number of A - sigma I. A shift very close to an eigenvalue makes that matrix ill conditioned, and the computed products lose accuracy in proportion.

The third is that the transformed problem’s extreme eigenvalues are clustered when the original eigenvalues near sigma are clustered. The transformation maps a cluster to a cluster, so the convergence rate for distinguishing eigenvalues within the cluster is not improved by the shift; only their separation from the rest of the spectrum is.

Restarting, locking, and the cost of memory

Lanczos needs three vectors regardless of k, but the Ritz vectors need the whole basis. Computing an eigenvector requires combining all k Lanczos vectors with the components of the eigenvector of the tridiagonal matrix, so the basis must be retained, at k vectors of length n.

For k in the hundreds and n in the millions, that storage is substantial. The standard remedy is to restart: keep a small number of converged or promising Ritz vectors, discard the rest of the basis, and continue from the reduced set. Thick restart keeps several vectors rather than one, which preserves the convergence of the remaining eigenvalues; keeping only the best one, as in a naive restart, forces the iteration to rediscover information it already had.

Locking is the complementary technique. Once a Ritz pair has converged to a tolerance, it is deflated: the vector is stored separately, and subsequent iterations are orthogonalized against it so that they cannot converge to it again. This is what prevents the duplicate eigenvalues that the unreorthogonalized recurrence produces.

Full reorthogonalization costs O(k n) work per step, which for large k dominates the O(nnz) cost of the product. Selective reorthogonalization, which orthogonalizes only when the loss of orthogonality exceeds a threshold, recovers most of the stability at a fraction of the cost. The threshold is compared against the norm of the current vector, and the criterion is a variant of the “twice is enough” rule used for Gram-Schmidt.

Singular values without the cross-product

The singular values of a rectangular matrix A are the square roots of the eigenvalues of A^T A, and computing them that way squares the condition number, exactly as forming the normal equations does. A matrix with kappa = 1e8 produces a cross-product with kappa = 1e16, and the small singular values are lost to rounding before the eigensolver sees them.

The stable route is Golub-Kahan bidiagonalization, which applies Lanczos implicitly to the symmetric matrix [0 A; A^T 0]. Its eigenvalues are plus and minus the singular values of A, and its eigenvectors combine the left and right singular vectors. The iteration works with two sequences of vectors, one of length m and one of length n, and produces a bidiagonal matrix whose singular values approximate those of A. The conditioning of the bidiagonal problem is that of A, not its square.

The cost per step is one product with A and one with A^T, which is twice the work of a symmetric product but still linear in the number of nonzeros. The storage is a small number of vectors of each length, and the same restart and locking machinery applies. This is the algorithm behind the standard sparse singular value routines, and it is the reason a request for the largest few singular values of a large sparse matrix does not require forming A^T A.

The choice between the two formulations is not a matter of taste. Whenever the small singular values matter, which is precisely when the matrix is ill conditioned and the answer is interesting, the cross-product formulation has already destroyed them. The bidiagonalization keeps them, at twice the arithmetic per step.