next up previous
Next: Inverse Power Method Up: NUMERICAL TECHNIQUES FOR EIGENVALUES Previous: NUMERICAL TECHNIQUES FOR EIGENVALUES

The Power Method

Used to estimate the dominant eigenvalue and its corresponding eigenvector.

There's two variants: the case when $A$ ($n\times n$ matrix) has

$
\left\{\begin{array}{l}
n \mbox { linearly independent eigenvectors }\\
\mbox {there's a single dominant eigenvalue.}
\end{array}\right.$

this one's easy.

The second variant is not discussed here and considers the cases not covered by above assumptions. This one's complicated. See reference to Demmel et al.

The first case: first, we order the eigenvalues

(i)
$\vert\lambda_1\vert>\vert\lambda_2\vert\ge \vert\lambda_3\vert\ge \cdots \ge\vert\lambda_n\vert$
(ii)
$\mbox { There's a basis } \left\{u^{(1)}, u^{(2)}, \cdots,
u^{(n)}\right\}\\
\mbox { for } {\mathbb{C}}^n \mbox { such that } $

(80) \begin{displaymath}
Au^{(i)}=\lambda_iu^{(i)}\quad (1\le i\le n)
\end{displaymath}

Let $x^{(0)}$ be any element of ${\mathbb{C}}^n$ such that when $x^{(0)}$ is expressed as a linear combination of the basis $u^{(1)}, u^{(2)},
\ldots, u^{(n)}4$ the coefficient of $u^{(1)}$ is not 0. Thus,
(81) \begin{displaymath}
x^{(0)}=a_1 u^{(1)}+a_2u^{(2)}+\cdots u_nu^{(n)}\quad a_1\ne 0
\end{displaymath}

We form

\begin{eqnarray*}
x^{(1)}&=&Ax^{(0)}, \\
x^{(2)}&=&Ax^{(1)}, \\
& & \cdots \\
x^{(k)}&=& Ax^{(k-1)}
\end{eqnarray*}



so that
(82) \begin{displaymath}
x^{(k)}=A^k x^{(0)}
\end{displaymath}

In what follows we'll absorb all the coefficients $a_j$ in the vectors $u^{(j)}$ that they multiply. Hence, write (82) as
(83) \begin{displaymath}
x^{(0)}=u^{(1)}+u^{(2)}+\ldots + u^{(n)}
\end{displaymath}

by (84) we can write (83) as

\begin{displaymath}
x^{(k)}=A^{(k)}u^{(1)}+A^{(k)}u^{(2)}+\ldots A^k u^{(n)}
\end{displaymath}

using (81)

\begin{displaymath}
x^{(k)}=\lambda^k_1 u^{(1)}+\lambda^k_2 u^{(2)}+\ldots
\lambda^k_n u^{(n)}
\end{displaymath}

rewrite this
(84) \begin{displaymath}
x^{(k)}=\lambda^k_1\left[u^{(1)}+\left(\frac{\lambda_2}{\la...
...ots \left(\frac{\lambda_n}{\lambda_1}\right)^k
u^{(n)}
\right]
\end{displaymath}

Since $\vert\lambda_1\vert>\vert\lambda_j\vert$ for $2\le j\le n$ we see that the coefficients $(\lambda_j/\lambda_1)^k$ tend to 0 and the quantity in brackets converges to $u^{(1)}$ as $k \rightarrow
\infty$.

Let's write (85) as

\begin{displaymath}
x^{(k)}=\lambda^k_1[u^{(1)}+\varepsilon ^{(k)}]
\end{displaymath}

where $\varepsilon ^{(k)}\rightarrow 0$ as $k \rightarrow
\infty$. In order to take ratios let $\theta$ be any linear functional on ${\mathbb{C}}^n$ for which $\theta (u^{(1)})\ne 0$. Note: $\theta (\alpha x+\beta
y)=\alpha\theta(x)+\beta \theta (y)$, since linear, where $\alpha,
\beta$ are scalars and $x$ and $y$ vectors. Then

\begin{displaymath}
\theta(x^{(k)})=\lambda^k_1[\theta(u^{(1)})+\theta (\varepsilon ^{(k)})]
\end{displaymath}

% latex2html id marker 16308
$\therefore$ the following ratios converge to $\lambda_1$ as $k \rightarrow
\infty$:

\begin{displaymath}
r_k\equiv\frac{\theta(x^{(k+1)})}{\theta(x^{(k)})}=\lambda_...
...(1)})+
\theta(\varepsilon ^{(k)})}\right]\rightarrow \lambda_1
\end{displaymath}

$\Box$

Remark: Since direction of $x^{(k)}$ aligns more and more with $u^{(1)}$ as $k \rightarrow
\infty$ the method also yields $u^{(1)}$ the eigenvector.

What to chose for $\theta$? It can simply evaluate the $j^{th}$ component of any given vector. In practice, the algorithm should normalize the $x^{(k)}$ otherwise they may converge to $0$ or become unbounded.

The algorithm goes like this:

input: $n, A, x, m$ % here $x$ is an initial guess

output: $0, x$

\begin{eqnarray*}
\mbox {do } k &=& 1, m\\
y &=& Ax\\
r &=& \theta (y)/\the...
...\vert\vert y\vert\vert\quad\% \mbox { use sup-norm for example.}
\end{eqnarray*}



end

Relative error: Can show that, since $\displaystyle \lim_{k\to\infty}r_k=\lambda_1,$

\begin{displaymath}
\frac{r_k-\lambda_1}{\lambda_1}=\left(\frac{\lambda_2}{\lambda_1}
\right)^k c_k
\end{displaymath}

where the numbers $c_k$ form a bounded sequence.

Can also show:

\begin{eqnarray*}
r_{k+1}-\lambda_1=(c+\delta k)(r_k-\lambda_1)\\
\vert c\vert<1 \mbox { and } \lim_{k\to \infty} \delta_k=0
\end{eqnarray*}



this implies $r_k$ converges linearly. Using this knowledge one can accelerate the procedure:

Theorem (Aitken Acceleration): Let $\{r_n\}$ be a sequence of numbers that converges to a limit $r$, then the sequence

\begin{displaymath}
s_n=\frac{r_n r_{n+2}-r^2_{n+1}}{r_{n+2}-2r_{n+1}+r_n}\quad n\ge
0
\end{displaymath}

converges to $r$ faster if $r_{n+1}-r=(c+\delta_n)(r_n-r)$ with $\vert c\vert<1$ and $\displaystyle \lim_{n\to \infty}\delta_n=0$. Indeed, $(s_n-r)/(r_n-r)\to 0$ as $n\to \infty $.

$\Box$


next up previous
Next: Inverse Power Method Up: NUMERICAL TECHNIQUES FOR EIGENVALUES Previous: NUMERICAL TECHNIQUES FOR EIGENVALUES
Juan Restrepo 2003-04-12