Power iteration, inverse iteration, and Rayleigh quotient iteration

Intuition Beginner

A matrix has a set of special directions, the eigenvectors, and each one comes with a number, its eigenvalue, that says how much the matrix stretches that direction. The largest eigenvalue belongs to the direction the matrix stretches the most. A natural question is how to find that dominant direction and its stretch factor without solving for every eigenvalue at once.

Power iteration answers with a one-line idea: pick any starting vector, multiply it by the matrix, and repeat. Each multiplication amplifies the dominant direction more than any other, because that direction is stretched the most. After enough rounds the vector points almost exactly along the dominant eigenvector, and you rescale it each step so it does not grow without bound. The leftover directions fade away at a rate set by how much the largest stretch factor beats the second largest.

Inverse iteration flips the trick to target any eigenvalue you like. You pick a guess near the eigenvalue you want, called a shift, and apply the inverse of the matrix minus that guess. This new matrix stretches most along the eigenvector whose eigenvalue sits closest to your guess, so the same repeated-multiplication idea now homes in on that one. Rayleigh quotient iteration takes the next step: after each round it uses the current vector to form a sharper guess, and for symmetric matrices this feedback makes the accuracy improve at a breathtaking pace.

Visual Beginner

The picture shows a starting vector being multiplied by a matrix again and again. At first the vector points in some arbitrary direction. After one multiplication it tilts toward the dominant eigendirection; after two it is closer still; after a handful of steps it sits almost exactly along that direction. The off-direction part shrinks by a constant factor every step, so the convergence looks like a vector swinging steadily into place.

Two facts are visible. The vector turns toward the dominant direction and the angle to it shrinks by the same ratio every step, the ratio of the second stretch factor to the first. If those two factors were close together the curve would drop slowly, and if the first dwarfs the second the curve would plunge fast.

Worked example Beginner

Take the symmetric matrix

$$A=\left(\begin{array}{cc}3& 1\\ 1& 3\end{array}\right),$$

whose eigenvalues are $4$ (along the direction $(1,1)$) and $2$ (along $(1,-1)$). Start power iteration from the vector $x_0=(1,0)$, which is not aligned with either eigendirection.

Step 1. Multiply: $A{x}_0=(3,1)$. Rescale by dividing by the larger entry $3$ to keep numbers tidy: $x_1=(1,0.333)$.

Step 2. Multiply: $A{x}_1=(3+0.333,1+1)=(3.333,2)$. Rescale by $3.333$: $x_2=(1,0.6)$.

Step 3. Multiply: $A{x}_2=(3+0.6,1+1.8)=(3.6,2.8)$. Rescale by $3.6$: $x_3=(1,0.778)$.

Step 4. Multiply: $A{x}_3=(3+0.778,1+2.333)=(3.778,3.333)$. Rescale by $3.778$: $x_4=(1,0.882)$.

The second coordinate is climbing toward $1$, which is the dominant eigendirection $(1,1)$. The gap to $1$ shrinks roughly by half each step, matching the ratio $2/4=0.5$ of the two eigenvalues.

Step 5. Estimate the eigenvalue. Take the Rayleigh quotient of $x_4=(1,0.882)$: compute $A{x}_4=(3.882,3.647)$, then the ratio $(A{x}_4)\cdot x_4/(x_4\cdot x_4)=(3.882+3.216)/(1+0.778)=7.098/1.778=3.99$.

What this tells us: a few rounds of plain matrix multiplication drove an arbitrary vector almost onto the dominant eigenvector, and the Rayleigh quotient turned that vector into the estimate $3.99$ for the dominant eigenvalue $4$. No characteristic polynomial was solved; the matrix did the work by stretching.

Check your understanding Beginner

Formal definition Intermediate+

Let $A\in M_n(\mathbb{F})$ with $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$, and suppose $A$ is diagonalisable with eigenvalues ${\lambda }_1,\dots ,{\lambda }_n$ and a basis of eigenvectors $u_1,\dots ,u_n$, where $A{u}_i={\lambda }_i{u}_i$ and $\Vert u_i\Vert =1$ in the Euclidean norm of 01.01.08. Order the eigenvalues by modulus, $|{\lambda }_1|\ge |{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$. The eigenvalue ${\lambda }_1$ is dominant when $|{\lambda }_1|>|{\lambda }_2|$; in that strict case the dominant eigenspace is one-dimensional and $u_1$ is determined up to a scalar.

Definition (power iteration). Given a start vector $x_0$ with $\Vert x_0\Vert =1$, power iteration generates the sequence

$$y_k=A{x}_{k-1},x_k=\frac{y_k}{\Vert y_k\Vert },k=1,2,\dots ,$$

together with the eigenvalue estimate by the Rayleigh quotient ${\rho }_k=R(x_k)=⟨A{x}_k,x_k⟩/⟨x_k,x_k⟩$, the optimal scalar estimate of the eigenvalue associated with the direction $x_k$, in the sense of 01.01.14. Each step is one matrix-vector product followed by a normalisation.

Definition (inverse iteration with a shift). Fix a shift $\mu \in \mathbb{F}$ that is not an eigenvalue, so $A-\mu I$ is invertible. Inverse iteration generates

$$(A-\mu I)y_k=x_{k-1},x_k=\frac{y_k}{\Vert y_k\Vert },$$

that is, power iteration applied to the matrix $(A-\mu I{)}^{-1}$. The vector $y_k$ is computed by solving a linear system with the fixed matrix $A-\mu I$ — a single factorisation reused across iterations — rather than by forming the inverse explicitly.

Definition (Rayleigh quotient iteration). Replace the fixed shift by the current Rayleigh quotient at every step. Starting from $x_0$ with $\Vert x_0\Vert =1$ and ${\mu }_0=R(x_0)$,

$$(A-{\mu }_{k-1}I)y_k=x_{k-1},x_k=\frac{y_k}{\Vert y_k\Vert },{\mu }_k=R(x_k).$$

The shift ${\mu }_{k-1}$ now varies, so each step solves a system with a different matrix, but the adaptively improving shift produces convergence far faster than a fixed shift.

Notation: $\angle (v,w)$ is the acute angle between the lines spanned by nonzero $v$ and $w$, with $\sin \angle (v,w)$ the standard measure of misalignment; $\sigma (A)$ is the spectrum; $R(v)=⟨Av,v⟩/⟨v,v⟩$ is the Rayleigh quotient. Throughout, "convergence of $x_k$ to $u$" means $\sin \angle (x_k,u)\to 0$, since an eigenvector is determined only up to a scalar.

Counterexamples to common slips

• Power iteration needs a strictly dominant eigenvalue. For $A=\mathrm{diag}(1,-1)$ the two eigenvalues have equal modulus, and a generic start vector oscillates between two directions without settling; the ratio $|{\lambda }_2/{\lambda }_1|=1$ gives no contraction.
• The method fails if the start vector has zero component along $u_1$. For $A=\mathrm{diag}(2,1)$ and $x_0=(0,1)$, every iterate stays $(0,1)$ and converges to the subdominant eigenvector. In exact arithmetic this is a measure-zero accident; in floating point, rounding usually reintroduces a $u_1$-component and the method recovers.
• A complex dominant pair $\lambda ,\bar{\lambda }$ of a real matrix has equal modulus, so real power iteration does not converge to a single vector; the dominant invariant subspace is two-dimensional and requires a block method.

Key theorem with proof Intermediate+

Theorem (linear convergence of power iteration; Trefethen-Bau Lecture 27 ^{[source pending]}; Golub-Van Loan §8.2 ^{[source pending]}). Let $A\in M_n(\mathbb{F})$ be diagonalisable with eigenpairs $({\lambda }_i,u_i)$, $\Vert u_i\Vert =1$, and a strictly dominant eigenvalue, $|{\lambda }_1|>|{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$. Write the start vector as $x_0={\sum }_{i=1}^n{c}_i{u}_i$ and assume $c_1\ne 0$. Then the power-iteration sequence $x_k=A^k{x}_0/\Vert A^k{x}_0\Vert$ satisfies

$$\sin \angle (x_k,u_1)=O\left({\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k\right)\text{as }k\to \infty ,$$

and the Rayleigh quotient ${\rho }_k=R(x_k)$ converges to ${\lambda }_1$. If in addition $A$ is self-adjoint, the eigenvalue error sharpens to ${\rho }_k-{\lambda }_1=O\left(|{\lambda }_2/{\lambda }_1{|}^{2k}\right)$.

Proof. Expand in the eigenbasis: $A^k{x}_0={\sum }_{i=1}^n{c}_i{\lambda }_i^k{u}_i$. Factor out the dominant term,

$$A^k{x}_0=c_1{\lambda }_1^k\left(u_1+\sum _{i=2}^n\frac{c_i}{c_1}{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{u}_i\right)=c_1{\lambda }_1^k\left(u_1+r_k\right),$$

where the remainder $r_k={\sum }_{i\ge 2}(c_i/c_1)({\lambda }_i/{\lambda }_1{)}^k{u}_i$ has norm bounded by

$$\Vert r_k\Vert \le \sum _{i\ge 2}\left|\frac{c_i}{c_1}\right|{\left|\frac{{\lambda }_i}{{\lambda }_1}\right|}^k\le \left(\sum _{i\ge 2}\left|\frac{c_i}{c_1}\right|\right){\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k=C{\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k,$$

because $|{\lambda }_i/{\lambda }_1|\le |{\lambda }_2/{\lambda }_1|$ for every $i\ge 2$. Normalising kills the scalar $c_1{\lambda }_1^k$, so $x_k$ is a unit vector in the direction of $u_1+r_k$. The sine of the angle to $u_1$ is the norm of the component of $x_k$ orthogonal to $u_1$, which is at most $\Vert r_k\Vert /\Vert u_1+r_k\Vert$; since $\Vert u_1+r_k\Vert \to 1$, this is $O(\Vert r_k\Vert )=O(|{\lambda }_2/{\lambda }_1{|}^k)$.

For the eigenvalue estimate, the Rayleigh quotient is continuous and $R(u_1)={\lambda }_1$, so ${\rho }_k=R(x_k)\to {\lambda }_1$ as $x_k\to u_1$ in direction. When $A=A^{\ast }$, the eigenbasis is orthonormal, so $r_k\perp u_1$ to leading order and $R(x_k)-{\lambda }_1$ is a quadratic form in the small quantity $r_k$: writing $x_k\propto u_1+r_k$ with $r_k\perp u_1$, one computes $R(x_k)-{\lambda }_1=⟨(A-{\lambda }_{1}I)r_k,r_k⟩/(1+\Vert r_k{\Vert }^2)$, which is $O(\Vert r_k{\Vert }^2)=O(|{\lambda }_2/{\lambda }_1{|}^{2k})$. This is the same second-order stationarity that makes the Rayleigh quotient the optimal eigenvalue estimate at an approximate eigenvector. $\square$

Bridge. This convergence analysis builds toward the inverse-iteration and Rayleigh-quotient-iteration theorems of the Advanced results, where the very same eigenbasis expansion is applied to $(A-\mu I{)}^{-1}$, whose eigenvalues $({\lambda }_i-\mu {)}^{-1}$ are dominated by the one nearest the shift; the ratio $|{\lambda }_2/{\lambda }_1|$ that governs the rate here appears again in 43.06.03 as the convergence engine of the QR algorithm, which is simultaneous power iteration on an orthonormal basis. The foundational reason the method works is the spectral expansion $A^k{x}_0={\sum }_i{c}_i{\lambda }_i^k{u}_i$: raising $A$ to a power raises each eigenvalue to that power, and the largest in modulus eventually swamps the rest. This is exactly the dynamical content of the dominant eigenvalue, and the quadratic sharpening in the self-adjoint case generalises the stationarity of the Rayleigh quotient established in 01.01.14. Putting these together, plain matrix multiplication is a convergence process whose speed is read off the spectral gap, and the bridge is the observation that a faster method must reshape the spectrum so that the targeted eigenvalue dominates more sharply — precisely what shifting and inverting accomplish.

Exercises Intermediate+

Advanced results Master

Theorem (inverse iteration converges to the eigenvalue nearest the shift; Wilkinson Ch. 9 ^{[source pending]}; Trefethen-Bau Lecture 27 ^{[source pending]}). Let $A\in M_n(\mathbb{F})$ be diagonalisable with eigenpairs $({\lambda }_i,u_i)$, and let $\mu \notin \sigma (A)$ be a shift such that a unique eigenvalue ${\lambda }_J$ minimises $|{\lambda }_i-\mu |$, with the runner-up distance $\delta ={\min }_{i\ne J}|{\lambda }_i-\mu |>|{\lambda }_J-\mu |$. For any start vector with nonzero $u_J$-component, inverse iteration $x_k=(A-\mu I{)}^{-1}{x}_{k-1}/\Vert \cdots \Vert$ satisfies

$$\sin \angle (x_k,u_J)=O\left({\left(\frac{|{\lambda }_J-\mu |}{\delta }\right)}^k\right).$$

The mechanism is the eigenvalue transformation ${\lambda }_i\mapsto ({\lambda }_i-\mu {)}^{-1}$, which makes the eigenvalue nearest $\mu$ the dominant eigenvalue of $(A-\mu I{)}^{-1}$; the convergence theorem of the Intermediate tier then applies verbatim to $B=(A-\mu I{)}^{-1}$. The closer the shift is to ${\lambda }_J$, the smaller the ratio $|{\lambda }_J-\mu |/\delta$ and the faster the convergence — the practical reason inverse iteration is the standard tool for refining an eigenvector once an approximate eigenvalue is known. The system $(A-\mu I)y_k=x_{k-1}$ becomes severely ill-conditioned as $\mu \to {\lambda }_J$, yet Wilkinson's analysis shows the computed $y_k$ points accurately in the direction of $u_J$: the large error introduced by the near-singular solve lies almost entirely along $u_J$ itself, so after normalisation it does no harm.

Theorem (cubic convergence of Rayleigh quotient iteration for self-adjoint $A$; Parlett §4.6 ^{[source pending]}). Let $A = A^$haveasimpleeigenvalue$\lambda$withuniteigenvector$u$.Thereisaneighbourhoodof$u$suchthat,startedfromany$x_0$initwith$\sin\angle(x_0, u)$ small, Rayleigh quotient iteration produces iterates with*

$$\sin \angle (x_{k+1},u)=O\left({\sin }^3\angle (x_k,u)\right).$$

For a general (non-normal) diagonalisable $A$ the convergence is quadratic, $\sin \angle (x_{k+1},u)=O({\sin }^2\angle (x_k,u))$.

Cubic convergence means the number of correct digits roughly triples per step; in practice Rayleigh quotient iteration on a symmetric matrix reaches machine precision in three or four iterations once it enters the basin of attraction. The cube — rather than the square — is bought by the second-order accuracy of the Rayleigh quotient: when $A=A^{\ast }$, the shift error $|{\mu }_k-\lambda |=O({\sin }^2\angle (x_k,u))$ is one power smaller than for a generic shift, and inverse iteration with that shift contracts the eigenvector error by an extra factor of the shift error, compounding the second-order shift accuracy with the first-order inverse-iteration contraction to give the third power.

Theorem (deflation and the full symmetric spectrum). Let $A = A^$withorthonormaleigenbasis$u_1, \ldots, u_n$andeigenvaluesorderedbymodulus.Havingcomputed$(\lambda_1, u_1)$,thedeflatedoperator$A^{(1)} = A - \lambda_1 u_1 u_1^$hasspectrum${0} \cup {\lambda_2, \ldots, \lambda_n}$ with the same eigenvectors; iterating the deflation recovers the entire spectrum one eigenpair at a time. The deflation is the Hermitian rank-one correction that annihilates the found eigenvalue while preserving orthogonal eigendirections, and it is the conceptual seed of the block and subspace methods that compute several eigenpairs at once. Numerically, explicit deflation accumulates rounding error across eigenpairs, which is why the production algorithm of the section is instead the QR algorithm of 43.06.03, computing the full spectrum simultaneously and stably.

Theorem (shifted power iteration and spectral relocation). For any scalar $s$, the eigenvalues of $A-sI$ are ${\lambda }_i-s$ with unchanged eigenvectors. Choosing $s$ to maximise the gap $|{\lambda }_1-s|/|{\lambda }_2-s|$ accelerates power iteration without changing the target eigenvector; for a self-adjoint $A$ with spectrum in $[\alpha ,\beta ]$ and dominant ${\lambda }_1=\beta$, the optimal real shift is $s=(\alpha +{\lambda }_2)/2$, centring the unwanted spectrum about zero. This spectral-relocation principle is the bridge from the single-vector methods of this unit to polynomial acceleration: replacing $A$ by a polynomial $p(A)$ chosen to amplify the wanted eigenvalue and damp the rest is the idea behind Chebyshev acceleration and, ultimately, the Krylov methods of chapter $43.07$.

Synthesis. The four algorithms of this unit are one idea viewed through successively sharper lenses, and the spectral expansion $A^k{x}_0={\sum }_i{c}_i{\lambda }_i^k{u}_i$ is the foundational reason they work: power iteration amplifies the eigenvalue of largest modulus, inverse iteration reshapes the spectrum by ${\lambda }_i\mapsto ({\lambda }_i-\mu {)}^{-1}$ so that the eigenvalue nearest a chosen shift becomes dominant, and Rayleigh quotient iteration closes the loop by feeding the current Rayleigh-quotient estimate back as the shift. This is exactly where the variational theory of 01.01.14 pays off: the Rayleigh quotient is the optimal eigenvalue estimate at an approximate eigenvector, stationary there, so for a self-adjoint operator its error is second order in the eigenvector error, and that second-order shift accuracy compounds with the first-order contraction of inverse iteration to produce cubic convergence. The convergence rate $|{\lambda }_2/{\lambda }_1|$ that governs power iteration generalises to the gap ratio of any reshaped spectrum, and the same ratio appears again as the engine of the QR algorithm, which is simultaneous orthonormal power iteration.

Putting these together, the single-vector eigenvalue algorithms occupy the place between the static spectral theory — eigenvalues as roots, as critical values of the Rayleigh quotient, as min-max over subspaces — and the production eigensolvers: the central insight is that a good eigenvalue algorithm is a dynamical system whose fixed points are the eigenvectors and whose contraction rate is set by how sharply the targeted eigenvalue can be made to dominate the rest. Deflation and shifting are the two levers on that domination, and the bridge is the recognition that the QR algorithm and the Krylov methods are these same levers applied to a whole basis at once.

Full proof set Master

Proposition (power iteration, full statement and rate). Under the hypotheses of the Intermediate theorem — $A$ diagonalisable, $|{\lambda }_1|>|{\lambda }_2|\ge \cdots$, start vector $x_0={\sum }_i{c}_i{u}_i$ with $c_1\ne 0$ — the normalised iterates satisfy $\sin \angle (x_k,u_1)\le C|{\lambda }_2/{\lambda }_1{|}^k$ for a constant $C$ depending only on the $c_i$ and the eigenbasis conditioning.

Proof. From $A^k{x}_0=c_1{\lambda }_1^k(u_1+r_k)$ with $r_k={\sum }_{i\ge 2}(c_i/c_1)({\lambda }_i/{\lambda }_1{)}^k{u}_i$, decompose any unit iterate $x_k$ as its $u_1$-component plus an orthogonal remainder. Writing $P_1$ for the spectral projection onto $u_1$ along the other eigenvectors, $x_k$ is the unit vector in the direction of $u_1+r_k$, so

$$\sin \angle (x_k,u_1)=\frac{\Vert (I-P_1)(u_1+r_k)\Vert }{\Vert u_1+r_k\Vert }=\frac{\Vert (I-P_1)r_k\Vert }{\Vert u_1+r_k\Vert },$$

since $(I-P_1)u_1=0$. The numerator is at most $\Vert I-P_1\Vert \Vert r_k\Vert$, and $\Vert r_k\Vert \le C_0|{\lambda }_2/{\lambda }_1{|}^k$ from the Intermediate bound. For $k$ large, $\Vert u_1+r_k\Vert \ge 1-\Vert r_k\Vert \ge 1/2$. Hence $\sin \angle (x_k,u_1)\le 2\Vert I-P_1\Vert C_0|{\lambda }_2/{\lambda }_1{|}^k=C|{\lambda }_2/{\lambda }_1{|}^k$. The projection norm $\Vert I-P_1\Vert$ is the eigenbasis conditioning, equal to $1$ when $A$ is normal and larger for non-normal $A$. $\square$

Proposition (inverse iteration rate). With $\mu \notin \sigma (A)$ and a unique nearest eigenvalue ${\lambda }_J$, inverse iteration converges to $u_J$ at rate $|{\lambda }_J-\mu |/\delta$, where $\delta ={\min }_{i\ne J}|{\lambda }_i-\mu |$.

Proof. Set $B=(A-\mu I{)}^{-1}$. By the computation $B{u}_i=({\lambda }_i-\mu {)}^{-1}{u}_i$, the eigenvalues of $B$ are ${\beta }_i=({\lambda }_i-\mu {)}^{-1}$ with the eigenvectors of $A$. The largest in modulus is ${\beta }_J$ because $|{\beta }_i|=1/|{\lambda }_i-\mu |$ is largest when $|{\lambda }_i-\mu |$ is smallest, which the hypothesis pins to $i=J$. The second-largest is $1/\delta$. Inverse iteration is power iteration on $B$, so the previous proposition gives $\sin \angle (x_k,u_J)\le C|{\beta }_2(B)/{\beta }_1(B){|}^k=C(|{\lambda }_J-\mu |/\delta {)}^k$, where ${\beta }_1(B)={\beta }_J$ and ${\beta }_2(B)$ is the runner-up. The uniqueness of the nearest eigenvalue ensures $B$ has a strictly dominant eigenvalue, the condition the power-iteration proposition requires. $\square$

Proposition (cubic convergence of Rayleigh quotient iteration, self-adjoint case). Let $A = A^$havesimpleeigenvalue$\lambda$withuniteigenvector$u$andspectralgap$\gamma = \mathrm{dist}(\lambda, \sigma(A) \setminus {\lambda}) > 0$.For$x_k$with$\theta_k = \angle(x_k, u)$small,theRayleighquotientiterationsatisfies$\sin\theta_{k+1} = O(\sin^3\theta_k)$.*

Proof. Write $x_k=(\cos {\theta }_k)u+(\sin {\theta }_k)z_k$ with $z_k\perp u$, $\Vert z_k\Vert =1$. By the self-adjoint Rayleigh-quotient expansion (Exercise 7), ${\mu }_k=R(x_k)=\lambda +(⟨A{z}_k,z_k⟩-\lambda ){\sin }^2{\theta }_k$, so the shift error is

$$|{\mu }_k-\lambda |\le (\Vert A\Vert +|\lambda |){\sin }^2{\theta }_k=O({\sin }^2{\theta }_k).$$

The next iterate is $x_{k+1}\propto (A-{\mu }_{k}I{)}^{-1}{x}_k$. Decompose along $u$ and $u^{\perp }$: applying $(A-{\mu }_{k}I{)}^{-1}$ scales the $u$-component $\cos {\theta }_k$ by $(\lambda -{\mu }_k{)}^{-1}$ and the $z_k$-component $\sin {\theta }_k$ by an operator of norm at most $(\gamma -|{\mu }_k-\lambda |{)}^{-1}=O(1)$, because $A$ restricted to the invariant subspace $u^{\perp }$ has spectrum at distance $\ge \gamma$ from $\lambda$ and hence $\ge \gamma -|{\mu }_k-\lambda |$ from ${\mu }_k$. Therefore

$$\tan {\theta }_{k+1}=\frac{\Vert \text{component}\perp u\Vert }{\Vert \text{component along }u\Vert }\le \frac{O(1)\sin {\theta }_k}{|\lambda -{\mu }_k{|}^{-1}\cos {\theta }_k}=O(|{\mu }_k-\lambda |)\tan {\theta }_k=O({\sin }^2{\theta }_k)\tan {\theta }_k.$$

Since $\tan {\theta }_k=O(\sin {\theta }_k)$ for small angles, $\tan {\theta }_{k+1}=O({\sin }^3{\theta }_k)$, and $\sin {\theta }_{k+1}\le \tan {\theta }_{k+1}=O({\sin }^3{\theta }_k)$. Self-adjointness enters twice: to make $u^{\perp }$ an $A$-invariant orthogonal complement, and to make ${\mu }_k-\lambda$ second order in ${\theta }_k$. For non-normal $A$ the shift error is only first order, $|{\mu }_k-\lambda |=O(\sin {\theta }_k)$, and the same argument yields the quadratic rate $\sin {\theta }_{k+1}=O({\sin }^2{\theta }_k)$. $\square$

Proposition (deflation preserves the orthogonal spectrum, self-adjoint case). For $A = A^$withuniteigenvector$u_1$andeigenvalue$\lambda_1$,theoperator$A' = A - \lambda_1 u_1 u_1^$isself-adjointwith$A' u_1 = 0$and$A' u_j = \lambda_j u_j$foreveryeigenvector$u_j \perp u_1$.

Proof. The correction ${\lambda }_1{u}_1{u}_1^{\ast }$ is Hermitian (${\lambda }_1$ real, $(u_1{u}_1^{\ast }{)}^{\ast }=u_1{u}_1^{\ast }$), so $A^{\prime }$ is self-adjoint. The eigenbasis of $A$ is orthonormal, so for $j\ge 2$, $u_1^{\ast }{u}_j=⟨u_j,u_1⟩=0$, giving $A^{\prime }{u}_j=A{u}_j-{\lambda }_1{u}_1(u_1^{\ast }{u}_j)={\lambda }_j{u}_j$. On $u_1$, $A^{\prime }{u}_1={\lambda }_1{u}_1-{\lambda }_1{u}_1(u_1^{\ast }{u}_1)={\lambda }_1{u}_1-{\lambda }_1{u}_1=0$ since $\Vert u_1\Vert =1$. The spectrum is therefore $\{0,{\lambda }_2,\dots ,{\lambda }_n\}$ with eigenvectors unchanged. $\square$

Connections Master

The single-vector eigenvalue algorithms are the dynamical counterpart of the variational eigenvalue theory of 01.01.14: where the Rayleigh quotient there is characterised statically as the optimal eigenvalue estimate and the critical values of $R$, here it is the adaptive shift that drives Rayleigh quotient iteration, and its second-order stationarity at an eigenvector is precisely what upgrades inverse iteration from linear to cubic convergence in the self-adjoint case.

Power iteration is the convergence engine of the QR algorithm of 43.06.03: the unshifted QR iteration is simultaneous power iteration applied to a full orthonormal basis rather than a single vector, and the Wilkinson shift that gives the QR algorithm its cubic local convergence on symmetric tridiagonals is the Rayleigh-quotient shift of this unit transplanted into the orthogonal-similarity setting.

The inverse-iteration step $(A-\mu I)y_k=x_{k-1}$ is a linear solve with a fixed, reused matrix, so its efficiency rests on the direct factorisations of the chapter — Gaussian elimination and the LU factorisation are computed once and applied at every iteration; the preliminary reduction to Hessenberg or tridiagonal form of 43.06.02 is what makes each of those solves cheap, which is why the practical pipeline reduces first and iterates second.

The spectral-relocation idea — replacing $A$ by $A-sI$ or by a polynomial $p(A)$ to reshape which eigenvalue dominates — is the seed of the Krylov-subspace methods, where the iterate lies in $\mathrm{span}\{b,Ab,\dots ,A^{m-1}b\}$ and a polynomial in $A$ is chosen to amplify the wanted part of the spectrum; the conjugate gradient and GMRES convergence bounds rest on exactly the polynomial-acceleration principle introduced here in its single-vector form.

Historical & philosophical context Master

The power method entered numerical practice through the 1929 paper of von Mises and Pollaczek-Geiringer, who set out the Potenzmethode as a systematic iterative procedure for the dominant eigenvalue of the matrices arising in elasticity and structural mechanics, where the characteristic polynomial was intractable to solve directly ^{[von Mises 1929]}. The idea of repeatedly applying a linear operator to expose its dominant mode was older in the analytic theory of integral equations, but von Mises gave it the matrix-iteration form that the digital era inherited.

The shifted inverse variant is due to Helmut Wielandt, who in a 1944 Göttingen aerodynamics report introduced gebrochene Iteration — fractional iteration — applying $(A-\mu I{)}^{-1}$ to target an eigenvalue near a chosen point $\mu$, the technique now called inverse iteration or Wielandt iteration ^{[Wielandt 1944]}. Wilkinson, in The Algebraic Eigenvalue Problem (1965), supplied the rounding-error analysis that resolved the apparent paradox of the method: the linear system $(A-\mu I)y=x$ becomes arbitrarily ill-conditioned as the shift approaches an eigenvalue, yet the computed solution is an excellent eigenvector, because the large error is aligned with the very eigenvector being sought ^{[Wilkinson 1965]}. The Rayleigh quotient as a self-correcting shift, and the resulting cubic convergence for symmetric matrices, was analysed in detail by Parlett, whose The Symmetric Eigenvalue Problem established the cubic local rate and its dependence on the second-order accuracy of the Rayleigh quotient ^{[Parlett 1998]}. These single-vector methods were the state of the art for the dominant eigenpair until the QR algorithm of Francis and Kublanovskaya (1961) made the full spectrum computable at comparable cost.

Bibliography Master

@article{vonMises1929,
  author  = {von Mises, Richard and Pollaczek-Geiringer, Hilda},
  title   = {Praktische Verfahren der Gleichungsaufl{\"o}sung},
  journal = {Zeitschrift f{\"u}r Angewandte Mathematik und Mechanik},
  volume  = {9},
  year    = {1929},
  pages   = {58--77, 152--164}
}

@techreport{Wielandt1944,
  author      = {Wielandt, Helmut},
  title       = {Beitr{\"a}ge zur mathematischen Behandlung komplexer Eigenwertprobleme},
  institution = {Aerodynamische Versuchsanstalt G{\"o}ttingen},
  number      = {Bericht B 44/J/37},
  year        = {1944}
}

@book{Wilkinson1965,
  author    = {Wilkinson, James H.},
  title     = {The Algebraic Eigenvalue Problem},
  publisher = {Oxford University Press},
  address   = {Oxford},
  year      = {1965}
}

@book{Parlett1998,
  author    = {Parlett, Beresford N.},
  title     = {The Symmetric Eigenvalue Problem},
  publisher = {SIAM},
  series    = {Classics in Applied Mathematics},
  year      = {1998}
}

@book{TrefethenBau1997,
  author    = {Trefethen, Lloyd N. and Bau, David},
  title     = {Numerical Linear Algebra},
  publisher = {SIAM},
  address   = {Philadelphia},
  year      = {1997}
}

@book{GolubVanLoan2013,
  author    = {Golub, Gene H. and Van Loan, Charles F.},
  title     = {Matrix Computations},
  edition   = {4th},
  publisher = {Johns Hopkins University Press},
  address   = {Baltimore},
  year      = {2013}
}

@article{Francis1961,
  author  = {Francis, John G. F.},
  title   = {The QR Transformation, I},
  journal = {The Computer Journal},
  volume  = {4},
  number  = {3},
  year    = {1961},
  pages   = {265--271}
}