Why Naive Gram-Schmidt + Power Iteration Fails

Here is a plan that sounds completely reasonable. You want the top few eigenvectors of a symmetric matrix. Power iteration finds the biggest one: multiply a random vector by $A$ repeatedly and it turns toward the dominant eigenvector. So: find $v_1$ that way, project it out of a fresh start vector (one Gram–Schmidt step — subtract the component along $v_1$), iterate again to get $v_2$, and repeat. Every piece is textbook. The assembly fails three separate ways, and every failure is measurable.

Failure 1: clustered eigenvalues make it crawl

Power iteration converges at the rate $|{\lambda }_2/{\lambda }_1|$ per step — each multiply shrinks the unwanted components by that factor. When eigenvalues are well separated, fine. When they cluster, arithmetic happens: at a ratio of 0.999, one correct digit costs $\ln 10/\ln (1/0.999)\approx 2301$ iterations. Ten digits, twenty-three thousand matrix multiplies — per eigenvector. Clustered spectra are not exotic; they are what large physical matrices look like.

Failure 2: classical Gram–Schmidt does not produce an orthogonal basis

Extend the plan to many vectors and you are running classical Gram–Schmidt (CGS): orthogonalize each new vector against all the accepted ones by subtracting projections computed from the original vector. In exact arithmetic, flawless. In floating point, each subtraction leaves a small residue of the directions you removed, and CGS never looks back to clean it up. Measured on a 40-column basis with condition number $10^{12}$ ($\Vert I-Q^{T}Q\Vert$ should be 0 for an orthonormal $Q$):

classical GS   ||I - QtQ|| = 1.2e+01     <- not even slightly orthogonal
modified GS    ||I - QtQ|| = 5.4e-05
Householder    ||I - QtQ|| = 5.2e-15     <- machine precision

Read the first line again: 12, not 12×10⁻⁶. The "orthonormal" basis CGS returns on a hard problem has columns pointing substantially along each other. Modified GS — mathematically identical, but subtracting each projection from the working vector immediately, so later projections see the corrected vector — buys ten orders of magnitude. Householder QR, which never forms projections at all, is exact to machine precision. Same equations on paper; the order of operations is the whole difference.

Failure 3: deflate-and-trust lets the old eigenvector grow back

The quiet killer. Suppose you found $v_1$ to accuracy $10^{-4}$ — good enough, surely — projected it out of your start vector once, and iterated toward $v_2$. The projection left a $10^{-4}$ sliver of true $v_1$ behind, and every multiply by $A$ grows that sliver by ${\lambda }_1/{\lambda }_2$ relative to everything else. With ${\lambda }_1/{\lambda }_2=1.0/0.98$, the sliver overtakes in roughly $\ln (10^4)/\ln (1.0204)\approx 455$ steps. Watched live:

k=   50   overlap with v1 = 0.0000   overlap with v2 = 1.0000
k=  200   overlap with v1 = 0.0005   overlap with v2 = 1.0000   <- looks converged
k=  500   overlap with v1 = 0.2061   overlap with v2 = 0.9785
k= 1000   overlap with v1 = 1.0000   overlap with v2 = 0.0002   <- wrong answer, confidently

At 200 iterations the method reports $v_2$ to four digits. Run it longer — the usual instinct when you want more accuracy — and it silently converges back to $v_1$. More iterations make the answer worse. The fix is not a better first estimate; it is re-projecting on every iteration, which keeps the sliver pinned at roundoff size instead of letting it compound.

What to do instead

Each failure has a standard repair: re-orthogonalize every step (and twice when it matters — "twice is enough" is a theorem-grade slogan in this business), prefer modified GS or Householder when building bases, and use shifts or subspace iteration to break through clustered eigenvalues. But assemble all the repairs and you discover you have built something that already has a name: Lanczos — a Krylov method that gets the deflation for free from three-term orthogonality and treats the floating-point decay of that orthogonality as a first-class concern rather than a surprise. The naive plan is not wrong to want what it wants; it is an early draft of the right algorithm, missing exactly the parts that fail above.

Reproduce it

import numpy as np
np.random.seed(1)

def cgs(A):
    Q = np.zeros_like(A)
    for j in range(A.shape[1]):
        v = A[:,j] - Q[:,:j] @ (Q[:,:j].T @ A[:,j])   # projections from ORIGINAL vector
        Q[:,j] = v / np.linalg.norm(v)
    return Q

def mgs(A):
    A = A.copy(); Q = np.zeros_like(A)
    for j in range(A.shape[1]):
        Q[:,j] = A[:,j] / np.linalg.norm(A[:,j])
        for k in range(j+1, A.shape[1]):
            A[:,k] -= (Q[:,j] @ A[:,k]) * Q[:,j]       # correct the WORKING vectors now
    return Q

n = 40
U,_ = np.linalg.qr(np.random.randn(n,n)); V,_ = np.linalg.qr(np.random.randn(n,n))
A = U @ np.diag(np.logspace(0, -12, n)) @ V.T          # condition number 1e12
for name, Q in [("CGS", cgs(A)), ("MGS", mgs(A)), ("Householder", np.linalg.qr(A)[0])]:
    print(name, np.linalg.norm(np.eye(n) - Q.T @ Q))