Multiply, normalize, repeat: the whole power method
The previous guide left us with a promise: since you cannot reliably solve the characteristic polynomial, you reach eigenvalues by iterating. The power method (also called power iteration) is the most stripped-down iteration that delivers, and its rule fits in one breath: pick any starting vector, multiply it by the matrix A, rescale it back to length 1, and repeat. Do that enough times and the vector stops moving — and where it settles is an eigenvector of A. There is no system to solve, no factorization, nothing but matrix-times-vector over and over.
choose q_0 (any nonzero vector, e.g. all ones)
for n = 0, 1, 2, ... :
z = A q_n # one matrix-vector product
q_{n+1} = z / ||z|| # rescale back to length 1
lambda_estimate = q_{n+1}^T A q_{n+1} # Rayleigh quotient
stop when q_{n+1} stops moving (residual small)Why on earth should multiplying by A again and again single out a special direction? Here is the picture. Suppose A has a full set of eigenvectors, and write your starting vector as a blend of them — some amount of each direction. Every time you multiply by A, each ingredient gets stretched by its own eigenvalue. The direction belonging to the largest eigenvalue (in absolute value) is stretched the most, so after many multiplications it overwhelms all the others, like the loudest voice eventually drowning out a choir. The rescaling step only keeps the length manageable; it never changes which direction is winning.
How fast, and when it stalls
The speed of the takeover is set by a single ratio. Call the two largest eigenvalues (by magnitude) lambda_1 and lambda_2. Each step, the unwanted lambda_2-direction shrinks relative to the wanted lambda_1-direction by the factor |lambda_2 / lambda_1|. So the error in the eigenvector behaves like (|lambda_2/lambda_1|)^n — geometric, or linear, convergence. If lambda_1 = 4 and lambda_2 = 1 the ratio is 1/4 and you gain roughly one decimal digit every two steps; lovely. But if lambda_1 = 4.01 and lambda_2 = 4.00, the ratio is 0.9975 and you crawl, needing thousands of steps for a few digits. This is the same spectral-ratio story that governed the stationary iterations of the linear-systems rung: convergence is linear, and its rate is a ratio of eigenvalues.
Be honest about the catches, because they are real and they bite. First, the power method finds only the single largest-magnitude eigenvalue — it is blind to all the others. Second, if the two largest are tied or nearly tied in magnitude, the ratio is near 1 and convergence is glacial. Third, a complex-conjugate pair of equal magnitude makes the iterate spin instead of settling, so it never converges at all. And finally, if you are unlucky and your starting vector happens to have zero component along the dominant eigenvector, in exact arithmetic the method could not see it — though in practice the tiny round-off of floating-point arithmetic sprinkles a nonzero component back in, and the method quietly recovers.
Inverse iteration: aim anywhere you like
The power method's blindness to all but the biggest eigenvalue is fixed by a clean piece of algebra. The eigenvectors of A are exactly the eigenvectors of (A - mu I)^{-1}, where mu is any number you choose as a target. But the eigenvalues are transformed: if A has eigenvalue lambda, then (A - mu I)^{-1} has eigenvalue 1/(lambda - mu). So the eigenvalue of A closest to your target mu becomes the LARGEST eigenvalue of (A - mu I)^{-1} — and the power method, which always grabs the largest, will now grab exactly the one you aimed at. This is inverse iteration (also inverse iteration with a shift), and it lets you fish out any eigenvalue, interior or smallest, just by choosing mu near it.
- Pick a target mu near the eigenvalue you want, and a starting vector q_0.
- Solve the linear system (A - mu I) z = q_n for z — this is the inverse multiply, but done by SOLVING, never by forming an actual inverse.
- Rescale: q_{n+1} = z / ||z||.
- Update the eigenvalue estimate with the Rayleigh quotient lambda = q_{n+1}^T A q_{n+1}, and stop when the residual ||A q - lambda q|| is tiny.
Two practical points keep this from going wrong. First, you do NOT compute (A - mu I)^{-1}; you SOLVE the system (A - mu I) z = q_n at each step. Factor A - mu I once with an LU factorization, then every iteration is a cheap forward-and-back substitution — the factor-once-solve-many habit from the linear-systems rung pays off again. Forming the explicit inverse would be slower and less accurate, the very inverse-avoidance lesson you already know. Second, the convergence speed is now |(lambda_target - mu)/(lambda_second - mu)|: pick mu very close to the eigenvalue you want and the numerator is tiny, so the method races.
An honest worry: A - mu I is nearly singular
Here is the objection that should be nagging you. To converge fast we want mu extremely close to a true eigenvalue lambda. But as mu approaches lambda, the matrix A - mu I approaches a singular matrix — its smallest singular value shrinks toward zero, so its condition number blows up. We spent the whole conditioning rung warning that solving a system with a huge condition number loses digits. Are we not sawing off the branch we sit on? It is a fair worry, and a beginner's instinct to back mu away from lambda 'for safety' is exactly the wrong move.
The resolution is one of the loveliest facts in numerical linear algebra. Yes, the system is ill-conditioned, so the solution z is computed with a large RELATIVE error. But the error points almost entirely along the very eigenvector you are chasing — because that is the direction the near-singular matrix amplifies most. So the 'wrong' part of z is mostly more of the right answer. After you rescale, the contamination washes out, and the inaccurate solve still hands you an excellent eigenvector. A backward-stable solver gives the exact answer to a nearby system, and for this special problem 'nearby' is good enough: the ill-conditioning helps rather than hurts.
Rayleigh quotient iteration: change your aim every step
Inverse iteration converges fast when mu is close to the eigenvalue — so why keep mu fixed? Each step already produces a better eigenvector, and from a better eigenvector the Rayleigh quotient gives a better eigenvalue estimate. So feed it back: after each step, set the new shift mu to the current Rayleigh quotient and re-aim. As your eigenvector improves, your target marches toward the true eigenvalue, the numerator of the convergence ratio collapses, and the iteration accelerates relentlessly. This feedback loop is Rayleigh quotient iteration (RQI), and it is the spiritual cousin of Newton's method for eigenproblems.
The payoff is spectacular. For a symmetric matrix, RQI converges CUBICALLY near the answer: the number of correct digits roughly triples every step. Watch a typical run — the error shrinks like 1e-2, then 1e-6, then 1e-18 — three steps and you are at machine precision. That cubic burst is exactly the eps-to-eps^2 magic of the Rayleigh quotient compounding with the inverse-iteration speed-up. The catch is the same as before in reverse: because mu changes each step, you can no longer factor A - mu I just once; you pay for a fresh LU factorization every iteration. But when only two or three iterations are needed, that is a bargain nobody refuses.