24.4 The Power Method

Chapter 19 computed a Perron vector three times and each time stopped at the same place. Section 4 bracketed $\rho(\mathbf{A})$ by the ratios of $\mathbf{A}\mathbf{x}$ to $\mathbf{x}$ and observed: “the iteration $\mathbf{x} \mapsto \mathbf{A}\mathbf{x}$ is the power method, and Chapter 24 studies how fast it converges.” Section 8 said the same of PageRank — “Chapter 24 studies it for general matrices, where the rate is governed by the ratio of the two largest eigenvalue moduli” — and Section 9’s ledger listed “the power method, the algorithm that computes a Perron vector at the size of the web” among the debts. Chapter 17 §01 and Chapter 20 §06 each promised, in the same words, that “Chapter 24 builds an iteration on” the Rayleigh quotient.

This section pays all five. It is about one idea used three times: multiply by a matrix repeatedly and the largest eigenvalue takes over. Applied to $\mathbf{A}$ itself, the idea finds the eigenvalue of largest modulus. Applied to $(\mathbf{A} - \mu\mathbf{I})^{-1}$ , it finds the eigenvalue nearest $\mu$ — any eigenvalue we choose to aim at. Applied to $(\mathbf{A} - \mu\mathbf{I})^{-1}$ with $\mu$ refreshed at every step from the Rayleigh quotient, it converges cubically, and that is the iteration Chapters 17 and 20 promised.

Throughout, $F = \mathbb{R}$ or $F = \mathbb{C}$ , matrices are in $M_n(F)$ read inside $M_n(\mathbb{C})$ , and $\left\| \cdot \right\| = \left\| \cdot \right\|_2$ . Iteration indices are superscripts in parentheses, so that subscripts stay free for coordinates.

24.4.1. Multiplying by the matrix

Take $\mathbf{A}$ diagonalizable with eigenbasis $\mathbf{v}_1, \dots, \mathbf{v}_n$ and eigenvalues $\lambda_1, \dots, \lambda_n$ , and expand a starting vector:

$\mathbf{x} = c_1\mathbf{v}_1 + \dots + c_n\mathbf{v}_n \quad\Longrightarrow\quad \mathbf{A}^k\mathbf{x} = c_1\lambda_1^k\mathbf{v}_1 + \dots + c_n\lambda_n^k\mathbf{v}_n .$

Divide by $\lambda_1^k$ , and suppose $\lvert\lambda_1\rvert$ is strictly the largest. Every coefficient except the first is multiplied by $(\lambda_i/\lambda_1)^k$ , a number of modulus less than $1$ raised to the $k$ -th power. So the sum collapses onto its first term. Nothing here needs to be discovered; it needs to be named, and its rate measured.

The power method is: multiply, rescale, repeat, and read off the dominant eigenvector.

Definition 24.4.1 (Power Iteration)

Let $\mathbf{A} \in M_n(F)$ and let $\mathbf{x}^{(0)} \in F^n$ with $\left\| \mathbf{x}^{(0)} \right\| = 1$ . The power iteration generated by $\mathbf{A}$ and $\mathbf{x}^{(0)}$ is the sequence defined for $k \ge 0$ by $\mathbf{x}^{(k+1)} = \frac{\mathbf{A}\mathbf{x}^{(k)}}{\left\| \mathbf{A}\mathbf{x}^{(k)} \right\|} ,$ which is defined as long as $\mathbf{A}\mathbf{x}^{(k)} \ne \mathbf{0}$ , together with the eigenvalue estimates $\mu_k = \bigl(\mathbf{x}^{(k)}\bigr)^{*}\mathbf{A}\,\mathbf{x}^{(k)} .$

Clause by clause. The vectors are unit vectors: the normalization has no mathematical content, since $\mathbf{x}^{(k)}$ spans the same line as $\mathbf{A}^k\mathbf{x}^{(0)}$ , and its only purpose is to stop the entries from growing or shrinking like $\lambda_1^k$ until the arithmetic can no longer hold them. The iteration is partial: if $\mathbf{A}\mathbf{x}^{(k)} = \mathbf{0}$ the sequence stops, and it stops exactly when $\mathbf{x}^{(k)} \in \operatorname{nul}(\mathbf{A})$ — which, as we shall see, a good starting vector avoids. The estimate $\mu_k$ is the natural number to attach to a unit vector $\mathbf{x}$ that is nearly an eigenvector: if $\mathbf{A}\mathbf{x} = \lambda\mathbf{x}$ exactly then $\mathbf{x}^{*}\mathbf{A}\mathbf{x} = \lambda$ . For Hermitian $\mathbf{A}$ it is the Rayleigh quotient $R_{\mathbf{A}}(\mathbf{x}^{(k)})$ of Definition 17.1.1, whose denominator is $1$ here; for a general $\mathbf{A}$ the book reserves that name for the Hermitian case, so we write the quadratic form out.

Examples.

The identity. $\mathbf{A} = \mathbf{I}$ : every unit vector is fixed, $\mathbf{x}^{(k)} = \mathbf{x}^{(0)}$ for all $k$ , and $\mu_k = 1$ . The degenerate case is worth a glance, because it shows that the method converges to an eigenvector, never to a distinguished one.
A diagonal matrix. $\mathbf{A} = \operatorname{diag}(3, 1)$ and $\mathbf{x}^{(0)} = \tfrac{1}{\sqrt2}(1,1)$ : then $\mathbf{A}^k\mathbf{x}^{(0)} \propto (3^k, 1)$ , and the normalized vector runs to $\mathbf{e}_1$ with second coordinate $3^{-k}$ . The rate $3^{-k}$ is $\lvert\lambda_2/\lambda_1\rvert^k$ .
A stochastic matrix. For the Google matrix $\mathbf{G}$ of Chapter 19 §08, $\lambda_1 = 1$ and every other eigenvalue has modulus at most $\alpha$ (Theorem 19.8.2), so the iteration $\mathbf{x}_{k+1} = \mathbf{G}\mathbf{x}_k$ of Proposition 19.8.3 is a power iteration with rate at most $\alpha$ . Chapter 19 proved the rate directly from stochasticity, with no hypothesis on $\mathbf{G}$ beyond that; the theorem below proves it from the eigenvalues, for any diagonalizable matrix.

Non-example by minimal change. Replace $\operatorname{diag}(3,1)$ by $\operatorname{diag}(3, -3)$ . The eigenvalues now have equal modulus, the collapse never happens, and from $\mathbf{x}^{(0)} = \tfrac{1}{\sqrt2}(1,1)$ the iterates alternate between $\tfrac{1}{\sqrt2}(1,1)$ and $\tfrac{1}{\sqrt2}(1,-1)$ forever. One clause of the hypothesis has been dropped — strictness in $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ — and the conclusion is gone with it.

Here is the theorem, with the rate.

Theorem 24.4.1 (The Power Method)

Let $\mathbf{A} \in M_n(F)$ be diagonalizable, with a basis $\mathbf{v}_1, \dots, \mathbf{v}_n$ of eigenvectors, $\mathbf{A}\mathbf{v}_i = \lambda_i\mathbf{v}_i$ , ordered by modulus so that $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert \ge \dots \ge \lvert\lambda_n\rvert .$ Let $\mathbf{x}^{(0)} = c_1\mathbf{v}_1 + \dots + c_n\mathbf{v}_n$ be a unit vector with $c_1 \ne 0$ , and put $\rho = \lvert\lambda_2/\lambda_1\rvert < 1$ , $\mathbf{q}_1 = \mathbf{v}_1/\left\| \mathbf{v}_1 \right\|$ and $C = \frac{1}{\left\| \mathbf{v}_1 \right\|}\sum_{i \ge 2}\Bigl\lvert\frac{c_i}{c_1}\Bigr\rvert\,\left\| \mathbf{v}_i \right\| .$ Then the power iteration of Definition 24.4.1 is defined for every $k$ , and for every $k$ with $C\rho^{k} \le \tfrac12$ there is a scalar $\omega_k$ of modulus $1$ with

(a) $\left\| \omega_k\mathbf{x}^{(k)} - \mathbf{q}_1 \right\| \le 4C\rho^{k}$ ;
(b) $\lvert\mu_k - \lambda_1\rvert \le 8\left\| \mathbf{A} \right\|\,C\rho^{k}$ .

In particular $\mu_k \to \lambda_1$ , and the iterates approach the eigenline $\operatorname{span}(\mathbf{v}_1)$ in the sense of (a), both at the rate $\rho^{k}$ .

Idea. The whole content is the display at the head of this subsection: divide $\mathbf{A}^k\mathbf{x}^{(0)}$ by $c_1\lambda_1^k$ and what is left is $\mathbf{v}_1$ plus a tail of size $\rho^k$ . Two pieces of bookkeeping then stand between that and the statement. First, normalizing: a bound on $\mathbf{y} - \mathbf{v}_1$ is not yet a bound on $\mathbf{y}/\left\| \mathbf{y} \right\| - \mathbf{v}_1/\left\| \mathbf{v}_1 \right\|$ , and the conversion costs a factor $2$ and the hypothesis $C\rho^k \le \tfrac12$ , which only says that we are past the point where the tail could swamp $\mathbf{v}_1$ . Second, phases: over $\mathbb{C}$ the scalar $c_1\lambda_1^k$ has a modulus and an argument, and normalizing kills the modulus but not the argument. That argument is the $\omega_k$ . It cancels out of $\mu_k$ , because a quadratic form does not see it.

Proof. The iterates are defined, and what they are. Since $c_1\lambda_1^k \ne 0$ , the vector $\mathbf{A}^k\mathbf{x}^{(0)} = \sum_i c_i\lambda_i^k\mathbf{v}_i$ has a non-zero coefficient on $\mathbf{v}_1$ , so $\mathbf{A}^k\mathbf{x}^{(0)} \ne \mathbf{0}$ for every $k$ , the eigenvectors being a basis. An induction on $k$ then gives $\mathbf{x}^{(k)} = \mathbf{A}^k\mathbf{x}^{(0)}/\left\| \mathbf{A}^k\mathbf{x}^{(0)} \right\|$ : true at $k = 0$ , and if it holds at $k$ then $\mathbf{A}\mathbf{x}^{(k)}$ is a positive multiple of $\mathbf{A}^{k+1}\mathbf{x}^{(0)} \ne \mathbf{0}$ , so normalizing it gives the same unit vector as normalizing $\mathbf{A}^{k+1}\mathbf{x}^{(0)}$ .

The tail. Put $\mathbf{y}_k = \frac{\mathbf{A}^k\mathbf{x}^{(0)}}{c_1\lambda_1^{k}} = \mathbf{v}_1 + \mathbf{e}_k, \qquad \mathbf{e}_k = \sum_{i \ge 2}\frac{c_i}{c_1}\Bigl(\frac{\lambda_i}{\lambda_1}\Bigr)^{k}\mathbf{v}_i .$ Since $\lvert\lambda_i/\lambda_1\rvert \le \rho$ for $i \ge 2$ , the triangle inequality (Corollary 11.1.6) gives $\left\| \mathbf{e}_k \right\| \le \rho^{k}\sum_{i\ge2}\Bigl\lvert\frac{c_i}{c_1}\Bigr\rvert\left\| \mathbf{v}_i \right\| = C\rho^{k}\left\| \mathbf{v}_1 \right\| .$

Normalizing. Let $\omega_k = \overline{\sigma_k}$ , where $\sigma_k = c_1\lambda_1^k/\lvert c_1\lambda_1^k\rvert$ has modulus $1$ . Then $\mathbf{A}^k\mathbf{x}^{(0)} = \lvert c_1\lambda_1^k\rvert\,\sigma_k\mathbf{y}_k$ , so $\mathbf{x}^{(k)} = \sigma_k\mathbf{y}_k/\left\| \mathbf{y}_k \right\|$ and $\omega_k\mathbf{x}^{(k)} = \mathbf{y}_k/\left\| \mathbf{y}_k \right\|$ . For any non-zero $\mathbf{a}, \mathbf{b}$ , $\begin{aligned} \frac{\mathbf{a}}{\left\| \mathbf{a} \right\|} - \frac{\mathbf{b}}{\left\| \mathbf{b} \right\|} &= \frac{\mathbf{a} - \mathbf{b}}{\left\| \mathbf{a} \right\|} + \mathbf{b}\Bigl(\frac{1}{\left\| \mathbf{a} \right\|} - \frac{1}{\left\| \mathbf{b} \right\|}\Bigr),\\ \Bigl\lvert\frac{1}{\left\| \mathbf{a} \right\|} - \frac{1}{\left\| \mathbf{b} \right\|}\Bigr\rvert\left\| \mathbf{b} \right\| &= \frac{\bigl\lvert\left\| \mathbf{b} \right\| - \left\| \mathbf{a} \right\|\bigr\rvert}{\left\| \mathbf{a} \right\|} \le \frac{\left\| \mathbf{a} - \mathbf{b} \right\|}{\left\| \mathbf{a} \right\|} , \end{aligned}$ the last step by the reverse triangle inequality (Lemma 16.1.1). Hence $\left\| \mathbf{a}/\left\| \mathbf{a} \right\| - \mathbf{b}/\left\| \mathbf{b} \right\| \right\| \le 2\left\| \mathbf{a}-\mathbf{b} \right\|/\left\| \mathbf{a} \right\|$ . Apply this with $\mathbf{a} = \mathbf{y}_k$ and $\mathbf{b} = \mathbf{v}_1$ . If $C\rho^k \le \tfrac12$ then $\left\| \mathbf{e}_k \right\| \le \tfrac12\left\| \mathbf{v}_1 \right\|$ , so $\left\| \mathbf{y}_k \right\| \ge \left\| \mathbf{v}_1 \right\| - \left\| \mathbf{e}_k \right\| \ge \tfrac12\left\| \mathbf{v}_1 \right\|$ , and $\left\| \omega_k\mathbf{x}^{(k)} - \mathbf{q}_1 \right\| \le \frac{2\left\| \mathbf{e}_k \right\|}{\left\| \mathbf{y}_k \right\|} \le \frac{4\left\| \mathbf{e}_k \right\|}{\left\| \mathbf{v}_1 \right\|} \le 4C\rho^{k} ,$ which is (a).

The estimates. Write $\mathbf{z} = \omega_k\mathbf{x}^{(k)}$ , a unit vector, and note $\mu_k = (\mathbf{x}^{(k)})^{*}\mathbf{A}\mathbf{x}^{(k)} = \mathbf{z}^{*}\mathbf{A}\mathbf{z}$ because $\overline{\omega_k}\omega_k = 1$ . Since $\mathbf{A}\mathbf{q}_1 = \lambda_1\mathbf{q}_1$ and $\left\| \mathbf{q}_1 \right\| = 1$ , we have $\mathbf{q}_1^{*}\mathbf{A}\mathbf{q}_1 = \lambda_1$ , so $\mu_k - \lambda_1 = \mathbf{z}^{*}\mathbf{A}(\mathbf{z} - \mathbf{q}_1) + (\mathbf{z} - \mathbf{q}_1)^{*}\mathbf{A}\mathbf{q}_1 .$ By Cauchy–Schwarz (Theorem 11.1.5) and Theorem 16.3.2 (a), each term is at most $\left\| \mathbf{A} \right\|\left\| \mathbf{z} - \mathbf{q}_1 \right\|$ , so $\lvert\mu_k - \lambda_1\rvert \le 2\left\| \mathbf{A} \right\|\left\| \mathbf{z} - \mathbf{q}_1 \right\| \le 8\left\| \mathbf{A} \right\|C\rho^k$ by (a). This proves the theorem.

Aftermath. One matrix–vector product per step: $2n^2$ flops for a dense $\mathbf{A}$ , and far fewer when $\mathbf{A}$ is sparse — for the hyperlink matrix of Chapter 19 §08, one operation per link. Nothing is factorized and nothing is stored but a vector, which is why the method survives at sizes where the matrix itself could never be written down. What it buys with that cheapness is a linear rate: to gain one decimal digit takes $1/\log_{10}(1/\rho)$ steps, and for $\rho = 0.99$ that is $230$ of them.

Example (Every iterate exactly)

Let $\mathbf{A} = \begin{psmallmatrix} 2 & 1 \\ 1 & 2\end{psmallmatrix}$ and $\mathbf{x}^{(0)} = \mathbf{e}_1$ . Compute $\mathbf{A}^k\mathbf{e}_1$ , the estimates $\mu_k$ , and the errors, in closed form.

Solution. $\mathbf{A}$ has eigenvectors $(1,1)$ and $(1,-1)$ with eigenvalues $3$ and $1$ , so $\rho = \tfrac13$ . Writing $\mathbf{e}_1 = \tfrac12(1,1) + \tfrac12(1,-1)$ , $\mathbf{A}^k\mathbf{e}_1 = \tfrac12 3^{k}(1,1) + \tfrac12(1,-1) = \Bigl(\frac{3^{k}+1}{2},\ \frac{3^{k}-1}{2}\Bigr) .$ So the iterates are proportional to $(1,0), (2,1), (5,4), (14,13), (41,40), (122,121), \dots$ , visibly running to the direction $(1,1)$ . Put $t = 3^{k}$ and $\mathbf{x} = (a,b)$ with $a = (t+1)/2$ , $b = (t-1)/2$ . Then $a^2+b^2 = (t^2+1)/2$ and $ab = (t^2-1)/4$ , so $\mu_k = \frac{2a^2 + 2ab + 2b^2}{a^2+b^2} = \frac{(t^2+1) + \tfrac12(t^2-1)}{\tfrac12(t^2+1)} = \frac{3t^2+1}{t^2+1} = 3 - \frac{2}{9^{k}+1} .$ The values are $2,\ \tfrac{14}{5},\ \tfrac{122}{41},\ \tfrac{1094}{365},\ \tfrac{9842}{3281}$ , with errors $1,\ \tfrac15,\ \tfrac1{41},\ \tfrac1{365},\ \tfrac1{3281}$ — falling by factors $5,\ 8.2,\ 8.9,\ 9.0$ , so by about a ninth in the limit, not a third.

The eigenvector error tells the other half. The component of $\mathbf{x} = (a,b)$ orthogonal to $(1,1)$ is $\tfrac12(a-b)(1,-1)$ , of length $\lvert a-b\rvert/\sqrt2 = 1/\sqrt2$ , and $\left\| \mathbf{x} \right\| = \sqrt{(t^2+1)/2}$ , so the sine of the angle between $\mathbf{x}^{(k)}$ and $(1,1)$ is $(9^{k}+1)^{-1/2}$ , which does shrink like $3^{-k} = \rho^{k}$ . Exactly: $3 - \mu_k = \frac{2}{9^k+1} = 2\sin^2\theta_k .$ The eigenvector converges at rate $\rho^{k}$ and the eigenvalue estimate at rate $\rho^{2k}$ , and the second is the square of the first. That is not an accident of this matrix.

The corollary explains it, and pays Chapter 17 §01 in its own words.

Corollary 24.4.2 (Hermitian Matrices: the Rayleigh Quotient Converges Twice as Fast)

Let $\mathbf{A} \in M_n(F)$ be Hermitian, with an orthonormal basis $\mathbf{q}_1, \dots, \mathbf{q}_n$ of eigenvectors, $\mathbf{A}\mathbf{q}_i = \lambda_i\mathbf{q}_i$ , ordered by modulus so that $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert \ge \dots \ge \lvert\lambda_n\rvert$ . Let $\mathbf{x}^{(0)} = \sum_i c_i\mathbf{q}_i$ be a unit vector with $c_1 \ne 0$ , and put $\rho = \lvert\lambda_2/\lambda_1\rvert$ . Then for every $k$ , $\bigl\lvert R_{\mathbf{A}}(\mathbf{x}^{(k)}) - \lambda_1 \bigr\rvert \ \le\ D\,\rho^{2k}, \qquad D = \frac{1}{\lvert c_1\rvert^{2}}\sum_{i \ge 2}\lvert c_i\rvert^{2}\,\lvert\lambda_i - \lambda_1\rvert .$

Idea. The Rayleigh quotient is unchanged by scaling, so evaluate it at the unnormalized $\mathbf{A}^k\mathbf{x}^{(0)}$ , where the coefficients are just $c_i\lambda_i^k$ . Both the numerator and the denominator are then sums of non-negative multiples of $\lambda_i^{2k}$ — the exponent is even, which is the whole point — so the denominator is bounded below by its first term alone, and the numerator, after $\lambda_1$ is subtracted, loses its first term entirely. What is left is $\lambda_i^{2k}$ against $\lambda_1^{2k}$ : one factor of $\rho^{2}$ per step.

Proof. Write $\mathbf{y} = \mathbf{A}^{k}\mathbf{x}^{(0)} = \sum_i c_i\lambda_i^{k}\mathbf{q}_i$ . Since $c_1\lambda_1^k \ne 0$ , $\mathbf{y} \ne \mathbf{0}$ , and $\mathbf{x}^{(k)}$ is a scalar multiple of $\mathbf{y}$ as in the proof of Theorem 24.4.1, so $R_{\mathbf{A}}(\mathbf{x}^{(k)}) = R_{\mathbf{A}}(\mathbf{y})$ by the scale-invariance of the Rayleigh quotient (Definition 17.1.1). The basis is orthonormal and $\mathbf{A}\mathbf{y} = \sum_i c_i\lambda_i^{k+1}\mathbf{q}_i$ , so by Theorem 11.2.1 (b) and (c), $R_{\mathbf{A}}(\mathbf{y}) = \frac{\langle \mathbf{A}\mathbf{y}, \mathbf{y} \rangle}{\langle \mathbf{y}, \mathbf{y} \rangle} = \frac{\sum_i \lvert c_i\rvert^{2}\lambda_i^{2k+1}}{\sum_i \lvert c_i\rvert^{2}\lambda_i^{2k}} ,$ all the $\lambda_i$ being real (Proposition 11.6.7). Subtract $\lambda_1$ , which is the same as subtracting $\lambda_1$ times the denominator from the numerator; the term $i = 1$ cancels, leaving $R_{\mathbf{A}}(\mathbf{y}) - \lambda_1 = \frac{\sum_{i \ge 2}\lvert c_i\rvert^{2}\lambda_i^{2k}(\lambda_i - \lambda_1)}{\sum_i \lvert c_i\rvert^{2}\lambda_i^{2k}} .$ Every $\lambda_i^{2k}$ is non-negative, so the denominator is at least its first term $\lvert c_1\rvert^2\lambda_1^{2k} = \lvert c_1\rvert^2\lvert\lambda_1\rvert^{2k}$ , and the numerator is at most $\lvert\lambda_2\rvert^{2k}\sum_{i\ge2}\lvert c_i\rvert^2\lvert\lambda_i - \lambda_1\rvert$ in modulus, since $\lvert\lambda_i\rvert \le \lvert\lambda_2\rvert$ for $i \ge 2$ . Dividing gives the stated bound with $\rho^{2k} = (\lvert\lambda_2\rvert/\lvert\lambda_1\rvert)^{2k}$ . This proves the corollary.

This is the “second-order accuracy” of Chapter 17 §01, seen along an iteration: the eigenvector error is $\rho^k$ and the eigenvalue error is its square. Theorem 20.6.2 says the same thing structurally and without any iteration — a residual of size $\varepsilon$ yields an eigenvalue error of at most $\varepsilon^2$ divided by the gap — and the corollary above is that estimate made explicit for the particular unit vectors the power method produces.

Quick check. In Example (Every iterate exactly) the estimates $\mu_k$ increase towards $3$ and never overshoot. Which earlier result says that they cannot?

Solution. Proposition 17.1.1: for a Hermitian $\mathbf{A}$ , the values of $R_{\mathbf{A}}$ fill exactly the interval $[\lambda_n(\mathbf{A}), \lambda_1(\mathbf{A})]$ between the smallest and largest eigenvalues. Here $\mathbf{A} = \begin{psmallmatrix}2&1\\mathbf{1}&2\end{psmallmatrix}$ has eigenvalues $3$ and $1$ , so every $\mu_k = R_{\mathbf{A}}(\mathbf{x}^{(k)})$ lies in $[1, 3]$ and no estimate can exceed $3$ . The power method on a Hermitian matrix with $\lambda_1 > 0$ dominant therefore produces a sequence of certified lower bounds for $\lambda_1$ , which is what Chapter 19 §04’s brackets were doing by hand.

Example (The method that never converges)

Let $\mathbf{P} \in M_3(\mathbb{R})$ be the cyclic permutation matrix of Example (The Cyclic Permutation), with columns $\mathbf{e}_2, \mathbf{e}_3, \mathbf{e}_1$ , so that $\mathbf{P}(x_1,x_2,x_3) = (x_3, x_1, x_2)$ . Run the power iteration from $\mathbf{x}^{(0)} = \mathbf{e}_1$ and from $\mathbf{x}^{(0)} = \tfrac{1}{\sqrt{14}}(1,2,3)$ .

Solution. By Example (The Cyclic Permutation) the eigenvalues of $\mathbf{P}$ are the three cube roots of unity, all of modulus $1$ : the hypothesis $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ fails as badly as it can. Every $\mathbf{P}^k$ is a permutation matrix, so no normalization is needed, and $\mathbf{P}^3 = \mathbf{I}$ . From $\mathbf{e}_1$ the iterates are $\mathbf{e}_1 \to \mathbf{e}_2 \to \mathbf{e}_3 \to \mathbf{e}_1 \to \dots$ and from $\tfrac{1}{\sqrt{14}}(1,2,3)$ they are $\tfrac{1}{\sqrt{14}}$ times $(1,2,3), (3,1,2), (2,3,1), (1,2,3), \dots$ . In both cases the sequence takes three distinct values over and over, so neither it nor the line it spans has a limit. The estimates behave worse than that: they are constant, $\mu_k = 0$ for every $k$ in the first case and $\mu_k = \tfrac{11}{14}$ in the second, so they converge perfectly well — to a number that is not an eigenvalue of $\mathbf{P}$ at all. Yet $\mathbf{P}$ does have the eigenvalue $1$ , with eigenvector $\mathbf{1}$ , and the method never finds it.

The failure is exactly the one Chapter 19 §04 warned of for the swap $\begin{psmallmatrix}0&1\\mathbf{1}&0\end{psmallmatrix}$ , whose brackets “never improve”. It is not a defect of the arithmetic and no amount of care repairs it: the hypothesis of Theorem 24.4.1 is false, and its conclusion is false too.

Warning. Two hypotheses, two different kinds of failure. If $\lvert\lambda_1\rvert = \lvert\lambda_2\rvert$ , as in Example (The method that never converges), the method fails permanently and no starting vector helps. If instead $c_1 = 0$ — the starting vector happens to miss the dominant eigendirection — the method converges, but to the wrong eigenvector. The second failure is fragile and the first is not: a starting vector with $c_1 = 0$ lies in a proper subspace, and in the arithmetic of Definition 24.1.1 the first rounding error generally puts a small non-zero $c_1$ back, after which Theorem 24.4.1 takes over, only $C$ being large. This is the one situation in the chapter where rounding helps, and it is not a reason to rely on it: a starting vector chosen without regard to $\mathbf{A}$ is a better plan than a lucky rounding.

24.4.2. Shifting and inverting

The power method finds the eigenvalue of largest modulus, and that is usually not the one anybody wants. The repair is a change of matrix. If $\mu$ is not an eigenvalue of $\mathbf{A}$ , then $\mathbf{A} - \mu\mathbf{I}$ is invertible, and the eigenvalues of its inverse are the numbers $(\lambda_i - \mu)^{-1}$ — large exactly for the $\lambda_i$ close to $\mu$ . Multiplying repeatedly by $(\mathbf{A} - \mu\mathbf{I})^{-1}$ therefore singles out the eigenvalue nearest to $\mu$ , and $\mu$ is ours to choose.

Definition 24.4.2 (Inverse Iteration)

Let $\mathbf{A} \in M_n(F)$ , let $\mu \in F$ be not an eigenvalue of $\mathbf{A}$ , and let $\mathbf{x}^{(0)}$ be a unit vector. Inverse iteration with shift $\mu$ is the power iteration of Definition 24.4.1 generated by $(\mathbf{A} - \mu\mathbf{I})^{-1}$ and $\mathbf{x}^{(0)}$ ; that is, $(\mathbf{A} - \mu\mathbf{I})\,\mathbf{z}^{(k)} = \mathbf{x}^{(k)}, \qquad \mathbf{x}^{(k+1)} = \frac{\mathbf{z}^{(k)}}{\left\| \mathbf{z}^{(k)} \right\|} .$ The eigenvalue estimates are taken from $\mathbf{A}$ itself: $\mu_k = (\mathbf{x}^{(k)})^{*}\mathbf{A}\mathbf{x}^{(k)}$ .

The shifted matrix is never inverted. One factorization of $\mathbf{A} - \mu\mathbf{I}$ by elimination with pivoting is computed once, at a cost of $\tfrac23n^3$ flops (Section 2), and each step then costs two triangular solves, $2n^2$ flops — the same price per step as the power method, after a one-off payment.

Theorem 24.4.3 (Inverse Iteration Converges to the Nearest Eigenvector)

Let $\mathbf{A} \in M_n(F)$ be diagonalizable with eigenbasis $\mathbf{v}_1, \dots, \mathbf{v}_n$ and eigenvalues $\lambda_1, \dots, \lambda_n$ , and let $\mu \in F$ be not an eigenvalue of $\mathbf{A}$ . Suppose there is an index $j$ with $\lvert\lambda_j - \mu\rvert < \lvert\lambda_i - \mu\rvert \qquad \text{for every } i \ne j ,$ so that $\lambda_j$ is the unique eigenvalue nearest $\mu$ . Let $\mathbf{x}^{(0)}$ be a unit vector whose $\mathbf{v}_j$ -coefficient is non-zero. Then inverse iteration with shift $\mu$ is defined for every $k$ , the iterates approach the eigenline $\operatorname{span}(\mathbf{v}_j)$ in the sense of Theorem 24.4.1 (a), and $\mu_k \to \lambda_j$ , both at the rate $\rho_{\mu}^{k}$ , where $\rho_{\mu} = \frac{\lvert\lambda_j - \mu\rvert}{\min_{i \ne j}\lvert\lambda_i - \mu\rvert} < 1 .$

Idea. There is nothing to prove that Theorem 24.4.1 has not proved. The only work is to check that $(\mathbf{A} - \mu\mathbf{I})^{-1}$ satisfies its hypotheses — same eigenvectors, reciprocal shifted eigenvalues, a strict winner in modulus — and to read off what the ratio $\rho$ of that theorem becomes after the reciprocals are taken: dividing by the largest reciprocal is multiplying by the smallest difference, so the ratio turns upside down.

Proof. Let $\mathbf{V} = (\mathbf{v}_1 \mid \dots \mid \mathbf{v}_n)$ , invertible because its columns are a basis, and $\boldsymbol{\Lambda} = \operatorname{diag}(\lambda_1, \dots, \lambda_n)$ , so that $\mathbf{A} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}$ (Theorem 9.4.1). Then $\mathbf{A} - \mu\mathbf{I} = \mathbf{V}(\boldsymbol{\Lambda} - \mu\mathbf{I})\mathbf{V}^{-1} ,$ and $\boldsymbol{\Lambda} - \mu\mathbf{I}$ is invertible because no $\lambda_i$ equals $\mu$ . Hence $\mathbf{B} \mathrel{:=} (\mathbf{A} - \mu\mathbf{I})^{-1} = \mathbf{V}(\boldsymbol{\Lambda} - \mu\mathbf{I})^{-1}\mathbf{V}^{-1}$ is diagonalizable, with the same eigenvectors $\mathbf{v}_i$ and eigenvalues $\beta_i = (\lambda_i - \mu)^{-1}$ .

Their moduli are $\lvert\beta_i\rvert = 1/\lvert\lambda_i - \mu\rvert$ , so the hypothesis on $j$ says exactly that $\lvert\beta_j\rvert > \lvert\beta_i\rvert$ for every $i \ne j$ : the eigenvalue $\beta_j$ is strictly dominant. Relabeling so that $j$ comes first, the second largest modulus is $\max_{i\ne j}\lvert\beta_i\rvert = 1/\min_{i\ne j}\lvert\lambda_i-\mu\rvert$ , so the ratio of Theorem 24.4.1 is $\frac{\max_{i\ne j}\lvert\beta_i\rvert}{\lvert\beta_j\rvert} = \frac{\lvert\lambda_j-\mu\rvert}{\min_{i\ne j}\lvert\lambda_i-\mu\rvert} = \rho_{\mu} .$ Since the $\mathbf{v}_j$ -coefficient of $\mathbf{x}^{(0)}$ is non-zero, Theorem 24.4.1 applies to $\mathbf{B}$ and $\mathbf{x}^{(0)}$ and gives unimodular $\omega_k$ with $\left\| \omega_k\mathbf{x}^{(k)} - \mathbf{q}_j \right\| \le 4C\rho_{\mu}^{k}$ , where $\mathbf{q}_j = \mathbf{v}_j/\left\| \mathbf{v}_j \right\|$ and $C$ is that theorem’s constant for $\mathbf{B}$ and $\mathbf{x}^{(0)}$ . That is the convergence of the line.

For the estimates, repeat the last paragraph of the proof of Theorem 24.4.1 with $\mathbf{A}$ in place of $\mathbf{B}$ : with $\mathbf{z} = \omega_k\mathbf{x}^{(k)}$ we have $\mu_k = \mathbf{z}^{*}\mathbf{A}\mathbf{z}$ and $\mathbf{q}_j^{*}\mathbf{A}\mathbf{q}_j = \lambda_j$ , so $\lvert\mu_k - \lambda_j\rvert \le 2\left\| \mathbf{A} \right\|\left\| \mathbf{z} - \mathbf{q}_j \right\| \le 8\left\| \mathbf{A} \right\|C\rho_{\mu}^{k}$ . This proves the theorem.

Example (Aiming at the smaller eigenvalue)

Let $\mathbf{A} = \begin{psmallmatrix}2&1\\mathbf{1}&2\end{psmallmatrix}$ again, and take the shift $\mu = \tfrac34$ , which is nearer to the eigenvalue $1$ than to $3$ . Run inverse iteration from $\mathbf{x}^{(0)} = \mathbf{e}_1$ , and compare with Example (Every iterate exactly).

Solution. $\mathbf{A} - \tfrac34\mathbf{I} = \begin{psmallmatrix} 5/4 & 1 \\ 1 & 5/4\end{psmallmatrix}$ has determinant $\tfrac{25}{16} - 1 = \tfrac9{16}$ , so $\bigl(\mathbf{A} - \tfrac34\mathbf{I}\bigr)^{-1} = \tfrac{16}{9}\begin{pmatrix} 5/4 & -1 \\ -1 & 5/4\end{pmatrix} = \tfrac19\begin{pmatrix} 20 & -16 \\ -16 & 20\end{pmatrix} .$ Ignoring positive scalars, each step is $(a, b) \mapsto (5a - 4b,\ -4a + 5b)$ , so the iterates are $(1, 0) \to (5, -4) \to (41, -40) \to (365, -364) \to (3281, -3280) \to \dots ,$ running to the direction $(1,-1)$ — the eigenvector of the eigenvalue $1$ . In closed form: the eigenvalues of $(\mathbf{A} - \tfrac34\mathbf{I})^{-1}$ are $(3 - \tfrac34)^{-1} = \tfrac49$ and $(1 - \tfrac34)^{-1} = 4$ , so expanding $\mathbf{e}_1 = \tfrac12(1,1) + \tfrac12(1,-1)$ gives $\bigl(\mathbf{A} - \tfrac34\mathbf{I}\bigr)^{-k}\mathbf{e}_1 = \tfrac12\bigl(\tfrac49\bigr)^{k}(1,1) + \tfrac12 4^{k}(1,-1) \ \propto\ \bigl(9^{k}+1,\ -(9^{k}-1)\bigr) ,$ the last step on multiplying by $2(9/4)^{k}$ . Halving, this is the list above. Writing $\mathbf{x} = (a, b)$ with $a + b = 1$ , we have $2ab = 1 - (a^2+b^2)$ and so $\mu_k = \frac{2a^2 + 2ab + 2b^2}{a^2+b^2} = 1 + \frac{1}{a^2+b^2} = 1 + \frac{2}{81^{k}+1} ,$ since $a^2 + b^2 = \tfrac12(81^{k}+1)$ . The errors are $1, \tfrac1{41}, \tfrac1{3281}, \tfrac1{265721}, \tfrac1{21523361}, \dots$ , each about $\tfrac1{81}$ of the last, while the sine of the angle to $(1,-1)$ is $(81^{k}+1)^{-1/2}$ , shrinking by $\tfrac19$ per step. The theorem predicts exactly that: $\rho_{\mu} = \frac{\lvert 1 - \tfrac34\rvert}{\lvert 3 - \tfrac34\rvert} = \frac{1/4}{9/4} = \frac19 ,$ with the estimates again going at the squared rate $\rho_{\mu}^{2k} = 81^{-k}$ , as in Example (Every iterate exactly) and for the same reason. Compare that example, where the same matrix and the same starting vector gave rate $\tfrac13$ towards the other eigenvector: one well-chosen shift has both redirected the method and made it gain twice as many digits per step.

The objection, and the answer. If $\mu$ is close to $\lambda_j$ , then $\mathbf{A} - \mu\mathbf{I}$ is close to singular and $\kappa_2(\mathbf{A} - \mu\mathbf{I})$ is enormous. For normal $\mathbf{A} = \mathbf{U}\mathbf{D}\mathbf{U}^{*}$ (Corollary 12.4.3) we have $\mathbf{A} - \mu\mathbf{I} = \mathbf{U}(\mathbf{D} - \mu\mathbf{I})\mathbf{U}^{*}$ , whose singular values are the numbers $\lvert\lambda_i - \mu\rvert$ , so by Proposition 16.8.2 (d) $\kappa_2(\mathbf{A} - \mu\mathbf{I}) = \frac{\max_i\lvert\lambda_i - \mu\rvert}{\lvert\lambda_j - \mu\rvert} ,$ which grows without bound as $\mu \to \lambda_j$ . Section 3 spent its length on how much a large condition number costs. And the better the shift, the worse the conditioning: the method is deliberately solving an ill-conditioned system, and it does so on purpose.

It does no harm, and here is why. Suppose the linear solve is backward stable in the sense of Definition 24.1.3, so that the computed $\widehat{\mathbf{z}}$ is the exact solution of $(\mathbf{A} - \mu\mathbf{I} + \boldsymbol{\Delta})\widehat{\mathbf{z}} = \mathbf{x}, \qquad \left\| \boldsymbol{\Delta} \right\| \le c\,u\,\left\| \mathbf{A} - \mu\mathbf{I} \right\| ,$ for a unit $\mathbf{x}$ . The computed $\widehat{\mathbf{z}}$ may be wrong in every digit — the forward error is bounded only by $\kappa_2(\mathbf{A}-\mu\mathbf{I})cu$ , which may exceed $1$ . But rearrange the same equation as a residual: $(\mathbf{A} - \mu\mathbf{I})\widehat{\mathbf{z}} - \mathbf{x} = -\boldsymbol{\Delta}\widehat{\mathbf{z}}$ , so the normalized vector $\widehat{\mathbf{q}} = \widehat{\mathbf{z}}/\left\| \widehat{\mathbf{z}} \right\|$ satisfies $\left\| \mathbf{A}\widehat{\mathbf{q}} - \mu\widehat{\mathbf{q}} \right\| \ \le\ \frac{1}{\left\| \widehat{\mathbf{z}} \right\|} + c\,u\,\left\| \mathbf{A} - \mu\mathbf{I} \right\| .$ The first term is small when $\left\| \widehat{\mathbf{z}} \right\|$ is large, and that takes an argument and a hypothesis. Near-singularity alone will not supply it: for $\mathbf{A} = \operatorname{diag}(1,2)$ , $\mu = 1 + 10^{-10}$ and $\mathbf{x} = \mathbf{e}_2$ , the matrix $\mathbf{A} - \mu\mathbf{I}$ is as near singular as one pleases and yet $\mathbf{z} = (\mathbf{A} - \mu\mathbf{I})^{-1}\mathbf{x}$ has norm $1$ to ten digits. What is needed is that $\mathbf{x}$ see the near-singular direction.

Take $\mathbf{A}$ normal with orthonormal eigenbasis $\mathbf{q}_1, \dots, \mathbf{q}_n$ (Corollary 12.4.3), so that $\mathbf{A}^{*}\mathbf{q}_j = \overline{\lambda_j}\mathbf{q}_j$ and hence $\mathbf{q}_j^{*}(\mathbf{A} - \mu\mathbf{I}) = (\lambda_j - \mu)\mathbf{q}_j^{*}$ . Multiplying the perturbed equation by $\mathbf{q}_j^{*}$ gives $\mathbf{q}_j^{*}\mathbf{x} = (\lambda_j - \mu)\mathbf{q}_j^{*}\widehat{\mathbf{z}} + \mathbf{q}_j^{*}\boldsymbol{\Delta}\widehat{\mathbf{z}}$ , so by Cauchy–Schwarz (Theorem 11.1.5) and Theorem 16.3.2 (a),

$\left\| \widehat{\mathbf{z}} \right\| \ \ge\ \frac{\lvert\mathbf{q}_j^{*}\mathbf{x}\rvert}{\lvert\lambda_j - \mu\rvert + \left\| \boldsymbol{\Delta} \right\|} .$

So a shift close to $\lambda_j$ does force $\left\| \widehat{\mathbf{z}} \right\|$ to be large, provided the component $\mathbf{q}_j^{*}\mathbf{x}$ is not small — which is the hypothesis Theorem 24.4.3 already makes on the starting vector. Granting it, the first term of the residual bound is small, and the second is at the rounding level. The residual of the normalized answer is small even though the answer itself is inaccurate, and a small residual is all that an eigenvector approximation has to have.

For Hermitian $\mathbf{A}$ one can say precisely how much it buys. Let $\widehat{\mathbf{q}}$ be a unit vector, $\widehat\mu = R_{\mathbf{A}}(\widehat{\mathbf{q}})$ its Rayleigh quotient and $\varepsilon = \left\| \mathbf{A}\widehat{\mathbf{q}} - \widehat\mu\widehat{\mathbf{q}} \right\|$ its residual. Expand $\widehat{\mathbf{q}} = \sum_i d_i\mathbf{q}_i$ in an orthonormal eigenbasis (Corollary 12.4.3); by Theorem 11.2.1 (c), $\varepsilon^{2} = \sum_{i}\lvert d_i\rvert^{2}(\lambda_i - \widehat\mu)^{2} .$ If $\lambda_j$ is the only eigenvalue within distance $g$ of $\widehat\mu$ , every term with $i \ne j$ has $(\lambda_i - \widehat\mu)^2 \ge g^2$ , so $\sum_{i\ne j}\lvert d_i\rvert^2 \le \varepsilon^2/g^2$ and the sine of the angle between $\widehat{\mathbf{q}}$ and $\mathbf{q}_j$ is at most $\varepsilon/g$ . Together with Theorem 20.6.1, which puts an eigenvalue of $\mathbf{A}$ within $\varepsilon$ of $\widehat\mu$ , and Theorem 20.6.2, which improves that to $\varepsilon^2/g$ , the small residual certifies both halves of the answer.

Warning. “Ill-conditioned” is a statement about the wrong question. The system $(\mathbf{A} - \mu\mathbf{I})\mathbf{z} = \mathbf{x}$ is indeed almost unsolvable if what is wanted is $\mathbf{z}$ . What is wanted is the direction of $\mathbf{z}$ , and the direction is precisely the part the ill-conditioning leaves alone: the huge factor $1/(\lambda_j - \mu)$ multiplies the component along $\mathbf{v}_j$ , which is the component being sought, and dividing by $\left\| \widehat{\mathbf{z}} \right\|$ removes it. The residual estimate above is the rigorous form of that sentence. The same reasoning does not license inverse iteration when the eigenvalue is not simple or two eigenvalues are nearly equal: then $g$ is small, the bound $\varepsilon/g$ says nothing, and the computed direction really is unreliable.

24.4.3. The shift that improves itself

Inverse iteration needs a shift, and a better shift converges faster. But the iteration is manufacturing better and better approximations to $\lambda_j$ as it goes, in the estimates $\mu_k$ . Use them.

Definition 24.4.3 (Rayleigh Quotient Iteration)

Let $\mathbf{A} \in M_n(F)$ be Hermitian and let $\mathbf{x}^{(0)}$ be a unit vector. Rayleigh quotient iteration is the sequence defined for $k \ge 0$ by $\mu_k = R_{\mathbf{A}}(\mathbf{x}^{(k)}) = (\mathbf{x}^{(k)})^{*}\mathbf{A}\mathbf{x}^{(k)}, \qquad (\mathbf{A} - \mu_k\mathbf{I})\mathbf{z}^{(k)} = \mathbf{x}^{(k)}, \qquad \mathbf{x}^{(k+1)} = \frac{\mathbf{z}^{(k)}}{\left\| \mathbf{z}^{(k)} \right\|} ,$ the step being undefined, and the iteration stopped, when $\mu_k$ is an eigenvalue of $\mathbf{A}$ — in which case $\mu_k$ is the answer, exactly.

It is inverse iteration with the shift refreshed at every step, and it is what Chapter 17 §01 and Chapter 20 §06 both promised. The price is a new factorization of $\mathbf{A} - \mu_k\mathbf{I}$ at every step, $\tfrac23n^3$ flops instead of $2n^2$ . What is bought is the following rate.

Theorem 24.4.4 (Rayleigh Quotient Iteration Converges Cubically)

Let $\mathbf{A} \in M_n(F)$ be Hermitian with an orthonormal basis $\mathbf{q}_1, \dots, \mathbf{q}_n$ of eigenvectors, $\mathbf{A}\mathbf{q}_i = \lambda_i\mathbf{q}_i$ , and fix an index $j$ whose eigenvalue is simple. Put $g = \min_{i \ne j}\lvert\lambda_i - \lambda_j\rvert > 0, \qquad D = \max_{i \ne j}\lvert\lambda_i - \lambda_j\rvert ,$ and for a unit vector $\mathbf{x}$ write $s(\mathbf{x}) = \left\| \mathbf{x} - \mathbf{q}_j\mathbf{q}_j^{*}\mathbf{x} \right\|$ , the length of the part of $\mathbf{x}$ orthogonal to $\mathbf{q}_j$ .

Let $\mathbf{x} = \mathbf{x}^{(k)}$ be a unit iterate with $s = s(\mathbf{x})$ satisfying $s \le \tfrac12$ and $Ds^{2} \le \tfrac{g}{2}$ . Then $\lvert\mu_k - \lambda_j\rvert \le Ds^{2}$ , and either $\mu_k = \lambda_j$ — the iteration has found the eigenvalue exactly — or the next iterate is defined and $s\bigl(\mathbf{x}^{(k+1)}\bigr) \ \le\ \frac{4D}{g}\,s^{3} .$ Consequently, if $s_0 = s(\mathbf{x}^{(0)})$ satisfies $s_0 \le \tfrac12$ and $s_0^{2} \le g/(8D)$ , then $a_k \mathrel{:=} \sqrt{4D/g}\;s(\mathbf{x}^{(k)})$ satisfies $a_k \le a_0^{3^{k}}$ with $a_0 \le 1/\sqrt2$ , so $s(\mathbf{x}^{(k)}) \to 0$ and $\mu_k \to \lambda_j$ .

Idea. Two facts, multiplied together. The Rayleigh quotient is accurate to second order in the eigenvector error, so a shift built from an iterate at distance $s$ is wrong by only $Ds^2$ — that is Corollary 24.4.2’s squaring again, in its non-iterative form. And one step of inverse iteration with a shift wrong by $\delta$ multiplies the error by about $\delta/g$ . Put $\delta = Ds^2$ and the step multiplies $s$ by $Ds^2/g$ : a factor $s^2$ per step, which is one cube. Everything else is the arithmetic of keeping the shift far enough from the other eigenvalues for $g$ to survive, and that is what the two smallness hypotheses are for.

Proof. Expand $\mathbf{x} = \sum_i c_i\mathbf{q}_i$ with $\sum_i\lvert c_i\rvert^2 = 1$ (Theorem 11.2.1 (c)). Then $\mathbf{q}_j\mathbf{q}_j^{*}\mathbf{x} = c_j\mathbf{q}_j$ , so $s^{2} = \sum_{i \ne j}\lvert c_i\rvert^{2}, \qquad \lvert c_j\rvert^{2} = 1 - s^{2} .$

Step 1: the shift is accurate to second order. By Theorem 11.2.1 (b), $\mu_k = \mathbf{x}^{*}\mathbf{A}\mathbf{x} = \sum_i\lvert c_i\rvert^2\lambda_i$ , and since the weights sum to $1$ , $\mu_k - \lambda_j = \sum_{i \ne j}\lvert c_i\rvert^{2}(\lambda_i - \lambda_j), \qquad \lvert\mu_k - \lambda_j\rvert \le D\sum_{i\ne j}\lvert c_i\rvert^2 = Ds^{2} .$

Step 2: the shift stays away from the other eigenvalues. By hypothesis $Ds^2 \le g/2$ , so $\lvert\mu_k - \lambda_j\rvert \le g/2$ and, for $i \ne j$ , $\lvert\lambda_i - \mu_k\rvert \ge \lvert\lambda_i - \lambda_j\rvert - \lvert\lambda_j - \mu_k\rvert \ge g - \tfrac g2 = \tfrac g2 > 0 .$ In particular the only eigenvalue $\mu_k$ can equal is $\lambda_j$ . If it does, the iteration stops with the exact answer. Otherwise $\mathbf{A} - \mu_k\mathbf{I}$ is invertible and the step is defined.

Step 3: one step. Since $\mathbf{A}\mathbf{q}_i = \lambda_i\mathbf{q}_i$ , $\mathbf{z}^{(k)} = (\mathbf{A} - \mu_k\mathbf{I})^{-1}\mathbf{x} = \sum_i \frac{c_i}{\lambda_i - \mu_k}\,\mathbf{q}_i .$ The quantity $s(\cdot)$ is unchanged by scaling, so $s(\mathbf{x}^{(k+1)}) = s(\mathbf{z}^{(k)}/\left\| \mathbf{z}^{(k)} \right\|)$ is the ratio of the length of the part of $\mathbf{z}^{(k)}$ orthogonal to $\mathbf{q}_j$ to the length of $\mathbf{z}^{(k)}$ : $s\bigl(\mathbf{x}^{(k+1)}\bigr)^{2} = \frac{\displaystyle\sum_{i \ne j}\frac{\lvert c_i\rvert^{2}}{\lvert\lambda_i - \mu_k\rvert^{2}}} {\displaystyle\sum_{i}\frac{\lvert c_i\rvert^{2}}{\lvert\lambda_i - \mu_k\rvert^{2}}} .$ Bound the numerator above using Step 2, and the denominator below by its $i = j$ term alone: $s\bigl(\mathbf{x}^{(k+1)}\bigr)^{2} \ \le\ \frac{(2/g)^{2}\,s^{2}}{\lvert c_j\rvert^{2}/\lvert\lambda_j - \mu_k\rvert^{2}} \ =\ \frac{4s^{2}\lvert\lambda_j - \mu_k\rvert^{2}}{g^{2}\,(1 - s^{2})} \ \le\ \frac{4s^{2}D^{2}s^{4}}{g^{2}(1-s^{2})} ,$ using Step 1 at the end. Taking square roots and using $s \le \tfrac12$ , so that $\sqrt{1-s^2} \ge \tfrac{\sqrt3}{2} \ge \tfrac12$ , $s\bigl(\mathbf{x}^{(k+1)}\bigr) \ \le\ \frac{2D\,s^{3}}{g\sqrt{1-s^{2}}} \ \le\ \frac{4D}{g}\,s^{3} .$

Step 4: the iteration. Put $C = 4D/g$ and $a_k = \sqrt{C}\,s(\mathbf{x}^{(k)})$ . The hypotheses $s_0 \le \tfrac12$ and $s_0^2 \le g/(8D)$ give $a_0^2 = Cs_0^2 \le \tfrac12$ , so $a_0 \le 1/\sqrt2 < 1$ , and also $Ds_0^2 \le g/2$ . Whenever the hypotheses hold at step $k$ , Step 3 gives $a_{k+1} \le \sqrt C\,C s_k^3 = a_k^{3}$ , so $a_{k+1} \le a_k^3 \le a_k \le a_0$ ; hence $s_{k+1} \le s_k \le s_0 \le \tfrac12$ and $Ds_{k+1}^2 \le Ds_0^2 \le g/2$ , and the hypotheses hold again. By induction $a_k \le a_0^{3^k}$ for every $k$ . Since $a_0 < 1$ , $a_k \to 0$ and so $s(\mathbf{x}^{(k)}) \to 0$ ; by Step 1, $\lvert\mu_k - \lambda_j\rvert \le Ds_k^2 \to 0$ . This proves the theorem.

Cubic convergence is a different animal from the linear convergence of Theorem 24.4.1. Each step roughly triples the number of correct digits, so an iterate with two correct digits has six after one step and eighteen after two — which is more than the arithmetic holds. In practice the iteration is over in three or four steps from any reasonable start, and the cost is dominated by those three or four factorizations.

Example (Three steps, exactly)

Run Rayleigh quotient iteration on $\mathbf{A} = \begin{psmallmatrix}2&1\\mathbf{1}&2\end{psmallmatrix}$ from $\mathbf{x}^{(0)} \propto (2,1)$ , in exact arithmetic.

Solution. Write each iterate as an integer vector, scaling being irrelevant. Take $\mathbf{q}_j = \tfrac{1}{\sqrt2}(1,1)$ , so $\lambda_j = 3$ , and here $g = D = 2$ , giving $4D/g = 4$ .

Step 0. $\mathbf{x} \propto (2,1)$ . Then $\mu_0 = (2\cdot 2^2 + 2\cdot2\cdot1 + 2\cdot 1^2)/5 = \tfrac{14}{5}$ , and $\mathbf{A} - \tfrac{14}5\mathbf{I} = \begin{psmallmatrix} -4/5 & 1 \\ 1 & -4/5\end{psmallmatrix}$ has determinant $-\tfrac9{25}$ . Solving gives $\mathbf{z} \propto (13, 14)$ .

Step 1. $\mathbf{x} \propto (13,14)$ , $\mu_1 = \tfrac{1094}{365}$ , and $3 - \mu_1 = \tfrac1{365}$ . Solving again gives $\mathbf{z} \propto (9842, 9841)$ .

Step 2. $\mu_2 = \tfrac{581130734}{193710245}$ , with $3 - \mu_2 = \tfrac{1}{193710245}$ .

Step 3. $\mathbf{x} \propto (3812798742493,\ 3812798742494)$ and $3 - \mu_3 = \tfrac1{29074868501520029845195085}$ .

The eigenvalue errors are $2\cdot10^{-1},\quad 2.7\cdot10^{-3},\quad 5.2\cdot10^{-9},\quad 3.4\cdot10^{-26} ,$ and the eigenvector errors $s_k$ are $3.16\cdot10^{-1},\quad 3.70\cdot10^{-2},\quad 5.08\cdot10^{-5},\quad 1.31\cdot10^{-13} .$ Check the theorem: $s_1/s_0^{3} = 1.17$ , $s_2/s_1^{3} = 1.00$ , $s_3/s_2^{3} = 1.00$ , all comfortably below the predicted constant $4D/g = 4$ . Four steps have produced twenty-five correct digits from a starting vector accurate to none.

Remark. The non-Hermitian rate is quadratic, and we do not prove it. For a general diagonalizable $\mathbf{A}$ , with $\mu_k$ taken to be $(\mathbf{x}^{(k)})^{*}\mathbf{A}\mathbf{x}^{(k)}$ , the same iteration converges locally at a quadratic rate rather than a cubic one. The reason is visible in Step 1 above: it used that $\mu_k$ is a weighted average of the eigenvalues with weights $\lvert c_i\rvert^2$ , so that the first-order term in $s$ is absent — and that identity is exactly the Hermitian spectral theorem at work. Without it the shift is accurate only to first order, $\lvert\mu_k - \lambda_j\rvert = O(s)$ , and the step multiplies $s$ by $O(s)$ . The general statement is classical and is not proved here, and nothing in this book depends on it.

24.4.4. Deflation, and what comes next

One eigenpair is rarely enough. For Hermitian $\mathbf{A}$ , once a unit eigenvector $\mathbf{q}_1$ and its eigenvalue $\lambda_1$ are known, the next one can be reached by removing the first: put

$\mathbf{A}' = \mathbf{A} - \lambda_1\mathbf{q}_1\mathbf{q}_1^{*} .$

Let $\mathbf{q}_1, \dots, \mathbf{q}_n$ be an orthonormal eigenbasis (Corollary 12.4.3). Then $\mathbf{A}'\mathbf{q}_1 = \lambda_1\mathbf{q}_1 - \lambda_1\mathbf{q}_1(\mathbf{q}_1^{*}\mathbf{q}_1) = \mathbf{0}$ , and for $i \ne 1$ , $\mathbf{A}'\mathbf{q}_i = \lambda_i\mathbf{q}_i - \lambda_1\mathbf{q}_1(\mathbf{q}_1^{*}\mathbf{q}_i) = \lambda_i\mathbf{q}_i$ by orthonormality. So $\mathbf{A}'$ is Hermitian with the same eigenvectors and the eigenvalue list $0, \lambda_2, \dots, \lambda_n$ : the power method applied to $\mathbf{A}'$ finds $\lambda_2$ , provided $\lambda_2$ is now dominant in modulus.

Warning. Deflation is exact only with an exact eigenvector. If the computed $\widehat\mathbf{q}_1$ differs from $\mathbf{q}_1$ by an angle $\eta$ , then the matrix actually formed differs from the intended $\mathbf{A}'$ by $\lambda_1(\mathbf{q}_1\mathbf{q}_1^{*} - \widehat\mathbf{q}_1\widehat\mathbf{q}_1^{*})$ , of norm of order $\lvert\lambda_1\rvert\eta$ , and by the perturbation bounds of Chapters 17 and 20 that error is inherited by every eigenvalue and eigenvector computed afterwards. Deflating $n - 1$ times accumulates $n - 1$ such errors. For one or two extreme eigenvalues deflation is the right tool; for the whole spectrum it is not, and the next section builds the algorithm that is.

Everything in this section has been a single vector chasing a single eigenvector. The next section iterates on a whole basis at once, and the result — the QR algorithm — turns out to be the power method in disguise, its accumulated factors furnishing a $\mathbf{Q}\mathbf{R}$ factorization of $\mathbf{A}^{k}$ , with the shifts of this section reappearing as the device that makes it fast.

24.4.5. Exercises

A. Check your understanding

Exercise (A1)

(a) State the two hypotheses of Theorem 24.4.1 on the eigenvalues and on the starting vector, and say what the rate $\rho$ is.
(b) For a Hermitian $\mathbf{A}$ , at what rate does the eigenvector converge and at what rate does the Rayleigh quotient? Which earlier result explains the relation between the two?
(c) Explain in one sentence why inverse iteration with shift $\mu$ finds the eigenvalue nearest $\mu$ .
(d) True or false: “the shifted matrix in inverse iteration is nearly singular, so the method is unreliable”. Justify your answer.
(e) Give the per-step cost of the power method, of inverse iteration and of Rayleigh quotient iteration, for a dense $n \times n$ matrix.

Solution.

(a) The dominant eigenvalue must be strictly dominant in modulus, $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ , and the starting vector must have a non-zero coefficient $c_1$ on the corresponding eigenvector. The rate is $\rho = \lvert\lambda_2/\lambda_1\rvert$ , the ratio of the two largest moduli.
(b) The eigenvector converges at rate $\rho^{k}$ (Theorem 24.4.1 (a)) and the Rayleigh quotient at rate $\rho^{2k}$ (Corollary 24.4.2). The relation is the second-order accuracy of the Rayleigh quotient at a Hermitian matrix, quantified by Theorem 20.6.2: a residual of size $\varepsilon$ gives an eigenvalue error of at most $\varepsilon^{2}$ divided by the gap.
(c) Because $(\mathbf{A} - \mu\mathbf{I})^{-1}$ has the same eigenvectors as $\mathbf{A}$ and the eigenvalues $(\lambda_i - \mu)^{-1}$ , whose moduli are largest precisely for the $\lambda_i$ nearest $\mu$ (Theorem 24.4.3).
(d) False. The forward error in $\mathbf{z}^{(k)}$ is indeed large, but only its direction is used, and normalizing removes the part of the error that the ill-conditioning creates; the residual $\left\| \mathbf{A}\widehat\mathbf{q} - \mu\widehat\mathbf{q} \right\|$ of the normalized answer is small, and that is what certifies an eigenvector. The method is unreliable for a different reason — a small gap $g$ — as the warning after Theorem 24.4.3 says.
(e) Power method: one matrix–vector product, $2n^2$ flops. Inverse iteration: one factorization of $\mathbf{A} - \mu\mathbf{I}$ by elimination with pivoting, $\tfrac23n^3$ , computed once, then $2n^2$ per step. Rayleigh quotient iteration: a new factorization every step, $\tfrac23n^3$ per step — paid for by the cubic rate, which needs only three or four steps.

B. Practice

Exercise (B1: A rate, read off the spectrum)

Let $\mathbf{A} = \begin{psmallmatrix} 4 & 1 \\ 2 & 3 \end{psmallmatrix}$ .

(a) Find the eigenvalues and eigenvectors of $\mathbf{A}$ , and the rate $\rho$ of the power method.
(b) Starting from $\mathbf{x}^{(0)} = \mathbf{e}_1$ , compute $\mathbf{A}^k\mathbf{e}_1$ in closed form and verify that the direction converges at rate $\rho^{k}$ .
(c) How many steps are needed to gain six decimal digits?

Solution.

(a) $\operatorname{tr}\mathbf{A} = 7$ and $\det\mathbf{A} = 10$ , so the eigenvalues are $5$ and $2$ . For $\lambda = 5$ : $(\mathbf{A} - 5\mathbf{I})\mathbf{v} = \mathbf{0}$ reads $-v_1 + v_2 = 0$ , giving $\mathbf{v}_1 = (1,1)$ . For $\lambda = 2$ : $2v_1 + v_2 = 0$ , giving $\mathbf{v}_2 = (1,-2)$ . Hence $\rho = \tfrac25$ .
(b) $\mathbf{e}_1 = \tfrac23(1,1) + \tfrac13(1,-2)$ , as one checks by adding. So $\mathbf{A}^{k}\mathbf{e}_1 = \tfrac23 5^{k}(1,1) + \tfrac13 2^{k}(1,-2) .$ Dividing by $\tfrac23 5^{k}$ gives $(1,1) + \tfrac12(\tfrac25)^{k}(1,-2)$ , so the distance from the direction $(1,1)$ is a constant times $(\tfrac25)^k = \rho^k$ , as Theorem 24.4.1 says.
(c) Each step multiplies the error by $\rho = 0.4$ , that is by $10^{-\log_{10}(2.5)} = 10^{-0.398}$ . Six digits need $6/0.398 = 15.1$ steps, so sixteen (and indeed $\rho^{15} = 1.1\cdot10^{-6}$ , $\rho^{16} = 4.3\cdot10^{-7}$ ).

Exercise (B2: Choosing a shift)

Let $\mathbf{A} = \operatorname{diag}(10, 3, 1)$ .

(a) What is the rate of the plain power method?
(b) Find the rate of inverse iteration with shift $\mu = 0$ , and say which eigenvector it finds.
(c) Find the rate of inverse iteration with shift $\mu = \tfrac52$ , and say which eigenvector it finds.
(d) Which shift would you choose to compute the middle eigenvalue as fast as possible, and what rate does it give? Explain why no shift gives a rate of $0$ .

Solution.

(a) By Theorem 24.4.1, $\rho = \lvert 3/10\rvert = 0.3$ , and the iterates converge to $\mathbf{e}_1$ .
(b) With $\mu = 0$ the eigenvalue nearest $0$ is $1$ , at distance $1$ ; the next nearest is $3$ , at distance $3$ . So $\rho_{\mu} = \tfrac13$ , and the method finds $\mathbf{e}_3$ . (Inverse iteration with shift $0$ is the power method on $\mathbf{A}^{-1}$ .)
(c) Distances from $\tfrac52$ are $7.5,\ 0.5,\ 1.5$ . The nearest eigenvalue is $3$ , the next nearest $1$ , so $\rho_{\mu} = 0.5/1.5 = \tfrac13$ , and the method finds $\mathbf{e}_2$ , the middle eigenvector — which the plain power method can never reach.
(d) Take $\mu = 3 + \delta$ with $0 < \delta < 1$ : the distances are $7 - \delta$ , $\delta$ , $2 + \delta$ , so $\rho_{\mu} = \delta/(2+\delta) \to 0$ as $\delta \to 0$ . Any $\mu$ very close to $3$ is excellent. No shift gives rate exactly $0$ , because $\mu = 3$ is an eigenvalue and Definition 24.4.2 then does not apply: $\mathbf{A} - 3\mathbf{I}$ is singular and there is nothing to solve with. The limiting case is not a fast method but no method — which is why Rayleigh quotient iteration, which walks the shift towards the eigenvalue, stops the moment it arrives.

Exercise (B3: Deflating a symmetric matrix)

Let $\mathbf{A} = \begin{psmallmatrix} 2 & 1 \\ 1 & 2\end{psmallmatrix}$ .

(a) Write down $\mathbf{A}' = \mathbf{A} - 3\mathbf{q}_1\mathbf{q}_1^{\top}$ for the unit eigenvector $\mathbf{q}_1$ of the eigenvalue $3$ , and verify its eigenvalues directly.
(b) What does the power method of Theorem 24.4.1, applied to $\mathbf{A}'$ from $\mathbf{x}^{(0)} = \mathbf{e}_1$ , produce, and at what rate?
(c) Explain why the answer to (b) makes deflation look better than it is.

Solution.

(a) $\mathbf{q}_1 = \tfrac{1}{\sqrt2}(1,1)$ and $\mathbf{q}_1\mathbf{q}_1^{\top} = \tfrac12\begin{psmallmatrix}1&1\\mathbf{1}&1\end{psmallmatrix}$ , so $\mathbf{A}' = \begin{pmatrix} 2 & 1 \\ 1 & 2\end{pmatrix} - \tfrac32\begin{pmatrix}1&1\\mathbf{1}&1\end{pmatrix} = \tfrac12\begin{pmatrix} 1 & -1 \\ -1 & 1\end{pmatrix} .$ Then $\mathbf{A}'(1,1) = \mathbf{0}$ and $\mathbf{A}'(1,-1) = (1,-1)$ : the eigenvalues are $0$ and $1$ , as predicted, with the same eigenvectors.
(b) $\mathbf{A}'\mathbf{e}_1 = \tfrac12(1,-1)$ , which is already an eigenvector, so the power iteration of Definition 24.4.1 gives $\mathbf{x}^{(1)} = \tfrac1{\sqrt2}(1,-1)$ and every later iterate equals it: it converges in one step, with $\mu_1 = 1$ exactly.
(c) Because the deflation was exact. The eigenvector $\mathbf{q}_1 = \tfrac1{\sqrt2}(1,1)$ was known in closed form, so $\mathbf{A}'$ had the eigenvalue $0$ exactly and the remaining eigenvector survived untouched. In a computation $\mathbf{q}_1$ is known only to within the accuracy of the first run, and by the warning after the deflation formula the error in $\mathbf{q}_1$ enters $\mathbf{A}'$ as a perturbation of size about $\lvert\lambda_1\rvert\eta$ , which every later eigenpair inherits. One clean step here is not evidence about $n - 1$ dirty ones.

C. Going deeper

Exercise (C1: What the squaring really needs)

Let $\mathbf{A} \in M_n(F)$ be Hermitian with orthonormal eigenbasis $\mathbf{q}_1, \dots, \mathbf{q}_n$ , and let $\mathbf{x}$ be a unit vector with $s = \left\| \mathbf{x} - \mathbf{q}_1\mathbf{q}_1^{*}\mathbf{x} \right\|$ .

(a) Prove that $\lvert R_{\mathbf{A}}(\mathbf{x}) - \lambda_1\rvert \le s^{2}\max_{i\ne1}\lvert\lambda_i - \lambda_1\rvert$ .
(b) Show by example that the corresponding statement for a non-Hermitian $\mathbf{A}$ is false: exhibit $\mathbf{A} \in M_2(\mathbb{C})$ , an eigenvalue $\lambda_1$ with unit eigenvector $\mathbf{q}_1$ , and unit vectors $\mathbf{x}$ with $s \to 0$ for which $\lvert\mathbf{x}^{*}\mathbf{A}\mathbf{x} - \lambda_1\rvert$ is not $O(s^2)$ .
(c) Hence say which hypothesis of Theorem 24.4.4 the example violates, and which step of the proof breaks.

Solution.

(a) This is Step 1 of the proof of Theorem 24.4.4 with $j = 1$ . Expanding $\mathbf{x} = \sum_i c_i\mathbf{q}_i$ with $\sum_i\lvert c_i\rvert^2 = 1$ , we get $s^2 = \sum_{i\ne1}\lvert c_i\rvert^2$ and, by Theorem 11.2.1 (b), $R_{\mathbf{A}}(\mathbf{x}) = \sum_i\lvert c_i\rvert^2\lambda_i$ . Since the weights sum to $1$ , $R_{\mathbf{A}}(\mathbf{x}) - \lambda_1 = \sum_{i\ne1}\lvert c_i\rvert^{2}(\lambda_i - \lambda_1),$ whose modulus is at most $s^2\max_{i\ne1}\lvert\lambda_i-\lambda_1\rvert$ .
(b) Take $\mathbf{A} = \begin{psmallmatrix} 0 & 1 \\ 0 & 0\end{psmallmatrix}$ , the standard non-diagonalizable example, with $\lambda_1 = 0$ and $\mathbf{q}_1 = \mathbf{e}_1$ . Let $\mathbf{x} = (\sqrt{1-s^2},\, s)$ for small $s > 0$ , a unit vector whose part orthogonal to $\mathbf{e}_1$ has length exactly $s$ . Then $\mathbf{x}^{*}\mathbf{A}\mathbf{x} = \overline{x_1}x_2 = s\sqrt{1-s^{2}} ,$ which is $s + O(s^3)$ : first order in $s$ , not second. So no bound of the form $Cs^2$ can hold as $s \to 0$ .
(c) The example violates the hypothesis that $\mathbf{A}$ be Hermitian. The step that breaks is Step 1: for a Hermitian $\mathbf{A}$ the quantity $\mathbf{x}^{*}\mathbf{A}\mathbf{x}$ is the weighted average $\sum_i\lvert c_i\rvert^2\lambda_i$ of the eigenvalues, in which only the squared moduli of the coefficients appear and the first-order term is therefore absent. For a general $\mathbf{A}$ the cross terms $\overline{c_i}c_k\,\mathbf{q}_i^{*}\mathbf{A}\mathbf{q}_k$ with $i \ne k$ survive, and they are of first order in the small coefficients. With Step 1 weakened to $O(s)$ , Step 3 multiplies $s$ by $O(s)$ and the rate is quadratic, not cubic — which is the remark after Theorem 24.4.4.

Exercise (C2: Two eigenvalues of equal modulus, over the reals)

Let $\mathbf{A} \in M_2(\mathbb{R})$ be the rotation matrix $\begin{psmallmatrix} 0 & -1 \\ 1 & 0\end{psmallmatrix}$ .

(a) Find the eigenvalues of $\mathbf{A}$ and explain why Theorem 24.4.1 does not apply.
(b) Describe the power iteration from an arbitrary unit $\mathbf{x}^{(0)} \in \mathbb{R}^2$ , and show that no line is a limit.
(c) Now take the shift $\mu = \tfrac12$ and run inverse iteration over $\mathbb{C}$ . Does the theorem apply? Compute $\rho_{\mu}$ . Hint: the eigenvalues of $\mathbf{A}$ are $\pm i$ .
(d) Deduce a one-sentence moral about which of the two failures of Example (The method that never converges) a shift can repair.

Solution.

(a) $p_{\mathbf{A}}(x) = x^2 + 1$ , so the eigenvalues are $i$ and $-i$ , both of modulus $1$ . The hypothesis $\lvert\lambda_1\rvert > \lvert\lambda_2\rvert$ fails.
(b) $\mathbf{A}$ is orthogonal, so $\left\| \mathbf{A}\mathbf{x} \right\| = \left\| \mathbf{x} \right\|$ and no normalization is needed: $\mathbf{x}^{(k)} = \mathbf{A}^k\mathbf{x}^{(0)}$ . Since $\mathbf{A}^2 = -\mathbf{I}$ , the iterates are $\mathbf{x}^{(0)}, \mathbf{A}\mathbf{x}^{(0)}, -\mathbf{x}^{(0)}, -\mathbf{A}\mathbf{x}^{(0)}, \mathbf{x}^{(0)}, \dots$ , of period $4$ , and the lines they span alternate between $\operatorname{span}(\mathbf{x}^{(0)})$ and $\operatorname{span}(\mathbf{A}\mathbf{x}^{(0)})$ , which are distinct because $\mathbf{A}$ has no real eigenvector. A sequence taking two distinct values infinitely often has no limit.
(c) Over $\mathbb{C}$ the matrix is diagonalizable, with eigenvectors $(1, i)$ and $(1, -i)$ . The shift $\mu = \tfrac12$ is not an eigenvalue, and the distances are $\lvert i - \tfrac12\rvert = \lvert -i - \tfrac12\rvert = \sqrt{5}/2$ — equal. So there is no unique nearest eigenvalue and Theorem 24.4.3 still does not apply; $\rho_{\mu} = 1$ . Any real shift fails in the same way, by symmetry in the conjugate pair; a shift such as $\mu = \tfrac{i}{2}$ does work, with distances $\tfrac12$ and $\tfrac32$ and $\rho_{\mu} = \tfrac13$ .
(d) A shift repairs a bad rate and a wrong target, by moving the origin from which moduli are measured; it cannot repair a tie, and over $\mathbb{R}$ a conjugate pair of eigenvalues is a tie that no real shift can break.

Exercise (C3: A small residual certifies an eigenpair)

Let $\mathbf{A} \in M_n(\mathbb{C})$ be Hermitian with orthonormal eigenbasis $\mathbf{q}_1, \dots, \mathbf{q}_n$ , let $\widehat{\mathbf{q}}$ be a unit vector, let $\widehat\mu = R_{\mathbf{A}}(\widehat{\mathbf{q}})$ and $\varepsilon = \left\| \mathbf{A}\widehat{\mathbf{q}} - \widehat\mu\widehat{\mathbf{q}} \right\|$ . Suppose $\lambda_j$ is the only eigenvalue of $\mathbf{A}$ in the interval $(\widehat\mu - g, \widehat\mu + g)$ .

(a) Prove that $\left\| \widehat{\mathbf{q}} - \mathbf{q}_j\mathbf{q}_j^{*}\widehat{\mathbf{q}} \right\| \le \varepsilon/g$ , assuming $\varepsilon < g$ .
(b) Deduce a bound on $\lvert\widehat\mu - \lambda_j\rvert$ from Theorem 20.6.2, and compare it with the bound from Theorem 20.6.1.
(c) In the inverse iteration of Definition 24.4.2 with a backward stable solve, $\varepsilon$ is about $1/\left\| \widehat{\mathbf{z}} \right\| + cu\left\| \mathbf{A} - \mu\mathbf{I} \right\|$ . Explain what (a) then says about the computed direction, and in particular why the accuracy improves as the shift improves.

Solution.

(a) Expand $\widehat{\mathbf{q}} = \sum_i d_i\mathbf{q}_i$ with $\sum_i\lvert d_i\rvert^2 = 1$ . Then $\mathbf{A}\widehat{\mathbf{q}} - \widehat\mu\widehat{\mathbf{q}} = \sum_i d_i(\lambda_i - \widehat\mu)\mathbf{q}_i$ , so by Theorem 11.2.1 (c), $\varepsilon^{2} = \sum_i \lvert d_i\rvert^{2}(\lambda_i - \widehat\mu)^{2} \ \ge\ \sum_{i \ne j}\lvert d_i\rvert^{2}(\lambda_i - \widehat\mu)^{2} \ \ge\ g^{2}\sum_{i\ne j}\lvert d_i\rvert^{2} ,$ because every $\lambda_i$ with $i \ne j$ lies outside $(\widehat\mu - g, \widehat\mu + g)$ . Since $\left\| \widehat{\mathbf{q}} - \mathbf{q}_j\mathbf{q}_j^{*}\widehat{\mathbf{q}} \right\|^2 = \sum_{i\ne j}\lvert d_i\rvert^2$ , taking square roots gives the bound.
(b) Theorem 20.6.1 gives an eigenvalue within $\varepsilon$ of $\widehat\mu$ , and by hypothesis the only candidate in that range (once $\varepsilon < g$ ) is $\lambda_j$ , so $\lvert\widehat\mu - \lambda_j\rvert \le \varepsilon$ . Theorem 20.6.2, applied with $\alpha = \widehat\mu - g$ and $\beta = \widehat\mu + g$ , improves this to $\lvert\widehat\mu - \lambda_j\rvert \le \frac{\varepsilon^{2}}{g} ,$ which is smaller by a factor $\varepsilon/g$ — the quadratic gain of the Rayleigh quotient once more.
(c) It says that the direction of the computed $\widehat{\mathbf{z}}$ is within $\varepsilon/g$ of $\mathbf{q}_j$ , whatever the relative error in $\widehat{\mathbf{z}}$ itself. As the shift $\mu$ approaches $\lambda_j$ , the matrix $\mathbf{A} - \mu\mathbf{I}$ becomes more nearly singular, and — provided the right-hand side $\mathbf{x}^{(k)}$ has a component along $\mathbf{q}_j$ that is not small, by the bound $\left\| \widehat{\mathbf{z}} \right\| \ge \lvert\mathbf{q}_j^{*}\mathbf{x}^{(k)}\rvert/(\lvert\lambda_j - \mu\rvert + \left\| \boldsymbol{\Delta} \right\|)$ established in the discussion following Theorem 24.4.3 — $\left\| \widehat{\mathbf{z}} \right\|$ grows like $1/\lvert\lambda_j - \mu\rvert$ . So the first part of $\varepsilon$ shrinks in proportion, and the second part is already at the rounding level. A worse-conditioned system gives a better answer, because the quantity that grows is the one in the denominator. What the bound never forgives is a small $g$ : two nearby eigenvalues leave $\varepsilon/g$ large no matter how good the shift is.