power_iteration.py

importnumpyasnp


defpower_iteration(
input_matrix: np.ndarray,
vector: np.ndarray,
error_tol: float=1e-12,
max_iterations: int=100,
) ->tuple[float, np.ndarray]:
"""
 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.
assertnp.shape(input_matrix)[0] ==np.shape(input_matrix)[1]
# Ensure proper dimensionality.
assertnp.shape(input_matrix)[0] ==np.shape(vector)[0]
# Ensure inputs are either both complex or both real
assertnp.iscomplexobj(input_matrix) ==np.iscomplexobj(vector)
is_complex=np.iscomplexobj(input_matrix)
ifis_complex:
# Ensure complex input_matrix is Hermitian
assertnp.array_equal(input_matrix, input_matrix.conj().T)

# 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
lambda_previous=0
iterations=0
error=1e12

whilenotconvergence:
# 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)
vector_h=vector.conj().Tifis_complexelsevector.T
lambda_=np.dot(vector_h, np.dot(input_matrix, vector))

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

iferror<=error_toloriterations>=max_iterations:
convergence=True

lambda_previous=lambda_

ifis_complex:
lambda_=np.real(lambda_)

returnfloat(lambda_), vector


deftest_power_iteration() ->None:
"""
 >>> test_power_iteration() # self running tests
 """
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)

forproblem_typein ["real", "complex"]:
ifproblem_type=="real":
input_matrix=real_input_matrix
vector=real_vector
elifproblem_type=="complex":
input_matrix=complex_input_matrix
vector=complex_vector

# Our implementation.
eigen_value, eigen_vector=power_iteration(input_matrix, vector)

# Numpy implementation.

# Get eigenvalues and eigenvectors using built-in numpy
# eigh (eigh used for symmetric or hermetian matrices).
eigen_values, eigen_vectors=np.linalg.eigh(input_matrix)
# Last eigenvalue is the maximum one.
eigen_value_max=eigen_values[-1]
# Last column in this matrix is eigenvector corresponding to largest eigenvalue.
eigen_vector_max=eigen_vectors[:, -1]

# Check our implementation and numpy gives close answers.
assertnp.abs(eigen_value-eigen_value_max) <=1e-6
# Take absolute values element wise of each eigenvector.
# as they are only unique to a minus sign.
assertnp.linalg.norm(np.abs(eigen_vector) -np.abs(eigen_vector_max)) <=1e-6


if__name__=="__main__":
importdoctest

doctest.testmod()
test_power_iteration()