Power iteration

File:Animation of the Power Iteration Algorithm.gif

The power iteration algorithm starts with a vector $b\_0$, which may be an approximation to the dominant eigenvector or a random vector. The method is described by the recurrence relation

:$b\_\{k+1\}\; =\; \backslash frac\{Ab\_k\}\{\backslash |Ab\_k\backslash |\}$

So, at every iteration, the vector $b\_k$ is multiplied by the matrix $A$ and normalized.

If we assume $A$ has an eigenvalue that is strictly greater in magnitude than its other eigenvalues and the starting vector $b\_0$ has a nonzero component in the direction of an eigenvector associated with the dominant eigenvalue, then a subsequence $\backslash left(\; b\_\{k\}\; \backslash right)$ converges to an eigenvector associated with the dominant eigenvalue.

Without the two assumptions above, the sequence $\backslash left(\; b\_\{k\}\; \backslash right)$ does not necessarily converge. In this sequence,

: $b\_k\; =\; e^\{i\; \backslash phi\_k\}\; v\_1\; +\; r\_k$,

where $v\_1$ is an eigenvector associated with the dominant eigenvalue, and $\backslash |\; r\_\{k\}\; \backslash |\; \backslash rightarrow\; 0$. The presence of the term $e^\{i\; \backslash phi\_\{k\}\}$ implies that $\backslash left(\; b\_\{k\}\; \backslash right)$ does not converge unless $e^\{i\; \backslash phi\_\{k\}\}\; =\; 1$. Under the two assumptions listed above, the sequence $\backslash left(\; \backslash mu\_\{k\}\; \backslash right)$ defined by

: $\backslash mu\_\{k\}\; =\; \backslash frac\{b\_\{k\}^\{*\}Ab\_\{k\}\}\{b\_\{k\}^\{*\}b\_\{k\}\}$

converges to the dominant eigenvalue (with Rayleigh quotient).{{Clarify|what's the relation between this sequence and the power method|date=October 2016}}

One may compute this with the following algorithm (shown in Python with NumPy):

1. !/usr/bin/env python3

import numpy as np

def power_iteration(A, num_iterations: int):

# Ideally choose a random vector

# To decrease the chance that our vector

# Is orthogonal to the eigenvector

b_k = np.random.rand(A.shape[1])

for _ in range(num_iterations):

# calculate the matrix-by-vector product Ab

b_k1 = np.dot(A, b_k)

# calculate the norm

b_k1_norm = np.linalg.norm(b_k1)

# re normalize the vector

b_k = b_k1 / b_k1_norm

return b_k

power_iteration(np.array(0.5, 0.5], [0.2, 0.8), 10)

The vector $b\_k$ converges to an associated eigenvector. Ideally, one should use the Rayleigh quotient in order to get the associated eigenvalue.

This algorithm is used to calculate the Google PageRank.

The method can also be used to calculate the spectral radius (the eigenvalue with the largest magnitude, for a square matrix) by computing the Rayleigh quotient

:$\backslash rho(A)\; =\; \backslash max\; \backslash left\; \backslash \{\; |\backslash lambda\_1|,\; \backslash dotsc,\; |\backslash lambda\_n|\; \backslash right\; \backslash \}\; =\; \backslash frac\{b\_k^\backslash top\; A\; b\_k\}\{b\_k^\backslash top\; b\_k\}.$

Let $A$ be decomposed into its Jordan canonical form: $A=VJV^\{-1\}$, where the first column of $V$ is an eigenvector of $A$ corresponding to the dominant eigenvalue $\backslash lambda\_\{1\}$. Since generically, the dominant eigenvalue of $A$ is unique, the first Jordan block of $J$ is the $1\; \backslash times\; 1$ matrix $[\backslash lambda\_1],$ where $\backslash lambda\_\{1\}$ is the largest eigenvalue of A in magnitude. The starting vector $b\_0$ can be written as a linear combination of the columns of V:

:$b\_\{0\}\; =\; c\_\{1\}v\_\{1\}\; +\; c\_\{2\}v\_\{2\}\; +\; \backslash cdots\; +\; c\_\{n\}v\_\{n\}.$

By assumption, $b\_\{0\}$ has a nonzero component in the direction of the dominant eigenvalue, so $c\_\{1\}\; \backslash ne\; 0$.

The computationally useful recurrence relation for $b\_\{k+1\}$ can be rewritten as:

:$b\_\{k+1\}=\backslash frac\{Ab\_\{k\}\}\{\backslash |Ab\_\{k\}\backslash |\}=\backslash frac\{A^\{k+1\}b\_\{0\}\}\{\backslash |A^\{k+1\}b\_\{0\}\backslash |\},$

where the expression: $\backslash frac\{A^\{k+1\}b\_\{0\}\}\{\backslash |A^\{k+1\}b\_\{0\}\backslash |\}$ is more amenable to the following analysis.

:$\backslash begin\{align\}$

b_k &= \frac{A^{k}b_{0}}{\| A^{k} b_{0} \|} \\

&= \frac{\left( VJV^{-1} \right)^{k} b_{0}}{\|\left( VJV^{-1} \right)^{k}b_{0}\|} \\

&= \frac{ VJ^{k}V^{-1} b_{0}}{\| V J^{k} V^{-1} b_{0}\|} \\

&= \frac{ VJ^{k}V^{-1} \left( c_{1}v_{1} + c_{2}v_{2} + \cdots + c_{n}v_{n} \right)}{\| V J^{k} V^{-1} \left( c_{1}v_{1} + c_{2}v_{2} + \cdots + c_{n}v_{n} \right)\|} \\

&= \frac{ VJ^{k}\left( c_{1}e_{1} + c_{2}e_{2} + \cdots + c_{n}e_{n} \right)}{\| V J^{k} \left( c_{1}e_{1} + c_{2}e_{2} + \cdots + c_{n}e_{n} \right) \|} \\

&= \left( \frac{\lambda_{1}}

\end{align}

The expression above simplifies as $k\; \backslash to\; \backslash infty$

:$\backslash left(\; \backslash frac\{1\}\{\backslash lambda\_\{1\}\}\; J\; \backslash right)^\{k\}\; =$

\begin{bmatrix}

[1] & & & & \\

& \left( \frac{1}{\lambda_{1}} J_{2} \right)^{k}& & & \\

& & \ddots & \\

& & & \left( \frac{1}{\lambda_{1}} J_{m} \right)^{k} \\

\end{bmatrix}

\rightarrow

\begin{bmatrix}

1 & & & & \\

& 0 & & & \\

& & \ddots & \\

& & & 0 \\

\end{bmatrix} \quad \text{as} \quad k \to \infty.

The limit follows from the fact that the eigenvalue of $\backslash frac\{1\}\{\backslash lambda\_\{1\}\}\; J\_\{i\}$ is less than 1 in magnitude, so

:$\backslash left(\; \backslash frac\{1\}\{\backslash lambda\_\{1\}\}\; J\_\{i\}\; \backslash right)^\{k\}\; \backslash to\; 0\; \backslash quad\; \backslash text\{as\}\; \backslash quad\; k\; \backslash to\; \backslash infty.$

It follows that:

:$\backslash frac\{1\}\{c\_\{1\}\}\; V\; \backslash left(\; \backslash frac\{1\}\{\backslash lambda\_1\}\; J\; \backslash right)^\{k\}\; \backslash left(\; c\_\{2\}e\_\{2\}\; +\; \backslash cdots\; +\; c\_\{n\}e\_\{n\}\; \backslash right)\; \backslash to\; 0\; \backslash quad\; \backslash text\{as\}\; \backslash quad\; k\; \backslash to\; \backslash infty$

Using this fact, $b\_\{k\}$ can be written in a form that emphasizes its relationship with $v\_\{1\}$ when k is large:

:$\backslash begin\{align\}$

b_k &= \left( \frac{\lambda_{1}}

&= e^{i \phi_{k}} \frac{c_{1}}

\end{align}

where $e^\{i\; \backslash phi\_\{k\}\}\; =\; \backslash left(\; \backslash lambda\_\{1\}\; /\; |\backslash lambda\_\{1\}|\; \backslash right)^\{k\}$ and $\backslash |\; r\_\{k\}\; \backslash |\; \backslash to\; 0$ as $k\; \backslash to\; \backslash infty$

The sequence $\backslash left(\; b\_\{k\}\; \backslash right)$ is bounded, so it contains a convergent subsequence. Note that the eigenvector corresponding to the dominant eigenvalue is only unique up to a scalar, so although the sequence $\backslash left(b\_\{k\}\backslash right)$ may not converge,

$b\_\{k\}$ is nearly an eigenvector of A for large k.

Alternatively, if A is diagonalizable, then the following proof yields the same result

Let λ₁, λ₂, ..., λ_m be the m eigenvalues (counted with multiplicity) of A and let v₁, v₂, ..., v_m be the corresponding eigenvectors. Suppose that $\backslash lambda\_1$ is the dominant eigenvalue, so that $|\backslash lambda\_1|\; >\; |\backslash lambda\_j|$ for $j>1$.

The initial vector $b\_0$ can be written:

:$b\_0\; =\; c\_\{1\}v\_\{1\}\; +\; c\_\{2\}v\_\{2\}\; +\; \backslash cdots\; +\; c\_\{m\}v\_\{m\}.$

If $b\_0$ is chosen randomly (with uniform probability), then c₁ ≠ 0 with probability 1. Now,

:$\backslash begin\{align\}$

A^{k}b_0 &= c_{1}A^{k}v_{1} + c_{2}A^{k}v_{2} + \cdots + c_{m}A^{k}v_{m} \\

&= c_{1}\lambda_{1}^{k}v_{1} + c_{2}\lambda_{2}^{k}v_{2} + \cdots + c_{m}\lambda_{m}^{k}v_{m} \\

&= c_{1}\lambda_{1}^{k} \left( v_{1} + \frac{c_{2}}{c_{1}}\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k}v_{2} + \cdots + \frac{c_{m}}{c_{1}}\left(\frac{\lambda_{m}}{\lambda_{1}}\right)^{k}v_{m}\right) \\

&\to c_{1}\lambda_{1}^{k} v_1 && \left |\frac{\lambda_j}{\lambda_1} \right | < 1 \text{ for } j>1

\end{align}

On the other hand:

:$b\_k\; =\; \backslash frac\{A^k\; b\_0\}\{\backslash |A^kb\_0\backslash |\}.$

Therefore, $b\_k$ converges to (a multiple of) the eigenvector $v\_1$. The convergence is geometric, with ratio

:$\backslash left|\; \backslash frac\{\backslash lambda\_2\}\{\backslash lambda\_1\}\; \backslash right|,$

where $\backslash lambda\_2$ denotes the second dominant eigenvalue. Thus, the method converges slowly if there is an eigenvalue close in magnitude to the dominant eigenvalue.

Although the power iteration method approximates only one eigenvalue of a matrix, it remains useful for certain computational problems. For instance, Google uses it to calculate the PageRank of documents in their search engine,{{cite news |author=Ipsen, Ilse, and Rebecca M. Wills |date=5–8 May 2005 |title=7th IMACS International Symposium on Iterative Methods in Scientific Computing |location=Fields Institute, Toronto, Canada |url=https://ipsen.math.ncsu.edu/ps/slides_imacs.pdf}} and Twitter uses it to show users recommendations of whom to follow.Pankaj Gupta, Ashish Goel, Jimmy Lin, Aneesh Sharma, Dong Wang, and Reza Bosagh Zadeh [http://dl.acm.org/citation.cfm?id=2488433 WTF: The who-to-follow system at Twitter], Proceedings of the 22nd international conference on World Wide Web The power iteration method is especially suitable for sparse matrices, such as the web matrix, or as the matrix-free method that does not require storing the coefficient matrix $A$ explicitly, but can instead access a function evaluating matrix-vector products $Ax$. For non-symmetric matrices that are well-conditioned the power iteration method can outperform more complex Arnoldi iteration. For symmetric matrices, the power iteration method is rarely used, since its convergence speed can be easily increased without sacrificing the small cost per iteration; see, e.g., Lanczos iteration and LOBPCG.

Some of the more advanced eigenvalue algorithms can be understood as variations of the power iteration. For instance, the inverse iteration method applies power iteration to the matrix $A^\{-1\}$. Other algorithms look at the whole subspace generated by the vectors $b\_k$. This subspace is known as the Krylov subspace. It can be computed by Arnoldi iteration or Lanczos iteration.

Gram iteration{{Citation

| last1=Delattre

| first1=B.

| last2=Barthélemy

| first2=Q.

| last3=Araujo

| first3=A.

| last4=Allauzen

| first4=A.

| title=Efficient Bound of Lipschitz Constant for Convolutional Layers by Gram Iteration

| journal=Proceedings of the 40th International Conference on Machine Learning

| year=2023

| pages=7513–7532

| url=https://proceedings.mlr.press/v202/delattre23a.html

}} is a super-linear and deterministic method to compute the largest eigenpair.

• Rayleigh quotient iteration
• Inverse iteration

{{Reflist}}

{{Numerical linear algebra}}

{{DEFAULTSORT:Power Iteration}}

Category:Numerical linear algebra

Category:Articles with example Python (programming language) code