# Power Method

## Basics

• A nonzero vector $\vec{x}$ is an eigenvector of a square matrix $\mathbf{A}$ if there exists a scalar $\lambda$, called an eigenvalue, such that $\mathbf{A}\vec{x} = \lambda\vec{x}$.
• Non-zero vectors in the eigenspace of the matrix $\mathbf{A}$ for the eigenvalue $\lambda$ are eigenvectors of $\mathbf{A}$

## Properties

• The power method is used to find a dominant (largest) eigenvalue (one with the largest absolute value), if one exists, and a corresponding eigenvector.
• It does not always guaranteed to converge if the given matrix is non-diagonalizable.
• Inverse power method (if convergent) calculates the eigenvalue with smallest absolute value.
• It is effective for a very large sparse matrix with appropriate implementation.

## Facts

• This method is known as Von Mises iteration.
• Especially suitable for sparse matrices
• It is used to calculate the Google page rank.

## Algorithm

1. To apply this method to a square matrix $\mathbf{A}$, start with an initial guess $\vec{u}_{0}$ for the eigenvector of the dominant eigenvalue
2. Multiply the most recently obtained vector on the left by $\mathbf{A}$ i.e. for $i \ge 1$, $(\vec{u}_{i})_{\text{un-normalize}} = \mathbf{A} \vec{u}_{i - 1}$

3. Normalize the result i.e. $\vec{u}_{i} = \frac{(\vec{u}_{i})_{\text{un-normalize}}}{|| \mathbf{A} \vec{u}_{i-1}||} = \frac{ \mathbf{A} \vec{u}_{i - 1}}{|| \mathbf{A} \vec{u}_{i-1}||}$

4. Repeat until consecutive vectors $\vec{u}_{i}$ are either identical or opposite.
5. If $\vec{u}_{k}$ denotes the last vector calculated in this process, then $\vec{u}_{k}$ is an approximate eigenvector of $\mathbf{A}$, and $|| \mathbf{A}\vec{u}_{k} ||$ is the absolute value of the dominant eigen vector for $\mathbf{A}$.

6. If $\vec{u}_{k}$ is an eigenvector of $\mathbf{A}$ then, its eigenvalue is given by Rayleigh quotient i.e. $\lambda = \frac{(\mathbf{A} \vec{u}_{k}). (\vec{u}_{k})^{\text{T}}}{\vec{u}_{k}. (\vec{u}_{k})^{\text{T}} }$

It is easy to prove point (6) as $\frac{(\mathbf{A} \vec{u}_{k}). (\vec{u}_{k})^{\text{T}}}{\vec{u}_{k}. (\vec{u}_{k})^{\text{T}} } = \frac{(\lambda \vec{u}_{k}). (\vec{u}_{k})^{\text{T}}}{\vec{u}_{k}. (\vec{u}_{k})^{\text{T}} } = \lambda$.

### Example

• Given a matrix $\mathbf{A}$ defined as

$$\mathbf{A} = \begin{pmatrix} 1 & 3 & 0 \\ 2 & 5 & 1 \\ -1 & 2 & 3 \end{pmatrix}.$$

Now, we need find the largest eigen value.
• Initial guess say

$$\vec{u}_{0} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ \end{pmatrix}$$

• Let $k = 1$,

$$\vec{u}_{1} = \mathbf{A} \vec{u}_{0} = \begin{pmatrix} 1 & 3 & 0 \\ 2 & 5 & 1 \\ -1 & 2 & 3 \end{pmatrix}. \begin{pmatrix} 1 \\ 1 \\ 1 \\ \end{pmatrix} = \begin{pmatrix} 4 \\ 8 \\ 4 \end{pmatrix}$$

So, we normalize $\vec{u}_{1}$ as

$$\vec{u}_{1} = \frac{\vec{u}_{1}}{|| \vec{u}_{1} ||} = \begin{pmatrix} 0.4082 \\ 0.8165 \\ 0.4082 \end{pmatrix}$$

where $|| \vec{u}_{1} || = \sqrt{4^{2} + 8^{2} + 4^{2}}$ i.e. an Euclidean norm.
• Similarly, $k = 2$,

$$\vec{u}_{2} = \mathbf{A} \vec{u}_{1} = \begin{pmatrix} 2.8557 \\ 5.3072 \\ 2.4495 \end{pmatrix}$$

So,

$$\vec{u}_{2} = \frac{\vec{u}_{2}}{|| \vec{u}_{2} ||} = \begin{pmatrix} 0.4392\\ 0.8157 \\ 0.3765 \end{pmatrix}$$

• For $k = 3$,

$$\vec{u}_{3} = \mathbf{A} \vec{u}_{2} = \begin{pmatrix} 2.8863 \\ 5.3334 \\ 2.3216 \end{pmatrix}$$

So,

$$\vec{u}_{3} = \frac{\vec{u}_{3}}{|| \vec{u}_{3} ||} = \begin{pmatrix} 0.4445\\ 0.8213 \\ 0.3575 \end{pmatrix}$$

• For $k = 4$,

$$\vec{u}_{4} = \mathbf{A} \vec{u}_{3} = \begin{pmatrix} 2.9085 \\ 5.3532 \\ 2.2708 \end{pmatrix}$$

So,

$$\vec{u}_{4} = \frac{\vec{u}_{4}}{|| \vec{u}_{4} ||} = \begin{pmatrix} 0.4473\\ 0.8233 \\ 0.3493 \end{pmatrix}$$

• For two digits accuracy of the eigen vector, we can stop the iteration when $k = 4$. But for higher accuracy, we need to run the iteration until it converges to desired digits of accuracy.

• Therefore, dominant eigenvalue is $\lambda = \frac{(\mathbf{A} \vec{u}_{4}). (\vec{u}_{4})^{\text{T}}}{\vec{u}_{4}. (\vec{u}_{4})^{\text{T}} } = 6.5036$.

### Implementation

#### Using Python

#!/usr/bin/env python3

import numpy as np

def power_iteration(A, num_simulations: 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_simulations):
# 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

b_kn = power_iteration(np.array([[1, 3, 0], [2, 5, 1], [-1, 2, 3]]), 4)
print(b_kn)


## Theorems

• If $\mathbf{A}$ is diagonalizable and has dominant eigenvalue then, normalized power iteration sequence ($\mathbf{A}\vec{u}, \mathbf{A}^{2}\vec{u}, \mathbf{A}^{3}\vec{u}, \ldots$) converges to the dominant eigenvector. You can notice matrix $\mathbf{A}$ get power on it. This is why this method is called power iteration method or simply power method.
• If $\mathbf{A}$ is a square matrix with dominant eigenvalue $\lambda_{\text{max}}$ and $\vec{x}$ is any vector then, $|| \mathbf{A} x|| \leq |\lambda_{\text{max}} |~|| x ||$ where $|| \cdot ||$ is an Euclidean norm. This is very useful in estimates. This is why we are more interested in dominant eigenvalue.