# ML Wiki

## Power Iteration

This is a linear algebra method for computing Eigenvalues and Eigenvectors

## Symmetric Matrices

Suppose $A$ is symmetric

• then Eigendecomposition of $A$ is $A = Q \Lambda Q^T$
• and $A^k = Q \Lambda^k Q^T$
• let $\mathbf q_i$ be the columns of $Q$

When $k \to \infty$

• suppose $\lambda_1 > \lambda_2$,
• then, as $k$ increases, $\lambda_1^k / \lambda_i^k \to 0$ for $i \geqslant 2$
• so $A^k \to \lambda_1^k \mathbf q_1 \mathbf q_1^T$
• so we can find $\lambda_1$ by powering $A$

Performance:

• computing $A \cdot A^{k-1}$ is costly
• instead, can compute $A \, A^{k-1} \mathbf x$ where $\mathbf x$ is some random unit vector
• $A^k \mathbf x \approx \lambda_1^k \, \mathbf q_1 \, (\mathbf q_1^T \mathbf x) = \alpha \cdot \mathbf q_1$
• so it's some scalar times a unit vector $\mathbf q_1$
• thus, to recover $\mathbf q_i$, use $A^k \mathbf x / \| A^k \mathbf x \| = \mathbf q_i$

## Finding Other Eigenvectors

• $A = Q \Lambda Q^T$, so $A = \sum_{i = 1}^n \lambda_i \mathbf q_i \mathbf q_i^T$
• use power iteration to find $\mathbf q_1$ and $\lambda_1$
• then let $A_2 \leftarrow A - \lambda_1 \mathbf q_1 \mathbf q_1^T$
• repeat power iteration on $A_2$ to find $\mathbf q_2$ and $\lambda_2$
• continue like this for $\lambda_3, \ ... \ , \lambda_n$

## Issues

• if $\lambda_1 - \lambda_2 \gg 0$, then it's fine, but if $\lambda_1 - \lambda_2$ is small, it will take a lot of time to converge

## Implementation

### NumPy

def eigenvalue(S, w):
return w.T.dot(S).dot(w)

def power_iteration(X, debug=True):
n, d = X.shape
converged = False

S = X.T.dot(X)
w = np.ones(d) / np.sqrt(d)
error = eigenvalue(S, w)

while not converged:
Sw = S.dot(w)
w_new = Sw / np.linalg.norm(Sw)

error_new = eigenvalue(S, w_new)
converged = np.abs(error - error_new) < 0.01

w = w_new
error = error_new

return w


### Mahout Samsara

var x: Vector = 1 to dim map (_ => 1.0 / Math.sqrt(dim))
var converged = false

while (!converged) {
val Ax = A %*% x
var x_new = Ax.collect(::, 0)
x_new = x_new / x_new.norm(2)

val diff = (x_new - x).norm(2)
converged = diff < 1e-6
x = x_new
}