Lecture 13: Power Iterations

1.Reduction to triangular form

The relevance of eigendecompositions $𝑨=𝑿𝚲{𝑿}^{-1}$ to repeated application of the linear operator $𝑨\in {ℂ}^{m\times m}$ as in

suggests that algorithms that construct powers of $𝑨$ might reveal eigenvalues. This is indeed the case and leads to a class of algorithms of wide applicability in scientific computation. First, observe that taking condition numbers gives

where ${\left|\lambda \right|}_{max}$, ${\left|\lambda \right|}_{min}$ are the eigenvalues of maximum and minimum absolute values. While these express an intrinsic property of the operator $𝑨$, the factor ${\mu }^2\left(𝑿\right)$ is associated with the conditioning of a change of coordinates, and a natural question is whether it is possible to avoid any ill-conditioning associated with a basis set $𝑿$ that is close to linear dependence. The answer to this line of inquiry is given by the following result.

Schur Theorem. For any $𝑨\in {ℂ}^{m\times m}$ there exists $𝑸$ unitary and $𝑻$ upper triangular such that $𝑨=𝑸𝑻{𝑸}^{\ast }$.

Proof. Proceed by induction, starting from an arbitrary eigenvalue $\lambda$ and eigenvector $𝒙$. Let ${𝒖}_1=𝒙/||𝒙||$, the first column vector of a unitary matrix $𝑼=\left[\begin{array}{ll}{𝒖}_1& 𝑽\end{array}\right]$. Then

with $𝑪\in {ℂ}^{(m-1)\times (m-1)}$ that by the inductive hypothesis can be written as $𝑪=𝑾𝑺{𝑾}^{\ast }$, with $𝑾$ unitary, $𝑺$ upper triangular. The matrix

is a product of unitary matrices, hence itself unitary. The computation

then shows that $𝑻$ is indeed triangular.$\square$

The eigenvalues of an upper triangular matrix are simply its diagonal elements, so the Schur factorization is an eigenvalue-revealing factorization.

2.Power iteration for real symmetric matrices

When the operator $𝑨$ expresses some physical phenomenon, the principle of action and reaction implies that $𝑨\in {ℝ}^{m\times m}$ is symmetric, $𝑨={𝑨}^T$ and has real eigenvalues. Componentwise, symmetry of $𝑨=\left[a_{ij}\right]$ implies $a_{ij}=a_{ji}$. Consider $𝑨𝒙=\lambda 𝒙$, and take the adjoint to obtain ${𝒙}^T{𝑨}^T=\bar{\lambda }{𝒙}^T$, or ${𝒙}^T𝑨=\bar{\lambda }{𝒙}^T$ since $𝑨$ is symmetric. Form scalar products ${𝒙}^T𝑨𝒙=\lambda {𝒙}^T𝒙$, ${𝒙}^T{𝑨}^T𝒙=\bar{\lambda }{𝒙}^T𝒙$, and subtract to obtain

since $𝒙\ne 𝟎$, an eigenvector.

Example. Consider a linear array of identical mass-springs. The $i^{\mathrm{th}}$ point mass obeys the dynamics

expressed in matrix form as $\ddot{𝒙}=𝑨𝒙$, with $𝑨$ symmetric.

For a real symmetric matrix the Schur theorem states that

and since a symmetric triangular matrix is diagonal, the Schur factorization is also an eigendecomposition, and the eigenvector matrix $𝑸$ is a basis, $C\left(𝑸\right)={ℝ}^m$.

2.1.The power iteration idea

Assume initially that the eigenvalues are distinct and ordered $\left|{\lambda }_1\right|>\left|{\lambda }_2\right|>\cdots >\left|{\lambda }_m\right|$. Repeated application of $𝑨$ on an arbitrary vector $𝒗=𝑸𝒄\in {ℝ}^m=C\left(𝑸\right)$ is expressed as

a linear combination of the columns of $𝑸$ (eigenvectors of $𝑨$) with coefficients ${𝚲}^n𝒄={\left[\begin{array}{llll}{\lambda }_1^n{c}_1& {\lambda }_2^n{c}_2& \dots & {\lambda }_m^n{c}_m\end{array}\right]}^T$.

∴ 

function genAQΛ(m)
  X=rand(m,m); Q=qr(X).Q; Λ=diagm(2.0 .^ (0:-1:1-m))
  A=Q*Λ*Q'
  return A,Q,Λ
end;

∴ 

∴ 

m=6; A,Q,Λ=genAQΛ(m);

∴ 

n=20; V=rand(m,n+1);

∴ 

for k=1:n
  global V
  V[:,k+1] = A*V[:,k]  
end

∴ 

∴ 

C=Q'*V;

∴ 

clf(); plot(log10.(abs.(C)));

∴ 

xlabel("Iteration"); ylabel("lg|c_i|");

∴ 

The sequence of normalized eigenvector approximants ${𝒗}_n={𝑨}^n𝒗/||{𝑨}^n𝒗||$ is linearly convergent at rate $r=|{\lambda }_2/{\lambda }_1|$.

2.2.Rayleigh quotient

To estimate the eigenvalue revealed by power iteration, formulate the least squares problem

that seeks the best approximation of one power iteration $𝑨𝒗$ as a linear combination of the initial vector $𝒗$. Of course, if $𝒗=𝒒$ is an eigenvector, then the solution would be $c=\lambda$, the associated eigenvalue. The projector onto $C\left(𝒗\right)$ is

that when applied to $𝑨𝒗$ gives the equation

The function $r:{ℝ}^m\to ℝ$,

is known as the Rayleigh quotient which, evaluated for an eigenvector, gives $r\left(𝒒\right)=\lambda$. To determine how well the eigenvalue is approximated, carry out a Taylor series in the vicinity of an eigenvector $𝒒$

where ${\nabla }_{𝒗}r$ is the gradient of $r\left(𝒗\right)$

Compute the gradient through differentiation of the Rayleigh quotient

Noting that ${\nabla }_{𝒗}{v}_i={𝒆}_i$, the $i^{\mathrm{th}}$ column of $𝑰$, the gradient of ${𝒗}^T𝒗$ is

To compute ${\nabla }_{𝒗}({𝒗}^T𝑨𝒗)$, let $𝒖=𝑨𝒗$, and since $𝑨$ is symmetric ${𝒖}^T={𝒗}^T{𝑨}^T={𝒗}^T𝑨$, leading to

Use $u_i={\sum }_{j=1}^m{a}_{ij}{v}_j$ also expressed as $u_j={\sum }_{i=1}^m{a}_{ji}{v}_i$ by swapping indices to obtain

and therefore

Use symmetry of $𝑨$ to write

and substitute above to obtain

Gathering the above results

gives the following gradient of the Rayleigh quotient

When evaluated at $𝒗=𝒒$, obtain ${\nabla }_{𝒗}r\left(𝒒\right)=𝟎$, implying that near an eigenvector the Rayleigh quotient approximation of an eigenvalue is of quadratic accuracy,

2.3.Refining the power iteration idea

Power iteration furnishes the largest eigenvalue. Further eigenvalues can be found by use of the following properties:

• $(\lambda ,𝒒)$ eigenpair of $𝑨\Rightarrow (\lambda -\mu ,𝒒)$ eigenpair of $𝑨-\mu 𝑰$;

• $(\lambda ,𝒒)$ eigenpair of $𝑨\Rightarrow (1/\lambda ,𝒒)$ eigenpair of ${𝑨}^{-1}$.

If $\mu$ is a known initial approximation of the eigenvalue then the inverse power iteration ${𝒗}_n={(𝑨-\mu 𝑰)}^{-1}{𝒗}_{n-1}$, actually implemented as successive solution of linear systems

leads to a sequence of Rayleigh quotients $𝒓\left({𝒗}_n\right)$ that converges quadratically to an eigenvalue close to $\mu$. An important refinement of the idea is to change the shift at each iteration which leads to cubic order of convergence

Power iteration can be applied simultaneously to multiple directions at once