Power’s iterative method - Executable Engineering: Interactive Numerical Methods in Python

We’ve seen that linear systems transform a vector into another vector by scaling and skewing it, and that the magnitude of the scale in a particular direction is given by the eigenvalues in that direction.

Largest (dominant) eigenvalue¶

The Power Method observes that a random vector is scaled along all eigenvectors, but it is scaled the most along the vector with the largest eigenvalue. Therefore, repeated application of $A$ followed by normalization of the vector will gradually rotate all vectors towards the eigenvector with the dominant eigenvalue!

The algorithm is:

Choose a random vector
Apply $A$
Normalize the transformed vector
Repeat until a subsequent vectors are the same (within a tolerance)

Note this involves copious matrix-matrix multiplications, which can be a) computationally intensive, and b) prone to accumulation of roundoff error.

Nonetheless, this algorithm is used extensively for situations where only the largest EV is needed. E.g.: The Google PageRank algorithm determines the dominant Eigenvector of a matrix associted withe the web’s link structure!

Example of the Power Method for finding the dominant eigenvalue - vector¶

# prompt: give me a function that implements Power's method for finding the largest eigenvalue. keep the vector from each iteration. Plot the vector indexed by a slider.

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

def power_method(A, initial_vector, iterations=100, tolerance=1e-6):
    """
    Implements the Power Method to find the dominant eigenvalue and eigenvector.

    Args:
        A: The square matrix.
        initial_vector: The initial guess for the eigenvector.
        iterations: The maximum number of iterations.
        tolerance: The convergence tolerance.

    Returns:
        A tuple containing:
            - The dominant eigenvalue.
            - The corresponding eigenvector.
            - A list of vectors from each iteration.
    """

    v = np.array(initial_vector, dtype=float)  # Ensure initial vector is a NumPy array and float type
    vectors = [v.copy()]  # Store the vectors from each iteration

    for _ in range(iterations):
        v_new = A @ v
        eigenvalue = np.linalg.norm(v_new) #np.max(np.abs(v_new))  # Approximate eigenvalue
        v = v_new / eigenvalue  # Normalize
        vectors.append(v.copy())

        if np.linalg.norm(v - vectors[-2]) < tolerance:
          break
    return eigenvalue, v, vectors

def plot_eigenvector(vectors, iteration):
  """Plots the eigenvector at a specific iteration."""
  plt.figure(figsize=(6, 6))
  plt.title(f"Eigenvector at Iteration {iteration}")
  plt.quiver(0, 0, vectors[iteration][0], vectors[iteration][1], angles='xy', scale_units='xy', scale=1, color='r')
  plt.xlim(-1.5, 1.5)
  plt.ylim(-1.5, 1.5)
  plt.gca().set_aspect('equal')
  plt.xlabel("X")
  plt.ylabel("Y")
  plt.grid(True)
  plt.show()

# Example usage:
A = np.array([[2, -1], [-1, 2]])
initial_vector = np.random.rand(2,1) # Example initial vector
print(initial_vector)
iterations = 10

eigenvalue, eigenvector, vectors = power_method(A, initial_vector, iterations)

print(f"Dominant eigenvalue: {eigenvalue}")
print(f"Eigenvector: {eigenvector}")

interact(plot_eigenvector, vectors=fixed(vectors), iteration=widgets.IntSlider(min=0, max=len(vectors)-1, step=1, value=0));

[[0.07024765]
 [0.31961613]]
Dominant eigenvalue: 2.9999999915880085
Eigenvector: [[-0.70708806]
 [ 0.7071255 ]]

Smallest eigenvalue¶

The smallest eigenvalue of $A$ is simply the largest eigenvalue of $A^{-1}$ , so all we have to do is calculate $A^{-1}$ and apply the algorithm above!

The smallest eigenvalue is of practical importance in a few places:

Stability analysis: The sign of the smallest EV determines the stability of systems (consider Von Neumann analysis!)
Condition number: The ratio of the largest and smallest EVs is a measure of the condition numeber of the matrix.
Optimization: The smallest EV is related to the convexivity of the function.

Intermediate eigenvalues¶

We can remove an eigenvector from the operator through deflation. The eigenvectors of symmetric matricies are necessarily orthogonal, so we can deflate the operator by:

$A^\prime = A - \lambda_1 x_1 x_1^T$

Therefore intermediate vectors can be found through iteratively removing the dominant one. This method is uncommon since roundoff errors are propagated between successive eigenvectors.

Example of deflation and finding the second largest eigenvalue/vector.¶

A_prime = A-eigenvalue*np.outer(eigenvector,eigenvector)
print('The deflated matrix is,\n', A_prime, '\n')

eigenvalue, eigenvector, vectors = power_method(A_prime, initial_vector, iterations)

print(f"Dominant eigenvalue: {eigenvalue}")
print(f"Eigenvector: {eigenvector}")

interact(plot_eigenvector, vectors=fixed(vectors), iteration=widgets.IntSlider(min=0, max=len(vectors)-1, step=1, value=0));

The deflated matrix is,
 [[0.50007943 0.49999999]
 [0.49999999 0.49992058]] 

Dominant eigenvalue: 1.000000004205996
Eigenvector: [[0.70716294]
 [0.70705061]]