Power Iteration for Dominant Eigenvector

To solve this problem, you mainly need to understand what eigenvalues/eigenvectors are and why repeated matrix–vector multiplication reveals the dominant one, then translate the math iteration into safe, numerically stable code.

---

1. Background Knowledge

An eigenvector of a square matrix $A \in \mathbb{R}^{n \times n}$ is a nonzero vector $\mathbf{v}$ such that

$$A\mathbf{v} = \lambda \mathbf{v},$$

where $\lambda$ is the corresponding eigenvalue. The dominant eigenvalue $\lambda_1$ is the eigenvalue with the largest magnitude $|$\lambda$_1|$; its eigenvector is the dominant eigenvector.

Power iteration relies on the fact that any vector can be written as a linear combination of eigenvectors (for diagonalizable $A$). If we write

$$\mathbf{v}_0 = c_1 \mathbf{u}_1 + c_2 \mathbf{u}_2 + \dots + c_n \mathbf{u}_n,$$

where $\mathbf{u}_i$ are eigenvectors and $\lambda_1$ is dominant, then after $k$ multiplications:

$$A^k \mathbf{v}_0 = c_1 \lambda_1^k \mathbf{u}_1 + c_2 \lambda_2^k \mathbf{u}_2 + \dots$$

The term with $\lambda_1^k$ dominates as $k$ grows (assuming $c_1 \neq 0$), so the direction of $A^k$\mathbf{v}$_0$ approaches the dominant eigenvector. Normalizing at each step prevents the vector from blowing up and keeps the iteration numerically stable.

The convergence rate depends on the eigenvalue gap: roughly proportional to $|$\lambda$_2 /$\lambda$_1|^k$, where $\lambda_2$ is the second-largest eigenvalue in magnitude. If this ratio is close to 1, convergence is slow; if it is small, convergence is fast.

---

2. Algorithm / General Approach

General pattern for power iteration:

1. Initialize with a random nonzero vector $\mathbf{v}_0$.
2. Repeat:

• Multiply by $A$ to get $\mathbf{w}_k = A$\mathbf{v}$_k$.
• Normalize: $\mathbf{v}_{k+1} =$\mathbf{w}$_k / \|$\mathbf{w}$_k\|_2$.

3. Check convergence by comparing $\mathbf{v}_{k+1}$ and $\mathbf{v}_k$.
4. After convergence, estimate the dominant eigenvalue $\lambda_1$ using the Rayleigh quotient:

$$\lambda_1 \approx \mathbf{v}^T A \mathbf{v}$$

when $\mathbf{v}$ is normalized.

You are essentially doing iterative matrix–vector multiplication plus normalization until the direction stabilizes.

---

3. Step-by-Step Strategy (Implementation-Oriented)

Assume you’re given a 2D array or tensor A (shape (n, n)), and you must return (eigenvalue, eigenvector).

1. Choose hyperparameters

• max_iters: maximum number of iterations (e.g., 1000).
• tol: tolerance for convergence (e.g., 1e-6).

2. Initialize vector

• Create a random vector:

v = np.random.randn(n)

• Normalize:

v = v / np.linalg.norm(v)

3. Iterative loop

for _ in range(max_iters):
w = A @ v # matrix-vector product
w_norm = np.linalg.norm(w)
if w_norm == 0:
# A sent v to zero; choose a new random vector or handle as special case
break
v_next = w / w_norm # normalize

# convergence check using vector difference
if np.linalg.norm(v_next - v) < tol:
v = v_next
break

v = v_next

Notes on convergence check:

• You can use np.linalg.norm(v_next - v) or check 1 - |v_next·v| (since eigenvectors are up to sign).
• Optionally store an error and stop when it’s below tol.

4. Compute dominant eigenvalue

• Using Rayleigh quotient with normalized v:

Av = A @ v
lambda_est = v.T @ Av

• Return lambda_est and v.

5. Package as a function

• Wrap the above into something like:

def power_iteration(A, max_iters=1000, tol=1e-6):
#... return lambda_est, v

---

4. Common Pitfalls

• No normalization each step:

• If you skip normalization, values can overflow or underflow, and the algorithm becomes unstable.

• Bad convergence check:

• Comparing eigenvalues too early is unstable; use vector change (‖v_{k+1} - v_k‖) or cosine similarity.

• Remember that eigenvectors are defined up to sign: v and -v are equivalent. A sign flip can make v_{k+1} - v_k large even if you’re converged. Using abs(v_next · v) or min of both norms can help:

if min(np.linalg.norm(v_next - v), np.linalg.norm(v_next + v)) < tol:
# converged

• Starting with an unlucky vector:

• If the initial vector is orthogonal to the dominant eigenvector, the method will fail. In practice, a random initialization makes this event probability zero in continuous spaces.

• Multiple eigenvalues with same magnitude:

• If the largest magnitude eigenvalue is repeated (e.g., ±same magnitude), the method may not converge to a unique direction, or may oscillate.

• Non-symmetric / pathological matrices:

• For general matrices, convergence can be slower or more complex, but for most typical CV/PCA-like problems (symmetric positive semidefinite matrices), it behaves well.

• Numerical precision:

• Using a very small tol (e.g., 1e-15) with float32 may never trigger; pick tolerance suited to your data type.

---

5. Time & Space Complexity

Let $n$ be the dimension of the matrix and $T$ the number of iterations until convergence.

• Per iteration:
• Matrix–vector multiplication $A\mathbf{v}$: $O(n^2)$ for dense $A$.
• Normalization and norms: $O(n)$.
• Total time complexity:

$$O(T \cdot n^2)$$

For sparse matrices with $m$ nonzero entries, the matvec cost is $O(m)$, so time is $O(T \cdot m)$.

• Space complexity:
• Storing $A$: $O(n^2)$ for dense, $O(m)$ for sparse.
• Storing a few vectors (v, w, Av): $O(n)$.
• Additional overhead is negligible.
• So extra space beyond the input matrix is $O(n)$.

This complexity is why power iteration is attractive in large-scale settings: it uses only matrix–vector products and $O(n)$ extra memory.