power iteration find the largest eigenvalue and corresponding eigenvector of matrix inputmatrix given a random vector in the same space will work so long as vector has component of largest eigenvector inputmatrix must be either real or hermitian input inputmatrix input matrix whose largest eigenvalue we will find numpy array np shapeinputmatrix n n vector random initial vector in same space as matrix numpy array np shapevector n or n 1 output largesteigenvalue largest eigenvalue of the matrix inputmatrix float scalar largesteigenvector eigenvector corresponding to largesteigenvalue numpy array np shapelargesteigenvector n or n 1 import numpy as np inputmatrix np array 41 4 20 4 26 30 20 30 50 vector np array41 4 20 poweriterationinputmatrix vector 79 66086378788381 array0 44472726 0 46209842 0 76725662 ensure matrix is square ensure proper dimensionality ensure inputs are either both complex or both real ensure complex inputmatrix is hermitian set convergence to false will define convergence when we exceed maxiterations or when we have small changes from one iteration to next multiple matrix by the vector normalize the resulting output vector find rayleigh quotient faster than usual bc we know vector is normalized already check convergence testpoweriteration self running tests our implementation numpy implementation get eigenvalues and eigenvectors using builtin numpy eigh eigh used for symmetric or hermetian matrices last eigenvalue is the maximum one last column in this matrix is eigenvector corresponding to largest eigenvalue check our implementation and numpy gives close answers take absolute values element wise of each eigenvector as they are only unique to a minus sign power iteration find the largest eigenvalue 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 either real or hermitian 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 input_matrix np array 41 4 20 4 26 30 20 30 50 vector np array 41 4 20 power_iteration input_matrix vector 79 66086378788381 array 0 44472726 0 46209842 0 76725662 ensure matrix is square ensure proper dimensionality ensure inputs are either both complex or both real ensure complex input_matrix is hermitian set convergence to false will define convergence when we exceed max_iterations or when we have small changes from one iteration to next multiple matrix by the vector normalize the resulting output vector find rayleigh quotient faster than usual b c we know vector is normalized already check convergence test_power_iteration self running tests our implementation numpy implementation get eigenvalues and eigenvectors using built in numpy eigh eigh used for symmetric or hermetian matrices last eigenvalue is the maximum one last column in this matrix is eigenvector corresponding to largest eigenvalue check our implementation and numpy gives close answers take absolute values element wise of each eigenvector as they are only unique to a minus sign
import numpy as np


def power_iteration(
    input_matrix: np.ndarray,
    vector: np.ndarray,
    error_tol: float = 1e-12,
    max_iterations: int = 100,
) -> tuple[float, np.ndarray]:
    

    
    assert np.shape(input_matrix)[0] == np.shape(input_matrix)[1]
    
    assert np.shape(input_matrix)[0] == np.shape(vector)[0]
    
    assert np.iscomplexobj(input_matrix) == np.iscomplexobj(vector)
    is_complex = np.iscomplexobj(input_matrix)
    if is_complex:
        
        assert np.array_equal(input_matrix, input_matrix.conj().T)

    
    

    convergence = False
    lambda_previous = 0
    iterations = 0
    error = 1e12

    while not convergence:
        
        w = np.dot(input_matrix, vector)
        
        vector = w / np.linalg.norm(w)
        
        
        vector_h = vector.conj().T if is_complex else vector.T
        lambda_ = np.dot(vector_h, np.dot(input_matrix, vector))

        
        error = np.abs(lambda_ - lambda_previous) / lambda_
        iterations += 1

        if error <= error_tol or iterations >= max_iterations:
            convergence = True

        lambda_previous = lambda_

    if is_complex:
        lambda_ = np.real(lambda_)

    return lambda_, vector


def test_power_iteration() -> None:
    
    real_input_matrix = np.array([[41, 4, 20], [4, 26, 30], [20, 30, 50]])
    real_vector = np.array([41, 4, 20])
    complex_input_matrix = real_input_matrix.astype(np.complex128)
    imag_matrix = np.triu(1j * complex_input_matrix, 1)
    complex_input_matrix += imag_matrix
    complex_input_matrix += -1 * imag_matrix.T
    complex_vector = np.array([41, 4, 20]).astype(np.complex128)

    for problem_type in ["real", "complex"]:
        if problem_type == "real":
            input_matrix = real_input_matrix
            vector = real_vector
        elif problem_type == "complex":
            input_matrix = complex_input_matrix
            vector = complex_vector

        
        eigen_value, eigen_vector = power_iteration(input_matrix, vector)

        

        
        
        eigen_values, eigen_vectors = np.linalg.eigh(input_matrix)
        
        eigen_value_max = eigen_values[-1]
        
        eigen_vector_max = eigen_vectors[:, -1]

        
        assert np.abs(eigen_value - eigen_value_max) <= 1e-6
        
        
        assert np.linalg.norm(np.abs(eigen_vector) - np.abs(eigen_vector_max)) <= 1e-6


if __name__ == "__main__":
    import doctest

    doctest.testmod()
    test_power_iteration()