From 25517683873ec6ad1cac022074eb3d20f32fec32 Mon Sep 17 00:00:00 2001 From: Zeyad Zaky Date: Fri, 25 Sep 2020 23:55:56 -0700 Subject: [PATCH 1/7] Initial commit of conjugate gradient method --- linear_algebra/src/conjugate_gradient.py | 172 +++++++++++++++++++++++ 1 file changed, 172 insertions(+) create mode 100644 linear_algebra/src/conjugate_gradient.py diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py new file mode 100644 index 000000000000..ad0075cdd620 --- /dev/null +++ b/linear_algebra/src/conjugate_gradient.py @@ -0,0 +1,172 @@ +import numpy as np + + +def _is_matrix_spd(A: np.array) -> bool: + + """ + Returns True if input matrix A is symmetric positive definite. + Returns False otherwise. + + For a matrix to be SPD, all eigenvalues must be positive. + + >>> import numpy as np + >>> A = np.array([ + ... [4.12401784, -5.01453636, -0.63865857], + ... [-5.01453636, 12.33347422, -3.40493586], + ... [-0.63865857, -3.40493586, 5.78591885]]) + >>> _is_matrix_spd(A) + True + >>> A = np.array([ + ... [0.34634879, 1.96165514, 2.18277744], + ... [0.74074469, -1.19648894, -1.34223498], + ... [-0.7687067 , 0.06018373, -1.16315631]]) + >>> _is_matrix_spd(A) + False + """ + # Ensure matrix is square. + assert np.shape(A)[0] == np.shape(A)[1] + + # Get eigenvalues and eignevectors for a symmetric matrix. + eigen_values, _ = np.linalg.eigh(A) + + # Check sign of all eigenvalues. + if np.all(eigen_values > 0): + return True + else: + return False + + +def _create_spd_matrix(N: np.int64) -> np.array: + """ + Returns a symmetric positive definite matrix given a dimension. + + Input: + N is an integer. + + Output: + A is an NxN symmetric positive definite (SPD) matrix. + + >>> import numpy as np + >>> N = 3 + >>> A = _create_spd_matrix(N) + >>> _is_matrix_spd(A) + True + """ + + A = np.random.randn(N, N) + A = np.dot(A, A.T) + + assert _is_matrix_spd(A) is True + + return A + + +def conjugate_gradient( + A: np.array, b: np.array, max_iterations=1000, tol=1e-8 +) -> np.array: + """ + Returns solution to the linear system Ax = b. + + Input: + A is an NxN Symmetric Positive Definite (SPD) matrix. + b is an Nx1 vector. + + Output: + x is an Nx1 vector. + + >>> import numpy as np + >>> A = np.array([ + ... [8.73256573, -5.02034289, -2.68709226], + ... [-5.02034289, 3.78188322, 0.91980451], + ... [-2.68709226, 0.91980451, 1.94746467]]) + >>> b = np.array([ + ... [-5.80872761], + ... [ 3.23807431], + ... [ 1.95381422]]) + >>> conjugate_gradient(A,b) + array([[-0.63114139], + [-0.01561498], + [ 0.13979294]]) + """ + # Ensure proper dimensionality. + assert np.shape(A)[0] == np.shape(A)[1] + assert np.shape(b)[0] == np.shape(A)[0] + assert _is_matrix_spd(A) + + N = np.shape(b)[0] + + # Initialize solution guess, residual, search direction. + x0 = np.zeros((N, 1)) + r0 = np.copy(b) + p0 = np.copy(r0) + + # Set initial errors in solution guess and residual. + error_residual = 1e9 + error_x_solution = 1e9 + error = 1e9 + + # Set iteration counter to threshold number of iterations. + iterations = 0 + + while error > tol: + + # Save this value so we only calculate the matrix-vector product once. + w = np.dot(A, p0) + + # The main algorithm. + + # Update search direction magnitude. + alpha = np.dot(r0.T, r0) / np.dot(p0.T, w) + # Update solution guess. + x = x0 + alpha * p0 + # Calculate new residual. + r = r0 - alpha * w + # Calculate new Krylov subspace scale. + beta = np.dot(r.T, r) / np.dot(r0.T, r0) + # Calculate new A conjuage search direction. + p = r + beta * p0 + + # Calculate errors. + error_residual = np.linalg.norm(r - r0) + error_x_solution = np.linalg.norm(x - x0) + error = np.maximum(error_residual, error_x_solution) + + # Update variables. + x0 = np.copy(x) + r0 = np.copy(r) + p0 = np.copy(p) + + # Update number of iterations. + iterations += 1 + + return x + + +def test_conjugate_gradient() -> None: + + """ + >>> test_conjugate_gradient() # self running tests + """ + + # Create linear system with SPD matrix and known solution x_true. + N = 3 + A = _create_spd_matrix(N) + x_true = np.random.randn(N, 1) + b = np.dot(A, x_true) + + # Numpy solution. + x_numpy = np.linalg.solve(A, b) + + # Our implementation. + x_conjugate_gradient = conjugate_gradient(A, b) + + # Ensure both solutions are close to x_true (and therefore one another). + assert np.linalg.norm(x_numpy - x_true) <= 1e-6 + assert np.linalg.norm(x_conjugate_gradient - x_true) <= 1e-6 + + +if __name__ == "__main__": + import doctest + + doctest.testmod() + test_conjugate_gradient() From 6d2c51c2b5a2bec2026d1150b9e1c58b5fb4b9c3 Mon Sep 17 00:00:00 2001 From: zakademic <67771932+zakademic@users.noreply.github.com> Date: Mon, 28 Sep 2020 23:28:11 -0700 Subject: [PATCH 2/7] Update linear_algebra/src/conjugate_gradient.py Co-authored-by: Christian Clauss --- linear_algebra/src/conjugate_gradient.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index ad0075cdd620..0d4c10533085 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -30,10 +30,7 @@ def _is_matrix_spd(A: np.array) -> bool: eigen_values, _ = np.linalg.eigh(A) # Check sign of all eigenvalues. - if np.all(eigen_values > 0): - return True - else: - return False + return np.all(eigen_values > 0) def _create_spd_matrix(N: np.int64) -> np.array: From 0613b60fe8a84cfd277e7ae0ecfc7a58375307d6 Mon Sep 17 00:00:00 2001 From: zakademic <67771932+zakademic@users.noreply.github.com> Date: Mon, 28 Sep 2020 23:28:21 -0700 Subject: [PATCH 3/7] Update linear_algebra/src/conjugate_gradient.py Co-authored-by: Christian Clauss --- linear_algebra/src/conjugate_gradient.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index 0d4c10533085..c555a1503952 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -53,7 +53,7 @@ def _create_spd_matrix(N: np.int64) -> np.array: A = np.random.randn(N, N) A = np.dot(A, A.T) - assert _is_matrix_spd(A) is True + assert _is_matrix_spd(A) return A From 322cfa4d8f9cb20cb054390125e83e3fd541be95 Mon Sep 17 00:00:00 2001 From: Zeyad Zaky Date: Mon, 28 Sep 2020 23:40:42 -0700 Subject: [PATCH 4/7] Added documentation links, chnaged variable names to lower case and more descriptive naming, added check for symmetry in _is_matrix_spd --- linear_algebra/src/conjugate_gradient.py | 69 +++++++++++++----------- 1 file changed, 38 insertions(+), 31 deletions(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index c555a1503952..39512ef787a1 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -1,33 +1,40 @@ import numpy as np +# https://en.wikipedia.org/wiki/Conjugate_gradient_method +# https://en.wikipedia.org/wiki/Definite_symmetric_matrix -def _is_matrix_spd(A: np.array) -> bool: + +def _is_matrix_spd(matrix: np.array) -> bool: """ - Returns True if input matrix A is symmetric positive definite. + Returns True if input matrix is symmetric positive definite. Returns False otherwise. For a matrix to be SPD, all eigenvalues must be positive. >>> import numpy as np - >>> A = np.array([ + >>> matrix = np.array([ ... [4.12401784, -5.01453636, -0.63865857], ... [-5.01453636, 12.33347422, -3.40493586], ... [-0.63865857, -3.40493586, 5.78591885]]) - >>> _is_matrix_spd(A) + >>> _is_matrix_spd(matrix) True - >>> A = np.array([ + >>> matrix = np.array([ ... [0.34634879, 1.96165514, 2.18277744], ... [0.74074469, -1.19648894, -1.34223498], ... [-0.7687067 , 0.06018373, -1.16315631]]) - >>> _is_matrix_spd(A) + >>> _is_matrix_spd(matrix) False """ # Ensure matrix is square. - assert np.shape(A)[0] == np.shape(A)[1] + assert np.shape(matrix)[0] == np.shape(matrix)[1] + + # If matrix not symmetric, exit right away. + if np.allclose(matrix, matrix.T) is False: + return False # Get eigenvalues and eignevectors for a symmetric matrix. - eigen_values, _ = np.linalg.eigh(A) + eigen_values, _ = np.linalg.eigh(matrix) # Check sign of all eigenvalues. return np.all(eigen_values > 0) @@ -38,41 +45,41 @@ def _create_spd_matrix(N: np.int64) -> np.array: Returns a symmetric positive definite matrix given a dimension. Input: - N is an integer. + N gives the square matrix dimension. Output: - A is an NxN symmetric positive definite (SPD) matrix. + spd_matrix is an NxN symmetric positive definite (SPD) matrix. >>> import numpy as np >>> N = 3 - >>> A = _create_spd_matrix(N) - >>> _is_matrix_spd(A) + >>> spd_matrix = _create_spd_matrix(N) + >>> _is_matrix_spd(spd_matrix) True """ - A = np.random.randn(N, N) - A = np.dot(A, A.T) + random_matrix = np.random.randn(N, N) + spd_matrix = np.dot(random_matrix, random_matrix.T) - assert _is_matrix_spd(A) + assert _is_matrix_spd(spd_matrix) - return A + return spd_matrix def conjugate_gradient( - A: np.array, b: np.array, max_iterations=1000, tol=1e-8 + spd_matrix: np.array, b: np.array, max_iterations=1000, tol=1e-8 ) -> np.array: """ - Returns solution to the linear system Ax = b. + Returns solution to the linear system np.dot(spd_matrix, x) = b. Input: - A is an NxN Symmetric Positive Definite (SPD) matrix. - b is an Nx1 vector. + spd_matrix is an NxN Symmetric Positive Definite (SPD) matrix. + b is an Nx1 vector that is the load vector. Output: - x is an Nx1 vector. + x is an Nx1 vector that is the solution vector. >>> import numpy as np - >>> A = np.array([ + >>> spd_matrix = np.array([ ... [8.73256573, -5.02034289, -2.68709226], ... [-5.02034289, 3.78188322, 0.91980451], ... [-2.68709226, 0.91980451, 1.94746467]]) @@ -80,15 +87,15 @@ def conjugate_gradient( ... [-5.80872761], ... [ 3.23807431], ... [ 1.95381422]]) - >>> conjugate_gradient(A,b) + >>> conjugate_gradient(spd_matrix,b) array([[-0.63114139], [-0.01561498], [ 0.13979294]]) """ # Ensure proper dimensionality. - assert np.shape(A)[0] == np.shape(A)[1] - assert np.shape(b)[0] == np.shape(A)[0] - assert _is_matrix_spd(A) + assert np.shape(spd_matrix)[0] == np.shape(spd_matrix)[1] + assert np.shape(b)[0] == np.shape(spd_matrix)[0] + assert _is_matrix_spd(spd_matrix) N = np.shape(b)[0] @@ -108,7 +115,7 @@ def conjugate_gradient( while error > tol: # Save this value so we only calculate the matrix-vector product once. - w = np.dot(A, p0) + w = np.dot(spd_matrix, p0) # The main algorithm. @@ -147,15 +154,15 @@ def test_conjugate_gradient() -> None: # Create linear system with SPD matrix and known solution x_true. N = 3 - A = _create_spd_matrix(N) + spd_matrix = _create_spd_matrix(N) x_true = np.random.randn(N, 1) - b = np.dot(A, x_true) + b = np.dot(spd_matrix, x_true) # Numpy solution. - x_numpy = np.linalg.solve(A, b) + x_numpy = np.linalg.solve(spd_matrix, b) # Our implementation. - x_conjugate_gradient = conjugate_gradient(A, b) + x_conjugate_gradient = conjugate_gradient(spd_matrix, b) # Ensure both solutions are close to x_true (and therefore one another). assert np.linalg.norm(x_numpy - x_true) <= 1e-6 From a602582789778a69b689d88ee964463e86db88b1 Mon Sep 17 00:00:00 2001 From: zakademic <67771932+zakademic@users.noreply.github.com> Date: Mon, 23 Nov 2020 11:43:23 -0800 Subject: [PATCH 5/7] Update linear_algebra/src/conjugate_gradient.py Co-authored-by: Christian Clauss --- linear_algebra/src/conjugate_gradient.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index 39512ef787a1..a44caebf1a8a 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -97,10 +97,8 @@ def conjugate_gradient( assert np.shape(b)[0] == np.shape(spd_matrix)[0] assert _is_matrix_spd(spd_matrix) - N = np.shape(b)[0] - # Initialize solution guess, residual, search direction. - x0 = np.zeros((N, 1)) + x0 = np.zeros((np.shape(b)[0], 1)) r0 = np.copy(b) p0 = np.copy(r0) From bb8f9003a2ce4d13ffb2583dfa4372b8739e0b0c Mon Sep 17 00:00:00 2001 From: Zeyad Zaky Date: Mon, 23 Nov 2020 11:50:53 -0800 Subject: [PATCH 6/7] Made changes to some variable naming to be more clear --- linear_algebra/src/conjugate_gradient.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index a44caebf1a8a..2443b897b60a 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -40,24 +40,24 @@ def _is_matrix_spd(matrix: np.array) -> bool: return np.all(eigen_values > 0) -def _create_spd_matrix(N: np.int64) -> np.array: +def _create_spd_matrix(dimension: np.int64) -> np.array: """ Returns a symmetric positive definite matrix given a dimension. Input: - N gives the square matrix dimension. + dimension gives the square matrix dimension. Output: - spd_matrix is an NxN symmetric positive definite (SPD) matrix. + spd_matrix is an diminesion x dimensions symmetric positive definite (SPD) matrix. >>> import numpy as np - >>> N = 3 - >>> spd_matrix = _create_spd_matrix(N) + >>> dimension = 3 + >>> spd_matrix = _create_spd_matrix(dimension) >>> _is_matrix_spd(spd_matrix) True """ - random_matrix = np.random.randn(N, N) + random_matrix = np.random.randn(dimension, dimension) spd_matrix = np.dot(random_matrix, random_matrix.T) assert _is_matrix_spd(spd_matrix) @@ -151,9 +151,9 @@ def test_conjugate_gradient() -> None: """ # Create linear system with SPD matrix and known solution x_true. - N = 3 - spd_matrix = _create_spd_matrix(N) - x_true = np.random.randn(N, 1) + dimension = 3 + spd_matrix = _create_spd_matrix(dimension) + x_true = np.random.randn(dimension, 1) b = np.dot(spd_matrix, x_true) # Numpy solution. From b3bc4d9f3c8da0a812a614a3553f69423cb02c36 Mon Sep 17 00:00:00 2001 From: Dhruv Manilawala Date: Thu, 10 Dec 2020 23:46:08 +0530 Subject: [PATCH 7/7] Update conjugate_gradient.py --- linear_algebra/src/conjugate_gradient.py | 29 ++++++++++++------------ 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/linear_algebra/src/conjugate_gradient.py b/linear_algebra/src/conjugate_gradient.py index 2443b897b60a..1a65b8ccf019 100644 --- a/linear_algebra/src/conjugate_gradient.py +++ b/linear_algebra/src/conjugate_gradient.py @@ -1,11 +1,12 @@ +""" +Resources: +- https://en.wikipedia.org/wiki/Conjugate_gradient_method +- https://en.wikipedia.org/wiki/Definite_symmetric_matrix +""" import numpy as np -# https://en.wikipedia.org/wiki/Conjugate_gradient_method -# https://en.wikipedia.org/wiki/Definite_symmetric_matrix - def _is_matrix_spd(matrix: np.array) -> bool: - """ Returns True if input matrix is symmetric positive definite. Returns False otherwise. @@ -56,24 +57,24 @@ def _create_spd_matrix(dimension: np.int64) -> np.array: >>> _is_matrix_spd(spd_matrix) True """ - random_matrix = np.random.randn(dimension, dimension) spd_matrix = np.dot(random_matrix, random_matrix.T) - assert _is_matrix_spd(spd_matrix) - return spd_matrix def conjugate_gradient( - spd_matrix: np.array, b: np.array, max_iterations=1000, tol=1e-8 + spd_matrix: np.array, + load_vector: np.array, + max_iterations: int = 1000, + tol: float = 1e-8, ) -> np.array: """ Returns solution to the linear system np.dot(spd_matrix, x) = b. Input: spd_matrix is an NxN Symmetric Positive Definite (SPD) matrix. - b is an Nx1 vector that is the load vector. + load_vector is an Nx1 vector. Output: x is an Nx1 vector that is the solution vector. @@ -87,19 +88,19 @@ def conjugate_gradient( ... [-5.80872761], ... [ 3.23807431], ... [ 1.95381422]]) - >>> conjugate_gradient(spd_matrix,b) + >>> conjugate_gradient(spd_matrix, b) array([[-0.63114139], [-0.01561498], [ 0.13979294]]) """ # Ensure proper dimensionality. assert np.shape(spd_matrix)[0] == np.shape(spd_matrix)[1] - assert np.shape(b)[0] == np.shape(spd_matrix)[0] + assert np.shape(load_vector)[0] == np.shape(spd_matrix)[0] assert _is_matrix_spd(spd_matrix) # Initialize solution guess, residual, search direction. - x0 = np.zeros((np.shape(b)[0], 1)) - r0 = np.copy(b) + x0 = np.zeros((np.shape(load_vector)[0], 1)) + r0 = np.copy(load_vector) p0 = np.copy(r0) # Set initial errors in solution guess and residual. @@ -145,11 +146,9 @@ def conjugate_gradient( def test_conjugate_gradient() -> None: - """ >>> test_conjugate_gradient() # self running tests """ - # Create linear system with SPD matrix and known solution x_true. dimension = 3 spd_matrix = _create_spd_matrix(dimension)