5.8 Iterative estimates for eigenvalues

Big idea

The characteristic polynomial is a definition, not an algorithm. For a matrix of any real size, expanding a symbolic determinant is out of reach, and even when the coefficients are available, polynomial roots move alarmingly under small changes in them. Practical eigenvalue software never forms the characteristic polynomial.

What it does instead is multiply. Take any starting vector, apply the matrix, apply it again, and keep going. Write the starting vector in an eigenvector basis and each term gets multiplied by its own eigenvalue at every step, so the term belonging to the largest eigenvalue grows fastest and eventually dominates. The direction of the iterate converges to the dominant eigenvector, and the amount the vector stretches at each step converges to the dominant eigenvalue.

Two adjustments make that idea usable. Rescale after every multiplication, or the numbers overflow or collapse. And to reach an eigenvalue other than the dominant one, shift and invert first: applying $(A - \alpha I)^{-1}$ makes the eigenvalue nearest $\alpha$ into the dominant one, so the same iteration finds it.

Decoder

The power method produces a sequence whose scaling factors converge to the strictly dominant eigenvalue, provided the initial vector has a nonzero component in the corresponding eigenspace.

Plainly: as long as your starting vector is not accidentally perpendicular to the answer, repeated multiplication finds the biggest eigenvalue. “Strictly dominant” means one eigenvalue is larger in absolute value than every other; a tie between $\lambda$ and $-\lambda$ breaks the method. The nonzero-component condition is nearly free in practice, since a random start has that component with probability one, and rounding error supplies it even when the start is unlucky.

Definitions and results

Strictly dominant eigenvalue. An eigenvalue $\lambda_1$ with $|\lambda_1| > |\lambda_j|$ for every other eigenvalue. The method needs one to exist, and the ratio $|\lambda_2/\lambda_1|$ controls the speed: small ratio, fast convergence.

The power method. Start with any $\mathbf{x}_0$ whose largest entry is $1$. Repeat:

1. Compute $A\mathbf{x}_k$.
2. Let $\mu_k$ be the entry of $A\mathbf{x}_k$ largest in absolute value.
3. Set $\mathbf{x}_{k+1} = (1/\mu_k)A\mathbf{x}_k$.

Then $\mu_k$ converges to the dominant eigenvalue and $\mathbf{x}_k$ converges to an eigenvector for it, scaled so its largest entry is $1$. The rescaling is what keeps the entries in a usable range; any other consistent normalization, such as dividing by the length, works as well.

Why it works. With $\mathbf{x}_0 = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$ and $c_1 \neq 0$,

$$ A^k\mathbf{x}_0 = \lambda_1^k\left(c_1\mathbf{v}_1 + c_2(\lambda_2/\lambda_1)^k\mathbf{v}_2 + \cdots + c_n(\lambda_n/\lambda_1)^k\mathbf{v}_n\right) $$

Every ratio in the parentheses has absolute value less than $1$, so every term after the first fades and the direction settles on $\mathbf{v}_1$.

The Rayleigh quotient. For a nonzero $\mathbf{x}$, the number

$$ R(\mathbf{x}) = \frac{\mathbf{x} \cdot A\mathbf{x}}{\mathbf{x} \cdot \mathbf{x}} $$

equals $\lambda$ exactly when $\mathbf{x}$ is an eigenvector for $\lambda$. Used on the iterates of the power method it gives an alternative estimate, and for symmetric matrices it is the estimate of choice.

The inverse power method. Pick an estimate $\alpha$ near the eigenvalue you want. The eigenvalues of $(A - \alpha I)^{-1}$ are $1/(\lambda - \alpha)$, so the eigenvalue of $A$ closest to $\alpha$ becomes the dominant eigenvalue of the inverted matrix, and the eigenvectors are unchanged. The iteration is:

1. Solve $(A - \alpha I)\mathbf{y} = \mathbf{x}_k$ for $\mathbf{y}$. Never form the inverse; factor the matrix once and reuse it.
2. Let $\mu_k$ be the entry of $\mathbf{y}$ largest in absolute value, and set $\nu_k = \alpha + 1/\mu_k$.
3. Set $\mathbf{x}_{k+1} = (1/\mu_k)\mathbf{y}$.

Then $\nu_k$ converges to the eigenvalue of $A$ nearest $\alpha$. The closer $\alpha$ starts, the faster it goes, which is why this method is usually run after a rough estimate from somewhere else.

Getting a starting estimate. Any bound on the eigenvalues works. Each eigenvalue lies within a distance of the sum of the absolute values of the off-diagonal entries of some row from that row’s diagonal entry, so the diagonal entries themselves make reasonable first guesses for $\alpha$.

When it fails. Two eigenvalues tied in absolute value leave the iterate oscillating or rotating instead of settling. A complex conjugate pair is the common case in real matrices, since the two have equal modulus. A starting vector with no component along the dominant eigenvector stalls in exact arithmetic. A ratio $|\lambda_2/\lambda_1|$ near $1$ converges so slowly that the iterate count becomes the cost.

Worked examples

The power method on a two by two

Let

$$ A = \begin{bmatrix} 6 & 5 \\ 1 & 2 \end{bmatrix}, \qquad \mathbf{x}_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix} $$

The trace is $8$ and the determinant is $12 - 5 = 7$, so the eigenvalues are the roots of $\lambda^2 - 8\lambda + 7$, namely $7$ and $1$. Knowing the answer lets you watch the method work.

Step one: $A\mathbf{x}_0 = (6,1)$. The largest entry is $\mu_0 = 6$, so $\mathbf{x}_1 = (1, \, 0.166667)$.

Step two: $A\mathbf{x}_1 = (6 + 0.833333, \; 1 + 0.333333) = (6.833333, \, 1.333333)$, so $\mu_1 = 6.833333$ and $\mathbf{x}_2 = (1, \, 0.195122)$.

Step three: $A\mathbf{x}_2 = (6.975610, \, 1.390244)$, so $\mu_2 = 6.975610$ and $\mathbf{x}_3 = (1, \, 0.199301)$.

Step four: $A\mathbf{x}_3 = (6.996503, \, 1.398601)$, so $\mu_3 = 6.996503$ and $\mathbf{x}_4 = (1, \, 0.199900)$.

The estimates $6$, $6.8333$, $6.9756$, $6.9965$ head for $7$, and the vectors head for $(1, 0.2)$. Check that limit: $A(1, 0.2) = (6 + 1, \; 1 + 0.4) = (7, 1.4) = 7(1, 0.2)$, so $(1,0.2)$ is genuinely an eigenvector for $7$, and clearing the fraction gives $(5,1)$. The error shrinks by a factor of about $|\lambda_2/\lambda_1| = 1/7$ per step, which is why four steps already give three correct digits.

The inverse power method for the other eigenvalue

Same $A$, and suppose a rough estimate says there is an eigenvalue somewhere near $0.5$. Take $\alpha = 0.5$, so

$$ A - 0.5I = \begin{bmatrix} 5.5 & 5 \\ 1 & 1.5 \end{bmatrix}, \qquad \det(A - 0.5I) = 8.25 - 5 = 3.25 $$

Start again at $\mathbf{x}_0 = (1,0)$ and solve $(A - 0.5I)\mathbf{y} = \mathbf{x}_0$. By Cramer’s rule or by hand, $\mathbf{y} = (1.5/3.25, \; -1/3.25) = (0.461538, \, -0.307692)$. The largest entry in absolute value is $\mu_0 = 0.461538$, so the eigenvalue estimate is $\nu_0 = 0.5 + 1/0.461538 = 2.666667$ and $\mathbf{x}_1 = (1, \, -0.666667)$.

Next solve $(A - 0.5I)\mathbf{y} = \mathbf{x}_1$, giving $\mathbf{y} = (1.487179, \, -1.435897)$. Then $\mu_1 = 1.487179$, $\nu_1 = 0.5 + 0.672414 = 1.172414$, and $\mathbf{x}_2 = (1, \, -0.965517)$.

Once more: $\mathbf{y} = (1.946950, \, -1.941645)$, so $\mu_2 = 1.946950$, $\nu_2 = 0.5 + 0.513624 = 1.013624$, and $\mathbf{x}_3 = (1, \, -0.997275)$.

The estimates $2.6667$, $1.1724$, $1.0136$ converge to $1$, and the vectors converge to $(1,-1)$. Check: $A(1,-1) = (6 - 5, \; 1 - 2) = (1,-1)$, so $1$ is an eigenvalue with eigenvector $(1,-1)$, exactly as the iteration claimed.

A matrix the power method cannot handle

Let

$$ A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \qquad \mathbf{x}_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix} $$

Its eigenvalues are $1$ and $-1$, equal in absolute value, so no eigenvalue is strictly dominant. Run the iteration: $A\mathbf{x}_0 = (0,1)$, so $\mu_0 = 1$ and $\mathbf{x}_1 = (0,1)$. Then $A\mathbf{x}_1 = (1,0)$, so $\mathbf{x}_2 = (1,0)$, which is where you started. The iterate cycles between two vectors forever and the direction never settles, even though every scaling factor is $1$.

The scaling factors alone would have you believe the method had converged. That is the practical warning: a stable $\mu_k$ means nothing unless $\mathbf{x}_k$ has also stopped moving, so check both.

Practice

The iteration is a chain of matrix-vector products, so keep that operation quick and exact.

Then check the iterative estimate against the eigenvalue found exactly from the characteristic polynomial, the way the worked examples do.

Quiz

Five items on matrix-vector products and eigenvalues of small matrices.