From 8019213943b1420c4f75c6724fb909dc7732a8c1 Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Mon, 7 Oct 2024 11:17:58 +0200 Subject: [PATCH 1/6] Add cholesky_decomposition --- maths/cholesky_decomposition.py | 102 ++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 maths/cholesky_decomposition.py diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py new file mode 100644 index 000000000000..fe0820f9a5c6 --- /dev/null +++ b/maths/cholesky_decomposition.py @@ -0,0 +1,102 @@ +import numpy as np + + +def cholesky_decomposition(a: np.ndarray) -> np.ndarray: + """Return a Cholesky decomposition of the matrix A. + + The Cholesky decomposition decomposes the square, positive definite matrix A + into a lower triangular matrix L such that A = L L^T. + + https://en.wikipedia.org/wiki/Cholesky_decomposition + + Arguments: + A -- a numpy.ndarray of shape (n, n) + + >>> A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float) + >>> L = cholesky_decomposition(A) + >>> np.allclose(L, np.array([[2, 0, 0], [6, 1, 0], [-8, 5, 3]])) + True + + >>> # check that the decomposition is correct + >>> np.allclose(L @ L.T, A) + True + + >>> # check that L is lower triangular + >>> np.allclose(np.tril(L), L) + True + + The Cholesky decomposition can be used to solve the system of equations A x = y. + + >>> x_true = np.array([1, 2, 3], dtype=float) + >>> y = A @ x_true + >>> x = solve_cholesky(L, y) + >>> np.allclose(x, x_true) + True + + It can also be used to solve multiple equations A X = Y simultaneously. + + >>> X_true = np.random.rand(3, 3) + >>> Y = A @ X_true + >>> X = solve_cholesky(L, Y) + >>> np.allclose(X, X_true) + True + """ + assert a.shape[0] == a.shape[1] + n = a.shape[0] + lo = np.tril(a) + + for i in range(n): + for j in range(i): + lo[i, j] = (lo[i, j] - np.sum(lo[i, :j] * lo[j, :j])) / lo[j, j] + + s = lo[i, i] - np.sum(lo[i, :i] * lo[i, :i]) + + if s <= 0: + raise ValueError("Matrix A is not positive definite") + + lo[i, i] = np.sqrt(s) + + return lo + + +def solve_cholesky(lo: np.ndarray, y: np.ndarray) -> np.ndarray: + """Given a Cholesky decomposition L L^T = A of a matrix A, solve the + system of equations A X = Y where B is either a matrix or a vector. + + >>> L = np.array([[2, 0], [3, 4]], dtype=float) + >>> Y = np.array([[22, 54], [81, 193]], dtype=float) + >>> X = solve_cholesky(L, Y) + >>> np.allclose(X, np.array([[1, 3], [3, 7]], dtype=float)) + True + """ + + # Handle vector case by reshaping to matrix and then flattening again + if len(y.shape) == 1: + return solve_cholesky(lo, y.reshape(-1, 1)).ravel() + + n, m = y.shape + + # Backsubstitute L X = B + x = y.copy() + for i in range(n): + for j in range(i): + x[i, :] -= lo[i, j] * x[j, :] + + for k in range(m): + x[i, k] /= lo[i, i] + + # Backsubstitute L^T + for i in reversed(range(n)): + for j in range(i + 1, n): + x[i, :] -= lo[j, i] * x[j, :] + + for k in range(m): + x[i, k] /= lo[i, i] + + return x + + +if __name__ == "__main__": + import doctest + + doctest.testmod() From 307ce1cd7f8cc999d76140574c6798df3d874b4f Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Wed, 9 Oct 2024 08:32:28 +0200 Subject: [PATCH 2/6] Simplify equations, rename variables --- maths/cholesky_decomposition.py | 58 ++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 26 deletions(-) diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py index fe0820f9a5c6..363436b40ac3 100644 --- a/maths/cholesky_decomposition.py +++ b/maths/cholesky_decomposition.py @@ -1,7 +1,8 @@ import numpy as np -def cholesky_decomposition(a: np.ndarray) -> np.ndarray: +# ruff: noqa: N803,N806 +def cholesky_decomposition(A: np.ndarray) -> np.ndarray: """Return a Cholesky decomposition of the matrix A. The Cholesky decomposition decomposes the square, positive definite matrix A @@ -41,25 +42,28 @@ def cholesky_decomposition(a: np.ndarray) -> np.ndarray: >>> np.allclose(X, X_true) True """ - assert a.shape[0] == a.shape[1] - n = a.shape[0] - lo = np.tril(a) - for i in range(n): - for j in range(i): - lo[i, j] = (lo[i, j] - np.sum(lo[i, :j] * lo[j, :j])) / lo[j, j] + assert A.shape[0] == A.shape[1], f"A is not square, {A.shape=}" - s = lo[i, i] - np.sum(lo[i, :i] * lo[i, :i]) + n = A.shape[0] + L = np.tril(A) - if s <= 0: - raise ValueError("Matrix A is not positive definite") + for i in range(n): + for j in range(i + 1): + L[i, j] -= np.sum(L[i, :j] * L[j, :j]) - lo[i, i] = np.sqrt(s) + if i == j: + if L[i, i] <= 0: + raise ValueError("Matrix A is not positive definite") - return lo + L[i, i] = np.sqrt(L[i, i]) + else: + L[i, j] /= L[j, j] + return L -def solve_cholesky(lo: np.ndarray, y: np.ndarray) -> np.ndarray: + +def solve_cholesky(L: np.ndarray, Y: np.ndarray) -> np.ndarray: """Given a Cholesky decomposition L L^T = A of a matrix A, solve the system of equations A X = Y where B is either a matrix or a vector. @@ -70,30 +74,32 @@ def solve_cholesky(lo: np.ndarray, y: np.ndarray) -> np.ndarray: True """ + assert L.shape[0] == L.shape[1], f"L is not square, {L.shape=}" + assert np.allclose(np.tril(L), L), "L is not lower triangular" + # Handle vector case by reshaping to matrix and then flattening again - if len(y.shape) == 1: - return solve_cholesky(lo, y.reshape(-1, 1)).ravel() + if len(Y.shape) == 1: + return solve_cholesky(L, Y.reshape(-1, 1)).ravel() - n, m = y.shape + n = Y.shape[0] - # Backsubstitute L X = B - x = y.copy() + # Solve L W = B for W + W = Y.copy() for i in range(n): for j in range(i): - x[i, :] -= lo[i, j] * x[j, :] + W[i] -= L[i, j] * W[j] - for k in range(m): - x[i, k] /= lo[i, i] + W[i] /= L[i, i] - # Backsubstitute L^T + # Solve L^T X = W for X + X = W for i in reversed(range(n)): for j in range(i + 1, n): - x[i, :] -= lo[j, i] * x[j, :] + X[i] -= L[j, i] * X[j] - for k in range(m): - x[i, k] /= lo[i, i] + X[i] /= L[i, i] - return x + return X if __name__ == "__main__": From 907b783668290ea76b08539bd811c721db2a1654 Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Wed, 9 Oct 2024 08:42:07 +0200 Subject: [PATCH 3/6] Enforce symmetry on A --- maths/cholesky_decomposition.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py index 363436b40ac3..a17a5c2a702f 100644 --- a/maths/cholesky_decomposition.py +++ b/maths/cholesky_decomposition.py @@ -43,7 +43,8 @@ def cholesky_decomposition(A: np.ndarray) -> np.ndarray: True """ - assert A.shape[0] == A.shape[1], f"A is not square, {A.shape=}" + assert A.shape[0] == A.shape[1], f"Matrix A is not square, {A.shape=}" + assert np.allclose(A, A.T), "Matrix A must be symmetric" n = A.shape[0] L = np.tril(A) @@ -74,8 +75,8 @@ def solve_cholesky(L: np.ndarray, Y: np.ndarray) -> np.ndarray: True """ - assert L.shape[0] == L.shape[1], f"L is not square, {L.shape=}" - assert np.allclose(np.tril(L), L), "L is not lower triangular" + assert L.shape[0] == L.shape[1], f"Matrix L is not square, {L.shape=}" + assert np.allclose(np.tril(L), L), "Matrix L is not lower triangular" # Handle vector case by reshaping to matrix and then flattening again if len(Y.shape) == 1: From 818448b05dd51f36c07eba15f87cc77ed55f375c Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Wed, 9 Oct 2024 09:09:27 +0200 Subject: [PATCH 4/6] Fix typo --- maths/cholesky_decomposition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py index a17a5c2a702f..d7fa5b45cba7 100644 --- a/maths/cholesky_decomposition.py +++ b/maths/cholesky_decomposition.py @@ -66,7 +66,7 @@ def cholesky_decomposition(A: np.ndarray) -> np.ndarray: def solve_cholesky(L: np.ndarray, Y: np.ndarray) -> np.ndarray: """Given a Cholesky decomposition L L^T = A of a matrix A, solve the - system of equations A X = Y where B is either a matrix or a vector. + system of equations A X = Y where Y is either a matrix or a vector. >>> L = np.array([[2, 0], [3, 4]], dtype=float) >>> Y = np.array([[22, 54], [81, 193]], dtype=float) From 45222589804ce2b7fd055d5d8e30cd7fc0c207e3 Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Thu, 10 Oct 2024 08:00:18 +0200 Subject: [PATCH 5/6] Rename variables --- maths/cholesky_decomposition.py | 55 ++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py index d7fa5b45cba7..db3b47511f40 100644 --- a/maths/cholesky_decomposition.py +++ b/maths/cholesky_decomposition.py @@ -1,8 +1,7 @@ import numpy as np -# ruff: noqa: N803,N806 -def cholesky_decomposition(A: np.ndarray) -> np.ndarray: +def cholesky_decomposition(a: np.ndarray) -> np.ndarray: """Return a Cholesky decomposition of the matrix A. The Cholesky decomposition decomposes the square, positive definite matrix A @@ -26,7 +25,7 @@ def cholesky_decomposition(A: np.ndarray) -> np.ndarray: >>> np.allclose(np.tril(L), L) True - The Cholesky decomposition can be used to solve the system of equations A x = y. + The Cholesky decomposition can be used to solve the linear system A x = y. >>> x_true = np.array([1, 2, 3], dtype=float) >>> y = A @ x_true @@ -43,28 +42,30 @@ def cholesky_decomposition(A: np.ndarray) -> np.ndarray: True """ - assert A.shape[0] == A.shape[1], f"Matrix A is not square, {A.shape=}" - assert np.allclose(A, A.T), "Matrix A must be symmetric" + assert a.shape[0] == a.shape[1], f"Matrix A is not square, {a.shape=}" + assert np.allclose(a, a.T), "Matrix A must be symmetric" - n = A.shape[0] - L = np.tril(A) + n = a.shape[0] + lower_triangle = np.tril(a) for i in range(n): for j in range(i + 1): - L[i, j] -= np.sum(L[i, :j] * L[j, :j]) + lower_triangle[i, j] -= np.sum( + lower_triangle[i, :j] * lower_triangle[j, :j] + ) if i == j: - if L[i, i] <= 0: + if lower_triangle[i, i] <= 0: raise ValueError("Matrix A is not positive definite") - L[i, i] = np.sqrt(L[i, i]) + lower_triangle[i, i] = np.sqrt(lower_triangle[i, i]) else: - L[i, j] /= L[j, j] + lower_triangle[i, j] /= lower_triangle[j, j] - return L + return lower_triangle -def solve_cholesky(L: np.ndarray, Y: np.ndarray) -> np.ndarray: +def solve_cholesky(lower_triangle: np.ndarray, y: np.ndarray) -> np.ndarray: """Given a Cholesky decomposition L L^T = A of a matrix A, solve the system of equations A X = Y where Y is either a matrix or a vector. @@ -75,32 +76,36 @@ def solve_cholesky(L: np.ndarray, Y: np.ndarray) -> np.ndarray: True """ - assert L.shape[0] == L.shape[1], f"Matrix L is not square, {L.shape=}" - assert np.allclose(np.tril(L), L), "Matrix L is not lower triangular" + assert ( + lower_triangle.shape[0] == lower_triangle.shape[1] + ), f"Matrix L is not square, {lower_triangle.shape=}" + assert np.allclose( + np.tril(lower_triangle), lower_triangle + ), "Matrix L is not lower triangular" # Handle vector case by reshaping to matrix and then flattening again - if len(Y.shape) == 1: - return solve_cholesky(L, Y.reshape(-1, 1)).ravel() + if len(y.shape) == 1: + return solve_cholesky(lower_triangle, y.reshape(-1, 1)).ravel() - n = Y.shape[0] + n = y.shape[0] # Solve L W = B for W - W = Y.copy() + w = y.copy() for i in range(n): for j in range(i): - W[i] -= L[i, j] * W[j] + w[i] -= lower_triangle[i, j] * w[j] - W[i] /= L[i, i] + w[i] /= lower_triangle[i, i] # Solve L^T X = W for X - X = W + x = w for i in reversed(range(n)): for j in range(i + 1, n): - X[i] -= L[j, i] * X[j] + x[i] -= lower_triangle[j, i] * x[j] - X[i] /= L[i, i] + x[i] /= lower_triangle[i, i] - return X + return x if __name__ == "__main__": From 7808f21e1cda505088712c0dd8cf50195a593a22 Mon Sep 17 00:00:00 2001 From: 99991 <99991@users.noreply.github.com> Date: Thu, 10 Oct 2024 08:46:05 +0200 Subject: [PATCH 6/6] Rename variables --- maths/cholesky_decomposition.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/maths/cholesky_decomposition.py b/maths/cholesky_decomposition.py index db3b47511f40..591a3e9dffd9 100644 --- a/maths/cholesky_decomposition.py +++ b/maths/cholesky_decomposition.py @@ -1,7 +1,7 @@ import numpy as np -def cholesky_decomposition(a: np.ndarray) -> np.ndarray: +def cholesky_decomposition(matrix: np.ndarray) -> np.ndarray: """Return a Cholesky decomposition of the matrix A. The Cholesky decomposition decomposes the square, positive definite matrix A @@ -42,11 +42,13 @@ def cholesky_decomposition(a: np.ndarray) -> np.ndarray: True """ - assert a.shape[0] == a.shape[1], f"Matrix A is not square, {a.shape=}" - assert np.allclose(a, a.T), "Matrix A must be symmetric" + assert ( + matrix.shape[0] == matrix.shape[1] + ), f"Input matrix is not square, {matrix.shape=}" + assert np.allclose(matrix, matrix.T), "Input matrix must be symmetric" - n = a.shape[0] - lower_triangle = np.tril(a) + n = matrix.shape[0] + lower_triangle = np.tril(matrix) for i in range(n): for j in range(i + 1): @@ -65,9 +67,13 @@ def cholesky_decomposition(a: np.ndarray) -> np.ndarray: return lower_triangle -def solve_cholesky(lower_triangle: np.ndarray, y: np.ndarray) -> np.ndarray: +def solve_cholesky( + lower_triangle: np.ndarray, + right_hand_side: np.ndarray, +) -> np.ndarray: """Given a Cholesky decomposition L L^T = A of a matrix A, solve the - system of equations A X = Y where Y is either a matrix or a vector. + system of equations A X = Y where the right-hand side Y is either + a matrix or a vector. >>> L = np.array([[2, 0], [3, 4]], dtype=float) >>> Y = np.array([[22, 54], [81, 193]], dtype=float) @@ -84,13 +90,13 @@ def solve_cholesky(lower_triangle: np.ndarray, y: np.ndarray) -> np.ndarray: ), "Matrix L is not lower triangular" # Handle vector case by reshaping to matrix and then flattening again - if len(y.shape) == 1: - return solve_cholesky(lower_triangle, y.reshape(-1, 1)).ravel() + if len(right_hand_side.shape) == 1: + return solve_cholesky(lower_triangle, right_hand_side.reshape(-1, 1)).ravel() - n = y.shape[0] + n = right_hand_side.shape[0] - # Solve L W = B for W - w = y.copy() + # Solve L W = Y for W + w = right_hand_side.copy() for i in range(n): for j in range(i): w[i] -= lower_triangle[i, j] * w[j]