From 78811cc29a31c633279ec064183ba93e54c18dc4 Mon Sep 17 00:00:00 2001 From: Ivan Radonov Date: Mon, 27 Sep 2021 17:54:54 +0300 Subject: [PATCH 1/5] added schur complement and tests to linear algebra --- linear_algebra/src/schur_complement.py | 72 ++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 linear_algebra/src/schur_complement.py diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py new file mode 100644 index 000000000000..e483f79d21f8 --- /dev/null +++ b/linear_algebra/src/schur_complement.py @@ -0,0 +1,72 @@ +import numpy as np + + +def schur_complement( + A: np.ndarray, B: np.ndarray, C: np.ndarray, pseudo_inv: np.ndarray = None +) -> np.ndarray: + """ + Schur complement of a symmetric matrix X given as a 2x2 block matrix + consisting of matrices A, B and C. + Matrix A must be quadratic and non-singular. + In case A is singular, a pseudo-inverse may be provided using + the pseudo_inv argument. + >>> import numpy as np + >>> A = np.array([[1, 2], [2, 1]]) + >>> B = np.array([[0, 3], [3, 0]]) + >>> C = np.array([[2, 1], [6, 3]]) + >>> schur_complement(A, B, C) + array([[ 5., -5.], + [ 0., 6.]]) + """ + shape_A = np.shape(A) + shape_B = np.shape(B) + shape_C = np.shape(C) + + if shape_A[0] != shape_B[0]: + raise ValueError( + f"Expected the same number of rows for A and B. \ + Instead found A of size {shape_A} and B of size {shape_B}" + ) + + if shape_B[1] != shape_C[1]: + raise ValueError( + f"Expected the same number of columns for B and C. \ + Instead found B of size {shape_B} and C of size {shape_C}" + ) + + A_inv = pseudo_inv + if A_inv is None: + try: + A_inv = np.linalg.inv(A) + except np.linalg.LinAlgError: + raise ValueError( + "Input matrix A is not invertible. Cannot compute Schur complement." + ) + + return C - B.T @ A_inv @ B + + +def test_schur_complement(): + """ + >>> test_schur_complement() # self running tests + """ + A = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) + B = np.array([[0, 3], [3, 0], [2, 3]]) + C = np.array([[2, 1], [6, 3]]) + + S = schur_complement(A, B, C) + + input_matrix = np.block([[A, B], [B.T, C]]) + + det_X = np.linalg.det(input_matrix) + det_A = np.linalg.det(A) + det_S = np.linalg.det(S) + + assert np.abs(det_X - det_A * det_S) <= 1e-6 + + +if __name__ == "__main__": + import doctest + + doctest.testmod() + test_schur_complement() From 8498cb73b34ccf555589a78529dc8958050ac65d Mon Sep 17 00:00:00 2001 From: Ivan Radonov Date: Mon, 27 Sep 2021 18:16:27 +0300 Subject: [PATCH 2/5] updated according to checklist --- linear_algebra/src/schur_complement.py | 53 ++++++++++++++------------ 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py index e483f79d21f8..b4ee7f2af1c7 100644 --- a/linear_algebra/src/schur_complement.py +++ b/linear_algebra/src/schur_complement.py @@ -2,7 +2,7 @@ def schur_complement( - A: np.ndarray, B: np.ndarray, C: np.ndarray, pseudo_inv: np.ndarray = None + a: np.ndarray, b: np.ndarray, c: np.ndarray, pseudo_inv: np.ndarray = None ) -> np.ndarray: """ Schur complement of a symmetric matrix X given as a 2x2 block matrix @@ -10,59 +10,62 @@ def schur_complement( Matrix A must be quadratic and non-singular. In case A is singular, a pseudo-inverse may be provided using the pseudo_inv argument. + + Link to Wiki: https://en.wikipedia.org/wiki/Schur_complement + See also Convex Optimization – Boyd and Vandenberghe, A.5.5 >>> import numpy as np - >>> A = np.array([[1, 2], [2, 1]]) - >>> B = np.array([[0, 3], [3, 0]]) - >>> C = np.array([[2, 1], [6, 3]]) - >>> schur_complement(A, B, C) + >>> a = np.array([[1, 2], [2, 1]]) + >>> b = np.array([[0, 3], [3, 0]]) + >>> c = np.array([[2, 1], [6, 3]]) + >>> schur_complement(a, b, c) array([[ 5., -5.], [ 0., 6.]]) """ - shape_A = np.shape(A) - shape_B = np.shape(B) - shape_C = np.shape(C) + shape_a = np.shape(a) + shape_b = np.shape(b) + shape_c = np.shape(c) - if shape_A[0] != shape_B[0]: + if shape_a[0] != shape_b[0]: raise ValueError( f"Expected the same number of rows for A and B. \ - Instead found A of size {shape_A} and B of size {shape_B}" + Instead found A of size {shape_a} and B of size {shape_b}" ) - if shape_B[1] != shape_C[1]: + if shape_b[1] != shape_c[1]: raise ValueError( f"Expected the same number of columns for B and C. \ - Instead found B of size {shape_B} and C of size {shape_C}" + Instead found B of size {shape_b} and C of size {shape_c}" ) - A_inv = pseudo_inv - if A_inv is None: + a_inv = pseudo_inv + if a_inv is None: try: - A_inv = np.linalg.inv(A) + a_inv = np.linalg.inv(a) except np.linalg.LinAlgError: raise ValueError( "Input matrix A is not invertible. Cannot compute Schur complement." ) - return C - B.T @ A_inv @ B + return c - b.T @ a_inv @ b def test_schur_complement(): """ >>> test_schur_complement() # self running tests """ - A = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) - B = np.array([[0, 3], [3, 0], [2, 3]]) - C = np.array([[2, 1], [6, 3]]) + a = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) + b = np.array([[0, 3], [3, 0], [2, 3]]) + c = np.array([[2, 1], [6, 3]]) - S = schur_complement(A, B, C) + s = schur_complement(a, b, c) - input_matrix = np.block([[A, B], [B.T, C]]) + input_matrix = np.block([[a, b], [b.T, c]]) - det_X = np.linalg.det(input_matrix) - det_A = np.linalg.det(A) - det_S = np.linalg.det(S) + det_x = np.linalg.det(input_matrix) + det_a = np.linalg.det(a) + det_s = np.linalg.det(s) - assert np.abs(det_X - det_A * det_S) <= 1e-6 + assert np.abs(det_x - det_a * det_s) <= 1e-6 if __name__ == "__main__": From ba3d8727a4bb8a9b6fcb135aed79628fc86e2f74 Mon Sep 17 00:00:00 2001 From: Ivan Radonov Date: Mon, 27 Sep 2021 18:24:23 +0300 Subject: [PATCH 3/5] updated variable names and typing --- linear_algebra/src/schur_complement.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py index b4ee7f2af1c7..7a2b9c944e4b 100644 --- a/linear_algebra/src/schur_complement.py +++ b/linear_algebra/src/schur_complement.py @@ -2,7 +2,10 @@ def schur_complement( - a: np.ndarray, b: np.ndarray, c: np.ndarray, pseudo_inv: np.ndarray = None + mat_a: np.ndarray, + mat_b: np.ndarray, + mat_c: np.ndarray, + pseudo_inv: np.ndarray = None, ) -> np.ndarray: """ Schur complement of a symmetric matrix X given as a 2x2 block matrix @@ -21,9 +24,9 @@ def schur_complement( array([[ 5., -5.], [ 0., 6.]]) """ - shape_a = np.shape(a) - shape_b = np.shape(b) - shape_c = np.shape(c) + shape_a = np.shape(mat_a) + shape_b = np.shape(mat_b) + shape_c = np.shape(mat_c) if shape_a[0] != shape_b[0]: raise ValueError( @@ -40,16 +43,16 @@ def schur_complement( a_inv = pseudo_inv if a_inv is None: try: - a_inv = np.linalg.inv(a) + a_inv = np.linalg.inv(mat_a) except np.linalg.LinAlgError: raise ValueError( "Input matrix A is not invertible. Cannot compute Schur complement." ) - return c - b.T @ a_inv @ b + return mat_c - mat_b.T @ a_inv @ mat_b -def test_schur_complement(): +def test_schur_complement() -> None: """ >>> test_schur_complement() # self running tests """ From bcffee3eb15c6389d47b296416fdf7ec8764cd51 Mon Sep 17 00:00:00 2001 From: Ivan Radonov Date: Thu, 30 Sep 2021 12:41:00 +0300 Subject: [PATCH 4/5] added two testcases for input validation --- linear_algebra/src/schur_complement.py | 43 +++++++++++++++++--------- 1 file changed, 29 insertions(+), 14 deletions(-) diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py index 7a2b9c944e4b..33e2acde9727 100644 --- a/linear_algebra/src/schur_complement.py +++ b/linear_algebra/src/schur_complement.py @@ -1,4 +1,5 @@ import numpy as np +import unittest def schur_complement( @@ -52,27 +53,41 @@ def schur_complement( return mat_c - mat_b.T @ a_inv @ mat_b -def test_schur_complement() -> None: - """ - >>> test_schur_complement() # self running tests - """ - a = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) - b = np.array([[0, 3], [3, 0], [2, 3]]) - c = np.array([[2, 1], [6, 3]]) +class TestSchurComplement(unittest.TestCase): + def test_schur_complement(self) -> None: + a = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) + b = np.array([[0, 3], [3, 0], [2, 3]]) + c = np.array([[2, 1], [6, 3]]) + + s = schur_complement(a, b, c) + + input_matrix = np.block([[a, b], [b.T, c]]) + + det_x = np.linalg.det(input_matrix) + det_a = np.linalg.det(a) + det_s = np.linalg.det(s) + + self.assertAlmostEqual(det_x, det_a * det_s) - s = schur_complement(a, b, c) + def test_improper_a_b_dimensions(self) -> None: + a = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) + b = np.array([[0, 3], [3, 0], [2, 3]]) + c = np.array([[2, 1], [6, 3]]) - input_matrix = np.block([[a, b], [b.T, c]]) + with self.assertRaises(ValueError): + schur_complement(a, b, c) - det_x = np.linalg.det(input_matrix) - det_a = np.linalg.det(a) - det_s = np.linalg.det(s) + def test_improper_b_c_dimensions(self) -> None: + a = np.array([[1, 2, 1], [2, 1, 2], [3, 2, 4]]) + b = np.array([[0, 3], [3, 0], [2, 3]]) + c = np.array([[2, 1, 3], [6, 3, 5]]) - assert np.abs(det_x - det_a * det_s) <= 1e-6 + with self.assertRaises(ValueError): + schur_complement(a, b, c) if __name__ == "__main__": import doctest doctest.testmod() - test_schur_complement() + unittest.main() From 7f975b238704fab34303747dbe531ebe914b016e Mon Sep 17 00:00:00 2001 From: Ivan Radonov Date: Wed, 13 Oct 2021 15:20:55 +0300 Subject: [PATCH 5/5] fixed import order --- linear_algebra/src/schur_complement.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/linear_algebra/src/schur_complement.py b/linear_algebra/src/schur_complement.py index 33e2acde9727..f3cb736d9084 100644 --- a/linear_algebra/src/schur_complement.py +++ b/linear_algebra/src/schur_complement.py @@ -1,6 +1,7 @@ -import numpy as np import unittest +import numpy as np + def schur_complement( mat_a: np.ndarray,