The Power Methods for Eigenvalue problems

In some applications, only the eigenvalue $\lambda_{max}$ or $\lambda_{min}$ of maximum or minimum absolution values (if real) or moduli (if commplex) is of interes (e.g., projecting data points in a high-dimensional space into a single dimension in the pricipal component direction). In such cases, we may be able to use the methods of power iteration to find $\lambda_{max}$ and $\lambda_{min}$ together with their corresponding eigenvectors, based on the assumption that ${\bf A}$ is a full-rank square matrix. In the following, we assume the eigenvalues are arranged in descending order based on their absolution values or modulus:

$$\vert\lambda_1\vert=\vert\lambda_{max}\vert>\vert\lambda_2\vert\ge\cdots\ge\vert\lambda_{N-1}\vert >\vert\lambda_N\vert=\vert\lambda_{min}\vert$$ (39)

• The basic power algorithm:

  We express the eigenvector matrix of ${\bf A}$ and its inverse as ${\bf V}=[{\bf v}_1,\cdots,{\bf v}_N]$ and ${\bf V}^{-1}=[{\bf v}'_1,\cdots,{\bf v}'_N]^T$, respectively, and write the eigenequation ${\bf A}{\bf V}={\bf V}{\bf\Lambda}$ as

  $${\bf A}={\bf V}{\bf\Lambda}{\bf V}^{-1} =[{\bf v}_1,\cdots,{\bf v... ...ime T} \end{array}\right] =\sum_{i=1}^N \lambda_i{\bf v}_i{\bf v}_i^{\prime T }$$ (40)

  As ${\bf A}^n{\bf V}={\bf V}{\bf\Lambda}^n\;(n=-1,\,0,\,1,\,2,\cdots)$ (see here), we also have

  $${\bf A}^n={\bf V}{\bf\Lambda}^n{\bf V}^{-1} =\sum_{i=1}^N \lambda_i^n{\bf v}_i{\bf v}_i^{\prime T}$$ (41)

  The power algorithm is a simple iteration based on arbitrary initial guess ${\bf u}_0$ of the eigenvector:

  $${\bf u}_n={\bf Au}_{n-1}={\bf A}^2{\bf u}_{n-2}=\cdots={\bf A}^n{... ... =\left( \sum_{i=1}^N \lambda_i^n{\bf v}_i{\bf v}_i^{\prime T} \right){\bf u}_0$$ (42)

  Pre-multiplying ${\bf u}^T_n$ on both sides we get

  $${\bf u}^T_n{\bf u}_n =\sum_{i=1}^N \lambda_i^n{\bf u}^T_n({\bf v}... ...=\sum_{i=1}^N \lambda_i^n({\bf u}^T_n{\bf v}_i)({\bf v}_i^{\prime T} {\bf u}_0)$$ (43)

  Dividing both sides by $\lambda_1^n$ and taking the limit $n\rightarrow\infty$, we get:

  $$\lim_{n\rightarrow\infty}\frac{{\bf u}^T_n{\bf u}_n}{\lambda _1^n... ...i^{\prime T} {\bf u}_0) =({\bf u}^T_n{\bf v}_1)({\bf v}_1^{\prime T} {\bf u}_0)$$ (44)

  The last equality is due to the fact that $\vert\lambda_i\vert/\vert\lambda_1\vert<1$ for all $i>1$. We also get the same result for $n+1$ iterations:

  $$\lim_{n\rightarrow\infty}\frac{{\bf u}^T_n{\bf u}_{n+1}}{\lambda _1^{n+1}} =({\bf u}^T_n{\bf v}_1)({\bf v}_1^{\prime T} {\bf u}_0)$$ (45)

  As the two equations above have the same right-hand side, we equate their left-hand sides to get:

  $$\lim_{n\rightarrow\infty}\frac{{\bf u}^T_n{\bf u}_{n+1}}{\lambda ... ...\lambda_1\lim_{n\rightarrow\infty}\frac{{\bf u}^T_n{\bf u}_n}{\lambda _1^{n+1}}$$ (46)

  solving for $\lambda_1$, we get the eigenvalue with the greatest absolute value:

  $$\lambda_1 =\lim_{n\rightarrow\infty}\frac{{\bf u}^T_n{\bf u}_{n+1... ...im_{n\rightarrow\infty}\frac{{\bf u}_n^T{\bf A}{\bf u}_n}{{\bf u}_n^T{\bf u}_n}$$ (47)

  This is the Rayleigh quotient.

  To find the eigenvector ${\bf v}_1$ corresponding to $\lambda_1$, we express the arbitrary initial vector as a linear combination ${\bf u}_0=\sum_{i=1}^Nc_i{\bf v}_i$ of the $N$ independent eigenvectors of ${\bf A}$, and get
  

  $${\bf u}_n = {\bf A}^n{\bf u}_0={\bf A}^n\left(\sum_{i=1}^Nc_i{\bf v}_i\right) =\sum_{i=1}^Nc_i{\bf A}^n{\bf v}_i$$

  $$= \sum_{i=1}^Nc_i\lambda_i^n{\bf v}_i =\lambda_1^n\left[c_1{\bf v}_1 +\sum_{i=2}^Nc_i\left(\frac{\lambda_i}{\lambda_1}\right)^n{\bf v}_i\right]$$ (48)

  
  At the limit $n\rightarrow\infty$, $(\lambda_i/\lambda_1)^n=0$ for all $i=2,\cdots,N$, and ${\bf u}_n$ approaches the first eigenvector ${\bf v}_1$ with some scaling coefficient:

  $$\lim_{n\rightarrow\infty}{\bf u}_n=\lambda_1^nc_1\;{\bf v}_1\propto {\bf v}_1$$ (49)

  By normalizing ${\bf u}_n$, we get the normalized eigenvector ${\bf v}_1$ corresponding to $\lambda_1$. One obvious requirement is that $c_1\ne 0$, i.e., the initial guess ${\bf u}_0$ must have a non-zero component in the direction of ${\bf v}_1$. Also, to avoid the possible overflow or underflow due to the iteration ${\bf u}_n={\bf A}{\bf u}_{n-1}$, we normalize ${\bf u}_n$ every time a new result is generated. In other words, at the end of the iteration, the normalized ${\bf u}_n$ is already the normalized eigenvector corresponding to $\lambda_1$.

  In summary, here are the steps of the power method:

  • Choose an artibrary ${\bf u}_0$ with $\vert\vert{\bf u}_0\vert\vert=1$;
  • Iterate for $n=1,\;2,\cdots$ until convergence, i.e., the difference between two consecutive eigenvalues $\vert\lambda_n-\lambda_{n-1}\vert$ is smaller than a preset error:

    ${\bf w}={\bf A}{\bf u}_{n-1},\;\;\;\;\;\;\; {\bf u}_n={\bf w}/\vert\vert{\bf w}\vert\vert,\;\;\;\; \lambda_n={\bf u}_n{\bf A}{\bf u}_n$

• Inverse iteration:

  The power method can be modified to find the smallest eigenvalue and the corresponding eigenvector, based on the fact that the eigenvalues of ${\bf A}^{-1}$ are the reciprocals of those of ${\bf A}$, i.e., if ${\bf A}{\bf V}={\bf V}{\bf\Lambda}$, then ${\bf A}^{-1}{\bf V}={\bf V}{\bf\Lambda}^{-1}$. To find the smallest eigenvalue $\lambda_N$ of ${\bf A}$, we can use the power method to find the largest eigenvalue of ${\bf A}^{-1}$ by iteration ${\bf u}_n={\bf A}^{-1}{\bf u}_{n-1}$, which is equivalent to solving ${\bf A}{\bf u}_n={\bf u}_{n-1}$ for ${\bf u}_n$. Here is the iteration step:

  • Choose an arbitrary ${\bf u}_0$ with $\vert\vert{\bf u}_0\vert\vert=1$;
  • Iterate for $n=1,\;2,\cdots$ until convergence:

    ${\bf w}={\bf A}^{-1}{\bf u}_{n-1},\;\;\;\;\;\;\; {\bf u}_n={\bf w}/\vert\vert{\bf w}\vert\vert,\;\;\;\;\;\;\; \lambda_n={\bf u}_n{\bf A}{\bf u}_n$

• Inverse iteration with shift

  The power method can be further modified to find the eigenvalues between the largest and the smallest ones. Let $\lambda$ and ${\bf v}$ be one of the eigenvalues and the corresponding eigenvector of ${\bf A}$, i.e., ${\bf Av}=\lambda{\bf v}$, then we have

  $$({\bf A}-\mu{\bf I}){\bf v}={\bf Av}-\mu{\bf v} =\lambda{\bf v}-\mu{\bf v}=(\lambda-\mu){\bf v}$$ (50)

  i.e., $\lambda-\mu$ is an eigenvalue of ${\bf A}-\mu{\bf I}$, and $1/(\lambda-\mu)$ is an eigenvalue of $({\bf A}-\mu{\bf I})^{-1}$. If $\mu$, called a shift, is closest to one of the eigenvalues $\lambda_i$ of ${\bf A}$, i.e., $\vert\lambda_i-\mu\vert\le\vert\lambda_j-\mu\vert$ for all $j=1,\cdots,N$, then $1/(\lambda_i-\mu)$ is the maximum eigenvalue of ${\bf A}-\mu{\bf I}$. Now we can find the eigenvalue $\lambda_i$ of ${\bf A}$ closest to $\mu$ and the corresponding eigenvector:
  • Choose an artitrary ${\bf u}_0$ with $\vert\vert{\bf u}_0\vert\vert=1$, and get an initial guess $\mu$ of a certain eigenvalue $\lambda_i$ (closest to $\lambda_i$);

  • Iterate for $n=1,\;2,\cdots$ until convergence:

    ${\bf w}=({\bf A}-\mu{\bf I})^{-1}{\bf u}_{n-1},\;\;\;\;\;\;\; {\bf u}_n={\bf w}/\vert\vert{\bf w}\vert\vert,\;\;\;\;\;\;\; \lambda_n={\bf u}_n{\bf A}{\bf u}_n$

We make a few observations about the power method.

• The rate of the convergence depends on the ratio $\vert\lambda_2/\lambda_1\vert<1$, the smaller the ratio, the faster the convergence.
• If a real matrix ${\bf A}$ has a complex conjugate eigenvalue pair $\lambda_{1,2}=a\pm jb$ of equal modulus $\vert\lambda_1\vert=\vert\lambda_2\vert=\sqrt{a^2+b^2}$, which happens to be the maximum (or minimum) among all eigenvalues, then $\vert\lambda_2/\lambda_1\vert=1$ and the power method (or the inverse iteration method) will not converge.
• If an eigenvalue is repeated $k$ times, there exist $k$ corresponding eigenvectors. However, the power method can find only one eigenvector, which is a linear combination of the $k$ eigenvectors.

For example, if the eigenvalues of a $4\times 4$ real matrix are $\{-9,\;1,\;3\pm j4\}$, then the power iteration methods above will find both $\lambda_{max}=-9$ and $\lambda_{min}=1$. But if the eigenvalues are $\{-2,\;1,\;3\pm j4\}$, then the power method for $\lambda_{max}$ will not converge because the modulus of the complex conjugate pair $\vert 3\pm j4\vert=5$ is the maximum of all eigenvalues. Similarly if the eigenvalues are $\{-9,\;6,\;3\pm j4\}$, then the inverse iteration method for $\lambda_{min}$ will not converge because $\vert 3\pm j4\vert=5$ is the minimum of all eigenvalues.

Example:

The eigenvalues of the following matrix

$${\bf A}=\left[\begin{array}{rrrrr} 3.7420 & -0.4750 & 2.5242 & -... ... -1.1276\\ 1.6295 & -0.6236 & 2.9369 & -2.0581 & -0.3322 \end{array}\right]$$

are $3,\;3,\;3,\;2,\;1$, and the corresponding eigenvectors are:

$${\bf V}=[{\bf v}_1,\cdots,{\bf v}_5] =\left[\begin{array}{rrrrr}... ...& -0.5262\\ -0.5063 & 0.6923 & 0.5607 & -0.5165 & 0.3986 \end{array}\right]$$

The maximum eigenvalue found by the power method is $\lambda_1=3$, and the corresponding eigenvector ${\bf u}_1$ is a linear combination of the three eigenvectors of this repeated eigenvalue:

$${\bf u}_1 =\left[\begin{array}{r}-0.4898\\ 0.1224\\ -0.5665\\ -0... ...\\ -0.5858\end{array}\right] =0.5923{\bf v}_1-0.4796{\bf v}_2-0.0275{\bf v}_3$$

The Matlab code for the power method is listed below:

% The power algorithm to find the largest and smallest eigenvalues and
% eigenvectors of a given matrix A. So, given an estimate mu of any of the
% eigenvalues, find the eigenvalue closest to mu and its eigenvector.

rng(2)
clear all
n=4;
a=3; b=4; 

A=eye(n);
A(1,1)=a; A(2,2)=a;
A(1,2)=b; A(2,1)=-b;
A(3,3)=8; A(4,4)=-2;
A
[v d]=eig(A)

V=rand(n);
A=V*A*inv(V);
A
[v d]=eig(A)
diag(d)'
abs(diag(d))
pause

u=rand(n,1);   u=u/norm(u);
u1=u;
u2=u;
u3=u;
mu=2.9+j*4.1
done=0;
k=0;
s1=1; s2=1; s3=1;
while ~done & k<99
    k=k+1;
    v1=u1;  v2=u2;  v3=u3;  
    r1=s1;  r2=s2;  r3=s3;
    w=A*u1;  u1=w/norm(w);  
    ev1=norm(v1-u1);
    s1=u1'*A*u1;  es1=abs(s1-r1);
    w=linsolve(A,u2); w=inv(A)*u2;  u2=w/norm(w);  
    ev2=norm(v2-u2);
    s2=u2'*A*u2;  es2=abs(s2-r2);
    w=linsolve(A-mu*eye(n),u3);  w=inv(A-mu*eye(n))*u3;   u3=w/norm(w); 
    ev3=norm(v3-u3);
    s3=u3'*A*u3;  es3=abs(s3-r3);
    fprintf('%d  (%e %e %e),  (%e %e %e)\n',k,ev1,ev2,ev3,es1,es2,es3)
%    if ev1+ev2+ev3<10^(-6)
    if es1+es2+es3 < 10^(-6)
        done=1;
    end
end
k
u1';  u1'*A*u1
u2';  u2'*A*u2
u3';  u3'*A*u3