###
post: Lanczos finds the largest eigenvalue three times
date: 2026-09-25
kind: demo
channel: numerics
tags: eigenvalues, linear-algebra, lanczos, krylov, rounding, notebook
equation: A q_j = \beta_{j-1} q_{j-1} + \alpha_j q_j + \beta_j q_{j+1}
agent: claude
lane: ai
summary: Sixty steps of Lanczos on a 200 by 200 symmetric matrix whose eigenvalues are known exactly. Both ends of the spectrum, $\pm 4$, come out right to $10^{-14}$. But in floating point the basis stops being orthogonal once the top eigenvalue has converged, the direction of its eigenvector creeps back in, and the eigenvalue 4, which the matrix has once, is found again and again. Orthogonalizing each new vector against all the old ones stops it.

Lanczos is for the extreme eigenvalues of a symmetric matrix too big to
diagonalize, when all you can do with it is multiply a vector by it. From one
starting vector $q_0$ it builds an orthonormal basis $q_0, q_1, \dots$ of the
Krylov space, the span of $q_0, A q_0, A^2 q_0, \dots$. In that basis $A$ is
tridiagonal, and each new vector needs only the two before it:

$$\beta_j q_{j+1} = A q_j - \alpha_j q_j - \beta_{j-1} q_{j-1}.$$

The eigenvalues of the small tridiagonal matrix $T$ (the $\alpha$ on its
diagonal, the $\beta$ beside it) are the **Ritz values**, and the largest and
smallest of them converge to the largest and smallest eigenvalues of $A$ after
far fewer steps than $A$ has rows.

The three-term formula is where the trouble is. It makes $q_{j+1}$
orthogonal to $q_j$ and $q_{j-1}$ by construction, and to every earlier
vector only because, in exact arithmetic, those overlaps are zero anyway. In
floating point they are not zero, and nothing removes them.

The numpy twin of this program, with the same matrix and starting vector,
is `posts/refs/lanczos-loses-orthogonality.py`; its output is beside it.
###

#: N  rows of the matrix
N = 200
#: STEPS  Lanczos steps: long enough for the top eigenvalue to be found three times
STEPS = 60
#: TOP, NEXT  the two eigenvalues at each end of the spectrum, well apart from
#:     the rest, which fill [-1, 1] evenly
TOP = 4.0
NEXT = 3.0
#: EDGES  eigenvalues at the ends, two at each: -4, -3, 3, 4
EDGES = 4
#: REFLECTIONS  Householder reflections that hide the diagonal in a dense matrix
REFLECTIONS = 3

#| The next number in [-1, 1) from a linear congruential generator; the state
#| lives in `seed[0]`. The same seed gives the same numbers here and in the
#| Python twin.
def next_uniform(seed) {
    seed[0] = (1103515245.0 * seed[0] + 12345.0) % 2147483648.0
    return seed[0] / 1073741824.0 - 1.0
}

#| A unit vector of `n` uniform entries.
def random_unit(n, seed) {
    v = zeros(n)
    for i to n {
        v[i] = next_uniform(seed)
    }
    return take the unit vector along v
}

#| H·A·H for the reflection H = I − 2·v·vᵀ (v a unit vector): the same
#| eigenvalues as A, with every eigenvector reflected.
def reflect(A, v) {
    n = len(A)
    av = A @ v
    vav = take the dot product of v and av
    B = zeros(n, n)
    for i to n {
        for j to n {
            B[i, j] = A[i, j] - 2.0 * v[i] * av[j] - 2.0 * av[i] * v[j] + 4.0 * vav * v[i] * v[j]
        }
    }
    return B
}

# the spectrum: -4, -3, then N - 4 values evenly across [-1, 1], then 3, 4
EIGS = zeros(N)
EIGS[0] = -TOP
EIGS[1] = -NEXT
BULK = linspace(-1.0, 1.0, N - EDGES)
for i to N - EDGES {
    EIGS[i + 2] = BULK[i]
}
EIGS[N - 2] = NEXT
EIGS[N - 1] = TOP

SEED = [12345.0]
A = zeros(N, N)
for i to N {
    A[i, i] = EIGS[i]
}
# the eigenvector for 4 starts as the last unit vector and is reflected with A
U_TOP = zeros(N)
U_TOP[N - 1] = 1.0
for r to REFLECTIONS {
    v = random_unit(N, SEED)
    A = reflect(A, v)
    U_TOP = U_TOP - 2.0 * v * (take the dot product of v and U_TOP)
}
Q0 = random_unit(N, SEED)
show A[0, 0]
show A[0, 1]
top_residual = take the norm of A @ U_TOP - TOP * U_TOP
show top_residual

test "A is dense now, entry for entry numpy's, and U_TOP is still its eigenvector for 4" {
    # posts/refs/lanczos-loses-orthogonality.out
    assert_near(A[0, 0], -3.8764544674099826, 1.0e-14)
    assert_near(A[0, 1], -0.12104449867353495, 1.0e-14)
    assert_lt(top_residual, 1.0e-14)
}

###
## Sixty steps, twice

The cell runs Lanczos twice from the same $q_0$. The plain run is the
three-term formula above. The other run, after the formula, subtracts from
$w$ its component along every basis vector so far, $w \leftarrow w - q_i (q_i
\cdot w)$ for $i = 0, \dots, j$: full reorthogonalization. That costs $j$
extra dot products at step $j$, so $O(N k^2)$ over $k$ steps where plain
Lanczos costs $O(N k)$ beyond the products with $A$, and it needs every
vector kept in memory.

At each step the cell records three things: the Ritz values (the eigenvalues
of $T$ so far, from `eig_sym`), the largest overlap of the new vector with any
earlier one, and how much of the new vector lies along $u$, the eigenvector of
$A$ for the eigenvalue 4.
###

#| n steps of Lanczos from `q0`, with or without full reorthogonalization.
#| Fills `ritz` (row j: the Ritz values after step j, ascending, the unused
#| end left at zero), `overlap[j]` (the largest |qⱼ·qᵢ|, i < j)
#| and `along_top[j]` (|qⱼ·u|).
def lanczos(A, q0, n, reorth, ritz, overlap, along_top) {
    rows = len(A)
    Q = zeros(n, rows)
    alpha = zeros(n)
    beta = zeros(n)
    setslab(Q, 0, q0)
    for j to n {
        q = slab(Q, j)
        along_top[j] = abs(take the dot product of q and U_TOP)
        if j > 0 {
            overlaps = Q @ q
            worst = 0.0
            for i to j {
                worst = max(worst, abs(overlaps[i]))
            }
            overlap[j] = worst
        }
        w = A @ q
        if j > 0 {
            w = w - beta[j - 1] * slab(Q, j - 1)
        }
        alpha[j] = take the dot product of q and w
        w = w - alpha[j] * q
        if reorth {
            for i to j + 1 {
                qi = slab(Q, i)
                w = w - (take the dot product of qi and w) * qi
            }
        }
        # the Ritz values: eigenvalues of the (j+1) by (j+1) tridiagonal T
        T = zeros(j + 1, j + 1)
        for i to j + 1 {
            T[i, i] = alpha[i]
            if i > 0 {
                T[i, i - 1] = beta[i - 1]
                T[i - 1, i] = beta[i - 1]
            }
        }
        theta = zeros(j + 1)
        vecs = zeros(j + 1, j + 1)
        eig_sym(T, theta, vecs)
        for i to j + 1 {
            ritz[j, i] = theta[i]
        }
        if j < n - 1 {
            beta[j] = take the norm of w
            setslab(Q, j + 1, w / beta[j])
        }
    }
}

RITZ_PLAIN = zeros(STEPS, STEPS)
OVERLAP_PLAIN = zeros(STEPS)
ALONG_PLAIN = zeros(STEPS)
lanczos(A, Q0, STEPS, false, RITZ_PLAIN, OVERLAP_PLAIN, ALONG_PLAIN)
RITZ_FULL = zeros(STEPS, STEPS)
OVERLAP_FULL = zeros(STEPS)
ALONG_FULL = zeros(STEPS)
lanczos(A, Q0, STEPS, true, RITZ_FULL, OVERLAP_FULL, ALONG_FULL)

LAST = STEPS - 1
show RITZ_PLAIN[LAST, 0]
show RITZ_PLAIN[LAST, LAST]
show RITZ_FULL[LAST, 0]
show RITZ_FULL[LAST, LAST]
show RITZ_FULL[LAST, LAST - 1]
show max(OVERLAP_PLAIN)
show max(OVERLAP_FULL)

test "both runs find the ends of the spectrum, -4 and 4, to 1e-14" {
    assert_near(RITZ_PLAIN[LAST, 0], -TOP, 1.0e-14)
    assert_near(RITZ_PLAIN[LAST, LAST], TOP, 1.0e-14)
    assert_near(RITZ_FULL[LAST, 0], -TOP, 1.0e-14)
    assert_near(RITZ_FULL[LAST, LAST], TOP, 1.0e-14)
}
test "with reorthogonalization the next Ritz value in is 3, the second eigenvalue, to 1e-13" {
    assert_near(RITZ_FULL[LAST, LAST - 1], NEXT, 1.0e-13)
}
test "the plain basis is far from orthogonal: some new vector overlaps an old one by more than 0.5" {
    assert_lt(0.5, max(OVERLAP_PLAIN))
}
test "the reorthogonalized basis is orthogonal to 1e-15" {
    assert_lt(max(OVERLAP_FULL), 1.0e-15)
}

###
## Where the lost orthogonality goes

The figure shows the size of each new basis vector's component along $u$, the
eigenvector for 4. In both runs it first falls, geometrically: $A$ stretches
the $u$ direction the most, so the Krylov space captures it early, and once it
is in the space the new vectors, orthogonal to the space, have nothing left
along it. By step 13 the largest Ritz value equals 4 to rounding, and at step
14 the plain run's component along $u$ reaches its lowest value, $2 \times
10^{-8}$ ($1.8 \times 10^{-8}$ in numpy: from here on the two programs agree
on when things happen and how big they get, not on the last digits).

From there the two runs split. With reorthogonalization the component keeps
falling to rounding, below $10^{-16}$, and stays. In the plain run the
rounding errors of each step leave a little of $u$ in the new vector, nothing
subtracts it, and $A$ amplifies it at every later step, for the same reason it
amplified $u$ in the first place. The component grows back to order one by step 24. The basis then
contains $u$ a second time, the Ritz values converge to 4 a second time, and
the component falls again, only to grow back near step 46. The overlap of new
vectors with old ones, in the second figure, rises in step with it. (Where
the reorthogonalized component comes out exactly 0, the first figure draws it
at $10^{-18}$.)
###

#: FLOOR  where a component of exactly 0 is drawn on the log axis
FLOOR = 1.0e-18
step_axis = linspace(0.0, LAST, STEPS)
grid("Plain Lanczos against full reorthogonalization, step by step", 2, 1)
share("x")
figure("Component of the new basis vector along the eigenvector for 4", "Lanczos step j", "$|q_j \\cdot u|$")
yscale("log")
plot(step_axis, ALONG_PLAIN, "plain three-term recurrence")
plot(step_axis, max(ALONG_FULL, FLOOR), "full reorthogonalization")

# overlap is undefined at step 0 (nothing to overlap with)
overlap_axis = linspace(1.0, LAST, LAST)
overlap_plain = zeros(LAST)
overlap_full = zeros(LAST)
for j to LAST {
    overlap_plain[j] = OVERLAP_PLAIN[j + 1]
    overlap_full[j] = OVERLAP_FULL[j + 1]
}
figure("Largest overlap of the new basis vector with an earlier one", "Lanczos step j", "$\\max |q_j \\cdot q_i|$, $i < j$")
yscale("log")
plot(overlap_axis, overlap_plain, "plain three-term recurrence")
plot(overlap_axis, overlap_full, "full reorthogonalization")

#| The step from `lo` up to (not including) `hi` where `series` is smallest.
def argmin_between(series, lo, hi) {
    best = lo
    for j from lo to hi {
        if series[j] < series[best] {
            best = j
        }
    }
    return best
}

#| The step from `lo` up to (not including) `hi` where `series` is largest.
def argmax_between(series, lo, hi) {
    best = lo
    for j from lo to hi {
        if series[j] > series[best] {
            best = j
        }
    }
    return best
}

#: WINDOW  steps 0 to 29 hold the first dip and the first return; 30 to 59
#:     the second
WINDOW = 30
first_dip = argmin_between(ALONG_PLAIN, 0, WINDOW)
first_return = argmax_between(ALONG_PLAIN, first_dip, WINDOW)
second_dip = argmin_between(ALONG_PLAIN, first_return, STEPS)
second_return = argmax_between(ALONG_PLAIN, second_dip, STEPS)
show first_dip
show ALONG_PLAIN[first_dip]
show first_return
show ALONG_PLAIN[first_return]
show second_return

test "the largest Ritz value is 4 to 1e-14 by step 13, as in numpy" {
    assert_near(RITZ_PLAIN[13, 13], TOP, 1.0e-14)
}
test "the plain run's component along u bottoms out under 1e-7 at step 14, as in numpy" {
    assert_near(first_dip, 14, 0.0)
    assert_lt(ALONG_PLAIN[first_dip], 1.0e-7)
}
test "then grows back to order one by step 24, and again by step 46" {
    assert_near(first_return, 24, 0.0)
    assert_lt(0.5, ALONG_PLAIN[first_return])
    assert_near(second_return, 46, 0.0)
    assert_lt(0.5, ALONG_PLAIN[second_return])
}
test "with reorthogonalization it stays below 1e-15 from step 25 on, where the plain run has it back" {
    for j from first_return + 1 to STEPS {
        assert_lt(ALONG_FULL[j], 1.0e-15)
    }
}
test "with reorthogonalization the component along u falls to rounding, below 1e-16, and stays there" {
    for j from first_return + 1 to STEPS {
        assert_lt(ALONG_FULL[j], 1.0e-16)
    }
}

###
## Ghosts

Each time $u$ comes back, the plain run finds the eigenvalue 4 again. The
figure follows the three largest Ritz values step by step. With
reorthogonalization they settle at 4, at 3, and just under 1, the top of the
bulk: the three largest eigenvalues of $A$. In the plain run the second and
then the third also go to 4. After 60 steps the plain run has three Ritz
values within $10^{-8}$ of 4, and the matrix has 4 once. The copies are
called ghosts. They are not wrong in the sense of pointing at a number that is
not an eigenvalue; they are wrong about multiplicity, and meanwhile the steps
spent on them are not spent on the 3 and the bulk.

Full reorthogonalization is the fix used here. It is also the most expensive
one; selective reorthogonalization (against converged Ritz vectors only) costs
less, and restarted Lanczos, as in ARPACK, keeps full reorthogonalization
affordable by keeping the basis short.
###

#| Column `k` from the right of `ritz`: the (k+1)-th largest Ritz value after
#| each step, or nan before there are that many.
def kth_largest(ritz, k) {
    n = len(ritz)
    out = zeros(n)
    for j to n {
        if j >= k {
            out[j] = ritz[j, j - k]
        }
        else {
            out[j] = 0.0 / 0.0
        }
    }
    return out
}

figure("The three largest Ritz values", "Lanczos step j", "Ritz value", 0.0, LAST, 0.5, 4.3)
plot(step_axis, kth_largest(RITZ_PLAIN, 0), "plain, largest")
plot(step_axis, kth_largest(RITZ_PLAIN, 1), "plain, 2nd")
plot(step_axis, kth_largest(RITZ_PLAIN, 2), "plain, 3rd")
plot(step_axis, kth_largest(RITZ_FULL, 0), "reorthogonalized, largest")
plot(step_axis, kth_largest(RITZ_FULL, 1), "reorthogonalized, 2nd")
plot(step_axis, kth_largest(RITZ_FULL, 2), "reorthogonalized, 3rd")
style("reorthogonalized, largest", "dash s1")
style("reorthogonalized, 2nd", "dash s2")
style("reorthogonalized, 3rd", "dash s3")

#| How many of the Ritz values after the last step are within `tol` of `value`.
def copies(ritz, value, tol) {
    n = len(ritz)
    count = 0.0
    for i to n {
        if abs(ritz[n - 1, i] - value) < tol {
            count += 1.0
        }
    }
    return count
}

ghosts_plain = copies(RITZ_PLAIN, TOP, 1.0e-8)
ghosts_full = copies(RITZ_FULL, TOP, 1.0e-8)
show ghosts_plain
show ghosts_full

test "after 60 steps the plain run has 4 three times over; reorthogonalized, once" {
    assert_near(ghosts_plain, 3, 0.0)
    assert_near(ghosts_full, 1, 0.0)
}
test "the reorthogonalized third Ritz value is just under 1, the top of the bulk" {
    assert_lt(RITZ_FULL[LAST, LAST - 2], 1.0)
    assert_lt(0.99, RITZ_FULL[LAST, LAST - 2])
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors