Numerical Analysis HT24, Power method
Flashcards
Describe the steps in the Power Method for computing the largest eigenvalue of a matrix $A \in \mathbb R^{n \times n}$, and how to establish the eigenvector / eigenvalue from the iterations.
- Select arbitrary $y \in \mathbb R^n$
- Set $x _ 0 = y/ \vert \vert y \vert \vert $
- Then, for $k = 0, 1, \cdots$
- Compute $y _ k = Ax _ k$
- Set $x _ {k + 1} = y _ k / \vert \vert y _ k \vert \vert $
Then $x _ k \to \pm v _ 1$ and $ \vert \vert y _ k \vert \vert \approx \vert \lambda _ 1 \vert $, where $\lambda _ 1$ is the largest eigenvalue of $A$ and $v _ 1$ is the corresponding eigenvector of unit length.
The Power Method for computing the largest eigenvalue of a matrix $A \in \mathbb R^{n \times n}$ works as follows:
- Select arbitrary $y \in \mathbb R^n$
- Set $x _ 0 = y/ \vert \vert y \vert \vert $
- Then, for $k = 0, 1, \cdots$
- Compute $y _ k = Ax _ k$
- Set $x _ {k + 1} = y _ k / \vert \vert y _ k \vert \vert $
Then $x _ k \to \pm v _ 1$ and $ \vert \vert y _ k \vert \vert \to \vert \lambda _ 1 \vert $ where $\lambda _ 1$ is the largest eigenvalue of $A$ and $v _ 1$ is the corresponding eigenvector of unit length. Under the assumptions that:
- The matrix is diagonalisable
- The modulus of the largest eigenvalue is strictly bigger than the modulus of the second largest
Quickly:
- Prove that this method works
- Give an expression that determines how fast the convergence is
- Justify why the same idea works for $A^{-1}$ (the “inverse power method”), and give another expression for how fast the convergence is
- Explain the idea of shifts and why they are more useful in the inverse power method than in the ordinary power method
- Compute $y _ k = Ax _ k$
- Set $x _ {k + 1} = y _ k / \vert \vert y _ k \vert \vert $
Proof it works:
Overall idea: Expand $x _ 0$ in the orthonormal basis of eigenvectors of $A$ (which exists since $A$ is diagonalisable). Then look at $A^k x _ 0$, then relate $x _ k$ and $y _ k$ to $A^k x _ 0$.
This algorithm computes a list $y _ k = \beta _ k A^k x _ 0$ for some set of coefficients $\beta _ k$. Suppose $\{v _ 1, \cdots, v _ n\}$ and $\{\lambda _ 1, \cdots, \lambda _ n\}$ are the eigenvalues of $A$ arranged in descending order of modulus and that the eigenvectors are of unit length. Then $\exists \alpha _ i$ such that
\[x _ 0 = \sum^n _ {i = 1} \alpha _ i v _ i\]Hence
\[\begin{aligned} A^k x _ 0 &= \sum^n _ {i = 1} \alpha _ i \lambda _ i^k v _ i \\ &= \lambda _ 1^k \left[ \alpha _ 1 v _ 1 + \sum^n _ {i = 2} \alpha _ i \left(\frac{\lambda _ i}{\lambda _ 1}\right)^k v _ i \right] \end{aligned}\]Since $\lambda _ 1 > \lambda _ i$ for $i \ge 2$, for $k$ large we have $A^k x _ 0 \approx \alpha _ 1 \lambda^k _ 1 v _ 1$. Hence
\[\begin{aligned} x _ {k} &= \frac{y _ {k-1} }{ \vert \vert y _ {k-1} \vert \vert } \\ &= \frac{\beta _ {k-1} }{ \vert \beta _ {k-1} \vert } \cdot \frac{A^{k-1} x _ 0}{ \vert \vert A^{k-1} x _ 0 \vert \vert } \\ &\approx \frac{\beta _ {k-1} }{ \vert \beta _ {k-1} \vert } \cdot \frac{\lambda _ 1^k \alpha _ 1}{ \vert \lambda _ 1^k \alpha _ 1 \vert } \cdot \frac{v _ 1}{ \vert \vert v \vert \vert _ 1} \\ &= \pm v _ 1 \end{aligned}\]and
\[\begin{aligned} y _ k &= Ax _ {k-1} \\ &\approx \pm A v _ 1 \\ &= \pm \lambda _ 1 v _ 1 \end{aligned}\]so $ \vert \vert y _ k \vert \vert = \vert \lambda _ 1 \vert $ (since the eigenvector basis elements have unit modulus).
(The assumption that the matrix is diagonalisable is not necessary, it just suffices that the dominant eigenvalue of $A$ is unique. The proof that doesn’t assume the matrix is diagonalisable instead uses the basis for the Jordan decomposition of $A$, and analyses the convergence that way).
How fast is convergence? The speed at which $A^k x _ 0$ tends to $\lambda^k _ 1 \alpha _ 1 v _ 1$ depends on the ratio
\[\frac{\lambda _ 2}{\lambda _ 1}\]where $\lambda _ 2$ is the second-largest eigenvalue. Ideally, this is as small as possible.
What about the inverse power method? If we instead use $A^{-1}$, then the eigenvalues are $1/\lambda _ i$ for each of the eigenvalues of $A$. Since the dominant eigenvalue of $A^{-1}$ is then $1/\lambda _ n$, applying the inverse power method to $A^{-1}$ means convergence to $1/\lambda _ n$ and corresponding eigenvector $v _ n$. In this case, the speed of convergence depends on
\[\frac{1/\lambda _ {n-1} }{1/\lambda _ n} = \frac{\lambda _ n}{\lambda _ {n-1}}\]Again, for the best convergence, this should be as small as possible.
Shifts: For any $s \in \mathbb C$, the eigenvalues of $A - sI$ are $\lambda - s$. So if we can pick the shift well and apply the power method to $A - sI$, we could make
\[\frac{\lambda _ {\sigma(2)} - s}{\lambda _ {\sigma(1)} - s}\]small, where $\sigma$ is a permutation such that $ \vert \lambda _ {\sigma(1)} - s \vert > \vert \lambda _ {\sigma(2)} - s \vert \ge \cdots \ge \vert \lambda _ {\sigma(n)} - s \vert $. But without knowing the eigenvalues ahead of time, this is difficult. And in some sense we are not going to be able to get crazy speedups because we’re affecting all of the eigenvalues uniformly, so making one really large will also make the others really large.
However, if we apply the inverse power method to $A - sI$, then the speed of convergence depends on
\[\frac{1/(\lambda _ {\sigma(n-1)} - s)}{1/(\lambda _ {\sigma(n)} - s)} = \frac{\lambda _ {\sigma(n)} - s}{\lambda _ {\sigma(n-1)} - s}\]Even if $s$ is only approximately close to an eigenvalue (by construction of $\sigma$, if $s$ is close to an eigenvalue, then this eigenvalue is $\lambda _ {\sigma(n)}$) then it is possible that $\frac{1}{\lambda _ {\sigma(n)} - s}$ is much bigger than the runner up $1/(\lambda _ {\sigma(n-1)} - s)$.
Bite-sized
When does the basic (unshifted) power method fail to converge?
When the dominant eigenvalue is not unique in modulus — e.g. if $A$ has eigenvalues $\pm \lambda$ with $ \vert \lambda _ 1 \vert = \vert \lambda _ 2 \vert $. Then $A^k x _ 0$ oscillates between two directions rather than aligning with a single eigenvector. The fix is shifts or block extensions; orthogonal/unitary $A$ (all $ \vert \lambda _ i \vert = 1$) is a notorious failure case.
Shifted inverse iteration uses a fixed shift $s$, so every iteration solves $(A - sI)\, y = x _ k$ where the coefficient matrix $A - sI$ stays the same and only the right-hand side $x _ k$ changes. One never forms $(A - sI)^{-1}$. Instead compute a $LU$ factorisation of $A - sI$ once at $O(n^3)$, then reuse it: each iteration is then just two triangular solves (forward with $L$, back with $U$) at $O(n^2)$. A single factorisation suffices precisely because the matrix is constant across iterations, so it can be reused for every (changing) right-hand side. (Rayleigh quotient iteration changes the shift $s _ k$ each step, so the matrix changes and it must refactorise, giving $O(n^3)$ per step in exchange for cubic convergence.)