Linear algebra/power iteration

linear_algebra/src/power_iteration.py

@@ -0,0 +1,96 @@
+import numpy as np
+
+
+def power_iteration(
+    input_matrix: np.array, vector: np.array, error_tol=1e-12, max_iterations=100
+) -> [float, np.array]:
+
+    """
+    Power Iteration.
+    Find the largest eignevalue and corresponding eigenvector
+    of matrix input_matrix given a random vector in the same space.
+    Will work so long as vector has component of largest eigenvector.
+    input_matrix must be symmetric.
+
+    Input
+    input_matrix: input matrix whose largest eigenvalue we will find.
+    Numpy array. np.shape(input_matrix) == (N,N).
+    vector: random initial vector in same space as matrix.
+    Numpy array. np.shape(vector) == (N,) or (N,1)
+
+    Output
+    largest_eigenvalue: largest eigenvalue of the matrix input_matrix.
+    Float. Scalar.
+    largest_eigenvector: eigenvector corresponding to largest_eigenvalue.
+    Numpy array. np.shape(largest_eigenvector) == (N,) or (N,1).
+
+    >>> import numpy as np
+    >>> A = np.array([
+    ... [41,  4, 20],
+    ... [ 4, 26, 30],
+    ... [20, 30, 50]
+    ... ])
+    >>> v = np.array([41,4,20])
+    >>> power_iteration(A,v)
+    (79.66086378788381, array([0.44472726, 0.46209842, 0.76725662]))
+    """
+
+    # Ensure matrix is square.
+    assert np.shape(input_matrix)[0] == np.shape(input_matrix)[1]
+    # Ensure proper dimensionality.
+    assert np.shape(input_matrix)[0] == np.shape(vector)[0]
+
+    # Set convergence to False. Will define convergence when we exceed max_iterations
+    # or when we have small changes from one iteration to next.
+
+    convergence = False
+    lamda_previous = 0
+    iterations = 0
+    error = 1e12
+
+    while not convergence:
+
+        # Multiple matrix by the vector.
+        w = np.dot(input_matrix, vector)
+        # Normalize the resulting output vector.
+        vector = w / np.linalg.norm(w)
+        # Find rayleigh quotient
+        # (faster than usual b/c we know vector is normalized already)
+        lamda = np.dot(vector.T, np.dot(input_matrix, vector))
+
+        # Check convergence.
+        error = np.abs(lamda - lamda_previous) / lamda
+        iterations += 1
+
+        if error <= error_tol or iterations >= max_iterations:
+            convergence = True
+
+        lamda_previous = lamda
+
+    return lamda, vector
+
+
+def tests() -> None:
+
+    # Our implementation.
+    A = np.array([[41, 4, 20], [4, 26, 30], [20, 30, 50]])
+    v = np.array([41, 4, 20])
+    eigen_value, eigen_vector = power_iteration(A, v)
+
+    # Numpy implementation.
+    eigs, eigvs = np.linalg.eigh(A)
+    eig_max = eigs[-1]
+    eigv_max = eigvs[:, -1]
+
+    # Check our implementation and numpy gives close answers.
+    assert np.abs(eigen_value - eig_max) <= 1e-6
+    # Take absolute values element wise of each eigenvector
+    # as they are only unique to a minus sign.
+    assert np.linalg.norm(np.abs(eigen_vector) - np.abs(eigv_max)) <= 1e-6
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
+    tests()