mathdatech/Python
1
0
Fork 0
mirror of https://github.com/TheAlgorithms/Python.git synced 2024-11-23 21:11:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

Extend power iteration to handle complex Hermitian input matrices (#5925)

Browse Source 
* works

python3 -m unittest discover --start-directory src --pattern "power*.py" --t . -v

* cleanup

* revert switch to unittest

* fix flake8
This commit is contained in:
Calvin McCarter 2022-02-02 15:05:05 -05:00 committed by GitHub
parent d28ac6483a
commit 54f765bdd0
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 44 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
linear_algebra/src/power_iteration.py
View File

@ -9,10 +9,10 @@ def power_iteration(
) -> tuple[float, np.ndarray]:
"""
Power Iteration.
Find the largest eignevalue and corresponding eigenvector
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 symmetric.
input_matrix must be either real or Hermitian.
Input
input_matrix: input matrix whose largest eigenvalue we will find.
@ -41,6 +41,12 @@ def power_iteration(
assert np.shape(input_matrix)[0] == np.shape(input_matrix)[1]
# Ensure proper dimensionality.
assert np.shape(input_matrix)[0] == np.shape(vector)[0]
# Ensure inputs are either both complex or both real
assert np.iscomplexobj(input_matrix) == np.iscomplexobj(vector)
is_complex = np.iscomplexobj(input_matrix)
if is_complex:
# Ensure complex input_matrix is Hermitian
assert np.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.
@ -57,7 +63,8 @@ def power_iteration(
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))
vectorH = vector.conj().T if is_complex else vector.T
lamda = np.dot(vectorH, np.dot(input_matrix, vector))
# Check convergence.
error = np.abs(lamda - lamda_previous) / lamda
@ -68,6 +75,9 @@ def power_iteration(
lamda_previous = lamda
if is_complex:
lamda = np.real(lamda)
return lamda, vector
@ -75,14 +85,28 @@ def test_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)
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
# Our implementation.
input_matrix = np.array([[41, 4, 20], [4, 26, 30], [20, 30, 50]])
vector = np.array([41, 4, 20])
eigen_value, eigen_vector = power_iteration(input_matrix, vector)
# Numpy implementation.
# Get eigen values and eigen vectors using built in numpy
# 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.