From f4886054b7a0bafb91fa139e588f69665338d1a4 Mon Sep 17 00:00:00 2001 From: Calvin McCarter Date: Sat, 15 Jan 2022 23:48:23 -0500 Subject: [PATCH 1/4] works python3 -m unittest discover --start-directory src --pattern "power*.py" --t . -v --- linear_algebra/src/power_iteration.py | 85 ++++++++++++++++++--------- 1 file changed, 57 insertions(+), 28 deletions(-) diff --git a/linear_algebra/src/power_iteration.py b/linear_algebra/src/power_iteration.py index 2cf22838e4a1..973cfcda3139 100644 --- a/linear_algebra/src/power_iteration.py +++ b/linear_algebra/src/power_iteration.py @@ -1,3 +1,5 @@ +import unittest + import numpy as np @@ -9,10 +11,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 +43,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 +65,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,37 +77,57 @@ def power_iteration( lamda_previous = lamda - return lamda, vector - + if is_complex: + lamda = np.real(lamda) -def test_power_iteration() -> None: - """ - >>> test_power_iteration() # self running tests - """ - # 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. + return lamda, vector - # Get eigen values and eigen vectors using built in numpy - # eigh (eigh used for symmetric or hermetian matrices). - eigen_values, eigen_vectors = np.linalg.eigh(input_matrix) - # Last eigen value is the maximum one. - eigen_value_max = eigen_values[-1] - # Last column in this matrix is eigen vector corresponding to largest eigen value. - eigen_vector_max = eigen_vectors[:, -1] - # Check our implementation and numpy gives close answers. - assert np.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. - assert np.linalg.norm(np.abs(eigen_vector) - np.abs(eigen_vector_max)) <= 1e-6 +class TestPowerIteration(unittest.TestCase): + def setUp(self) -> None: + self.real_input_matrix = np.array([[41, 4, 20], [4, 26, 30], [20, 30, 50]]) + self.real_vector = np.array([41, 4, 20]) + self.complex_input_matrix = self.real_input_matrix.astype(np.complex128) + imag_matrix = np.triu(1j * self.complex_input_matrix, 1) + self.complex_input_matrix += imag_matrix + self.complex_input_matrix += -1 * imag_matrix.T + self.complex_vector = np.array([41, 4, 20]).astype(np.complex128) + + + def test_power_iteration(self) -> None: + """ + >>> test_power_iteration() # self running tests + """ + for problem_type in ["real", "complex"]: + if problem_type == "real": + input_matrix = self.real_input_matrix + vector = self.real_vector + elif problem_type == "complex": + input_matrix = self.complex_input_matrix + vector = self.complex_vector + + # Our implementation. + eigen_value, eigen_vector = power_iteration(input_matrix, vector) + + # Numpy implementation. + + # Get eigen values and eigen vectors using built in numpy + # eigh (eigh used for symmetric or hermetian matrices). + eigen_values, eigen_vectors = np.linalg.eigh(input_matrix) + # Last eigen value is the maximum one. + eigen_value_max = eigen_values[-1] + # Last column in this matrix is eigen vector corresponding to largest eigen value. + eigen_vector_max = eigen_vectors[:, -1] + + # Check our implementation and numpy gives close answers. + assert np.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. + 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() + unittest.main() From 2532aa3d7876c21af10e85a9b517bf242a03d095 Mon Sep 17 00:00:00 2001 From: Calvin McCarter Date: Sat, 15 Jan 2022 23:48:56 -0500 Subject: [PATCH 2/4] cleanup --- linear_algebra/src/power_iteration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linear_algebra/src/power_iteration.py b/linear_algebra/src/power_iteration.py index 973cfcda3139..c3d403a8b9c0 100644 --- a/linear_algebra/src/power_iteration.py +++ b/linear_algebra/src/power_iteration.py @@ -130,4 +130,4 @@ def test_power_iteration(self) -> None: import doctest doctest.testmod() - unittest.main() + test_power_iteration() From 9fff131b7a5a598ecad54646f635b74101217c2d Mon Sep 17 00:00:00 2001 From: Calvin McCarter Date: Sun, 30 Jan 2022 15:47:58 -0500 Subject: [PATCH 3/4] revert switch to unittest --- linear_algebra/src/power_iteration.py | 81 +++++++++++++-------------- 1 file changed, 38 insertions(+), 43 deletions(-) diff --git a/linear_algebra/src/power_iteration.py b/linear_algebra/src/power_iteration.py index c3d403a8b9c0..12c1fb3d1994 100644 --- a/linear_algebra/src/power_iteration.py +++ b/linear_algebra/src/power_iteration.py @@ -1,5 +1,3 @@ -import unittest - import numpy as np @@ -83,47 +81,44 @@ def power_iteration( return lamda, vector -class TestPowerIteration(unittest.TestCase): - def setUp(self) -> None: - self.real_input_matrix = np.array([[41, 4, 20], [4, 26, 30], [20, 30, 50]]) - self.real_vector = np.array([41, 4, 20]) - self.complex_input_matrix = self.real_input_matrix.astype(np.complex128) - imag_matrix = np.triu(1j * self.complex_input_matrix, 1) - self.complex_input_matrix += imag_matrix - self.complex_input_matrix += -1 * imag_matrix.T - self.complex_vector = np.array([41, 4, 20]).astype(np.complex128) - - - def test_power_iteration(self) -> None: - """ - >>> test_power_iteration() # self running tests - """ - for problem_type in ["real", "complex"]: - if problem_type == "real": - input_matrix = self.real_input_matrix - vector = self.real_vector - elif problem_type == "complex": - input_matrix = self.complex_input_matrix - vector = self.complex_vector - - # Our implementation. - eigen_value, eigen_vector = power_iteration(input_matrix, vector) - - # Numpy implementation. - - # Get eigen values and eigen vectors using built in numpy - # eigh (eigh used for symmetric or hermetian matrices). - eigen_values, eigen_vectors = np.linalg.eigh(input_matrix) - # Last eigen value is the maximum one. - eigen_value_max = eigen_values[-1] - # Last column in this matrix is eigen vector corresponding to largest eigen value. - eigen_vector_max = eigen_vectors[:, -1] - - # Check our implementation and numpy gives close answers. - assert np.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. - assert np.linalg.norm(np.abs(eigen_vector) - np.abs(eigen_vector_max)) <= 1e-6 +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. + eigen_value, eigen_vector = power_iteration(input_matrix, vector) + + # Numpy implementation. + + # Get eigen values and eigen vectors using built in numpy + # eigh (eigh used for symmetric or hermetian matrices). + eigen_values, eigen_vectors = np.linalg.eigh(input_matrix) + # Last eigen value is the maximum one. + eigen_value_max = eigen_values[-1] + # Last column in this matrix is eigen vector corresponding to largest eigen value. + eigen_vector_max = eigen_vectors[:, -1] + + # Check our implementation and numpy gives close answers. + assert np.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. + assert np.linalg.norm(np.abs(eigen_vector) - np.abs(eigen_vector_max)) <= 1e-6 if __name__ == "__main__": From 75380544b5b0834a02ce102a412498c77e0372a9 Mon Sep 17 00:00:00 2001 From: Calvin McCarter Date: Mon, 31 Jan 2022 21:28:19 -0500 Subject: [PATCH 4/4] fix flake8 --- linear_algebra/src/power_iteration.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/linear_algebra/src/power_iteration.py b/linear_algebra/src/power_iteration.py index 12c1fb3d1994..4c6525b6e4af 100644 --- a/linear_algebra/src/power_iteration.py +++ b/linear_algebra/src/power_iteration.py @@ -106,12 +106,12 @@ def test_power_iteration() -> None: # 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 eigen value is the maximum one. + # Last eigenvalue is the maximum one. eigen_value_max = eigen_values[-1] - # Last column in this matrix is eigen vector corresponding to largest eigen value. + # 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.