From bf4fcef68caaf0c0309287c041439b7f92855383 Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Tue, 16 Jun 2020 08:11:21 -0600 Subject: [PATCH 1/9] Add check on eigenvalues to determine semidefiniteness --- sympy/matrices/eigen.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 48b9317b60cd..0235ebaff2f0 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -5,6 +5,7 @@ from mpmath import mp, workprec from mpmath.libmp.libmpf import prec_to_dps +from sympy import Symbol from sympy.core.compatibility import default_sort_key from sympy.core.evalf import DEFAULT_MAXPREC, PrecisionExhausted from sympy.core.logic import fuzzy_and, fuzzy_or @@ -700,6 +701,12 @@ def _is_positive_semidefinite(M): if nonnegative_diagonals and M.is_weakly_diagonally_dominant: return True + if M.rows < 5 or all(a.func != Symbol for a in M.atoms()): + try: + return _is_positive_semidefinite_by_eigenvalues(M) + except Exception: + pass + return _is_positive_semidefinite_by_minors(M) @@ -752,6 +759,12 @@ def _is_positive_semidefinite_by_minors(M): ) +def _is_positive_semidefinite_by_eigenvalues(M): + """Determines if a matrix is positive semidefinite by checking + if all the eigenvalues are nonnegative""" + return all(eigenvalue.is_nonnegative for eigenvalue in M.eigenvals()) + + _doc_positive_definite = \ r"""Finds out the definiteness of a matrix. From 1bf89e8f36b5ef9b9037e4ad7a9c6625cd561edd Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Tue, 16 Jun 2020 23:34:32 -0600 Subject: [PATCH 2/9] Got cholesky algorithm to pass existing tests --- sympy/matrices/eigen.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 0235ebaff2f0..fef3a9dc2599 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -701,6 +701,11 @@ def _is_positive_semidefinite(M): if nonnegative_diagonals and M.is_weakly_diagonally_dominant: return True + try: + return _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M) + except Exception: + pass + if M.rows < 5 or all(a.func != Symbol for a in M.atoms()): try: return _is_positive_semidefinite_by_eigenvalues(M) @@ -765,6 +770,25 @@ def _is_positive_semidefinite_by_eigenvalues(M): return all(eigenvalue.is_nonnegative for eigenvalue in M.eigenvals()) +def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): + M = M.copy() + for k in range(M.rows): + pivot = max(range(k, M.rows), key=lambda i: M[i, i]) + if pivot > k: + M[[k, pivot], :] = M[[pivot, k], :] + M[:, [k, pivot]] = M[:, [pivot, k]] + if M[k, k].is_negative: + return False + M[k, k] = sqrt(M[k, k]) + for j in range(k+1, M.rows): + M[k, j] /= M[k, k] + for j in range(k+1, M.rows): + for i in range(k+1, j+1): + M[i, j] -= M[k, i] * M[k, j] + return M[-1, -1].is_nonnegative + + + _doc_positive_definite = \ r"""Finds out the definiteness of a matrix. From 50fbb841fd3539d8bb2fbf6f3ee6d2eb96e050d9 Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Tue, 16 Jun 2020 23:38:23 -0600 Subject: [PATCH 3/9] vectorize inner loops --- sympy/matrices/eigen.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index fef3a9dc2599..5479c568c510 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -780,11 +780,9 @@ def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): if M[k, k].is_negative: return False M[k, k] = sqrt(M[k, k]) + M[k, (k+1):] /= M[k, k] for j in range(k+1, M.rows): - M[k, j] /= M[k, k] - for j in range(k+1, M.rows): - for i in range(k+1, j+1): - M[i, j] -= M[k, i] * M[k, j] + M[(k+1):(j+1), j] = M[k, (k+1):(j+1)].T * M[k, j] return M[-1, -1].is_nonnegative From b07e43fc261b453f57d718cfb565d0a88eac9880 Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Tue, 16 Jun 2020 23:47:59 -0600 Subject: [PATCH 4/9] Add reference --- sympy/matrices/eigen.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 5479c568c510..246ac5505ed5 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -771,6 +771,24 @@ def _is_positive_semidefinite_by_eigenvalues(M): def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): + """Determines if a matrix is positive semidefinite by computing its + Cholesky decomposition with complete pivoting + + Any symmetric positive semidefinite matrix has a factorization + P.T * A * P = R.T * R, where + + R = [R_11 R_22] + [0 0 ], + + and R_11 is rxr upper triangular with positive diagonal elements, + rank(A) = r, and P is a permutation matrix. + + An algorithm for computing R is described in [1]. If at any point in + the algorithm a negative diagonal term is obtained, then A is not + positive semidefinite, so the algorithm exits early and returns False. + + [1] http://eprints.ma.man.ac.uk/1199/1/covered/MIMS_ep2008_116.pdf + """ M = M.copy() for k in range(M.rows): pivot = max(range(k, M.rows), key=lambda i: M[i, i]) @@ -786,7 +804,6 @@ def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): return M[-1, -1].is_nonnegative - _doc_positive_definite = \ r"""Finds out the definiteness of a matrix. From 9e80ddf82f9a94f895c1bd244b4f90df16bbc3cf Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Wed, 17 Jun 2020 00:32:40 -0600 Subject: [PATCH 5/9] Add individual tests for subroutines --- sympy/matrices/eigen.py | 16 ++++++------ sympy/matrices/tests/test_eigen.py | 40 ++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 8 deletions(-) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 246ac5505ed5..90ed5946076c 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -793,14 +793,14 @@ def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): for k in range(M.rows): pivot = max(range(k, M.rows), key=lambda i: M[i, i]) if pivot > k: - M[[k, pivot], :] = M[[pivot, k], :] - M[:, [k, pivot]] = M[:, [pivot, k]] - if M[k, k].is_negative: - return False - M[k, k] = sqrt(M[k, k]) - M[k, (k+1):] /= M[k, k] - for j in range(k+1, M.rows): - M[(k+1):(j+1), j] = M[k, (k+1):(j+1)].T * M[k, j] + M.col_swap(k, pivot) + M.row_swap(k, pivot) + if M[k, k].is_negative: + return False + M[k, k] = sqrt(M[k, k]) + M[k, (k+1):] /= M[k, k] + for j in range(k+1, M.rows): + M[(k+1):(j+1), j] -= M[k, (k+1):(j+1)].T * M[k, j] return M[-1, -1].is_nonnegative diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py index 5feaa5a733e3..e878cef8eb48 100644 --- a/sympy/matrices/tests/test_eigen.py +++ b/sympy/matrices/tests/test_eigen.py @@ -4,6 +4,11 @@ from sympy.matrices import eye, Matrix from sympy.core.singleton import S from sympy.testing.pytest import raises, XFAIL +from sympy.matrices.eigen import ( + _is_positive_semidefinite_by_eigenvalues, + _is_positive_semidefinite_by_cholesky_factorization_with_pivots, + _is_positive_semidefinite_by_minors +) from sympy.matrices.matrices import NonSquareMatrixError, MatrixError from sympy.simplify.simplify import simplify from sympy.matrices.immutable import ImmutableMatrix @@ -493,6 +498,9 @@ def test_definite(): m = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert _is_positive_semidefinite_by_eigenvalues(m) + assert _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -500,6 +508,9 @@ def test_definite(): m = Matrix([[5, 4], [4, 5]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert _is_positive_semidefinite_by_eigenvalues(m) + assert _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -508,6 +519,9 @@ def test_definite(): m = Matrix([[2, -1, -1], [-1, 2, -1], [-1, -1, 2]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert _is_positive_semidefinite_by_eigenvalues(m) + assert _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -515,6 +529,9 @@ def test_definite(): m = Matrix([[1, 2], [2, 4]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) + assert _is_positive_semidefinite_by_eigenvalues(m.T + m) + assert _is_positive_semidefinite_by_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -524,13 +541,20 @@ def test_definite(): m = Matrix([[2, 3], [4, 8]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) + assert _is_positive_semidefinite_by_eigenvalues(m.T + m) + assert _is_positive_semidefinite_by_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False + # Hermetian matrices m = Matrix([[1, 2*I], [-I, 4]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.H + m) + assert _is_positive_semidefinite_by_eigenvalues(m.H + m) + assert _is_positive_semidefinite_by_minors(m.H + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -541,6 +565,9 @@ def test_definite(): m = Matrix([[a, 0, 0], [0, a, 0], [0, 0, a]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert _is_positive_semidefinite_by_eigenvalues(m) + assert _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -548,6 +575,9 @@ def test_definite(): m = Matrix([[b, 0, 0], [0, b, 0], [0, 0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False + assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert not _is_positive_semidefinite_by_eigenvalues(m) + assert not _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == True assert m.is_negative_semidefinite == True assert m.is_indefinite == False @@ -555,6 +585,9 @@ def test_definite(): m = Matrix([[a, 0], [0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False + assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert not _is_positive_semidefinite_by_eigenvalues(m) + assert not _is_positive_semidefinite_by_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == True @@ -571,6 +604,9 @@ def test_definite(): ]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True + assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) + assert _is_positive_semidefinite_by_eigenvalues(m.T + m) + assert _is_positive_semidefinite_by_minors(m.T + m) assert m.is_indefinite == False # test for issue 19547: https://github.com/sympy/sympy/issues/19547 @@ -581,3 +617,7 @@ def test_definite(): ]) assert not m.is_positive_definite assert not m.is_positive_semidefinite + assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) + assert not _is_positive_semidefinite_by_eigenvalues(m) + assert not _is_positive_semidefinite_by_minors(m) + From 1682ac3e710e950c0374a53839b75cda9f2682c8 Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Wed, 17 Jun 2020 07:58:32 -0600 Subject: [PATCH 6/9] Use shorter function names --- sympy/matrices/eigen.py | 12 ++--- sympy/matrices/tests/test_eigen.py | 72 +++++++++++++++--------------- 2 files changed, 42 insertions(+), 42 deletions(-) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 90ed5946076c..78c1dbf61cd8 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -702,17 +702,17 @@ def _is_positive_semidefinite(M): return True try: - return _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M) + return _is_positive_semidefinite_cholesky(M) except Exception: pass if M.rows < 5 or all(a.func != Symbol for a in M.atoms()): try: - return _is_positive_semidefinite_by_eigenvalues(M) + return _is_positive_semidefinite_evals(M) except Exception: pass - return _is_positive_semidefinite_by_minors(M) + return _is_positive_semidefinite_minors(M) def _is_negative_definite(M): @@ -754,7 +754,7 @@ def _is_positive_definite_GE(M): return True -def _is_positive_semidefinite_by_minors(M): +def _is_positive_semidefinite_minors(M): """A method to evaluate all principal minors for testing positive-semidefiniteness by using Sylvestre's crieterion.""" return all( @@ -764,13 +764,13 @@ def _is_positive_semidefinite_by_minors(M): ) -def _is_positive_semidefinite_by_eigenvalues(M): +def _is_positive_semidefinite_evals(M): """Determines if a matrix is positive semidefinite by checking if all the eigenvalues are nonnegative""" return all(eigenvalue.is_nonnegative for eigenvalue in M.eigenvals()) -def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M): +def _is_positive_semidefinite_cholesky(M): """Determines if a matrix is positive semidefinite by computing its Cholesky decomposition with complete pivoting diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py index e878cef8eb48..d8521115bd52 100644 --- a/sympy/matrices/tests/test_eigen.py +++ b/sympy/matrices/tests/test_eigen.py @@ -5,9 +5,9 @@ from sympy.core.singleton import S from sympy.testing.pytest import raises, XFAIL from sympy.matrices.eigen import ( - _is_positive_semidefinite_by_eigenvalues, - _is_positive_semidefinite_by_cholesky_factorization_with_pivots, - _is_positive_semidefinite_by_minors + _is_positive_semidefinite_cholesky, + _is_positive_semidefinite_evals, + _is_positive_semidefinite_minors ) from sympy.matrices.matrices import NonSquareMatrixError, MatrixError from sympy.simplify.simplify import simplify @@ -498,9 +498,9 @@ def test_definite(): m = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert _is_positive_semidefinite_by_eigenvalues(m) - assert _is_positive_semidefinite_by_minors(m) + assert _is_positive_semidefinite_cholesky(m) + assert _is_positive_semidefinite_evals(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -508,9 +508,9 @@ def test_definite(): m = Matrix([[5, 4], [4, 5]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert _is_positive_semidefinite_by_eigenvalues(m) - assert _is_positive_semidefinite_by_minors(m) + assert _is_positive_semidefinite_cholesky(m) + assert _is_positive_semidefinite_evals(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -519,9 +519,9 @@ def test_definite(): m = Matrix([[2, -1, -1], [-1, 2, -1], [-1, -1, 2]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert _is_positive_semidefinite_by_eigenvalues(m) - assert _is_positive_semidefinite_by_minors(m) + assert _is_positive_semidefinite_cholesky(m) + assert _is_positive_semidefinite_evals(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -529,9 +529,9 @@ def test_definite(): m = Matrix([[1, 2], [2, 4]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) - assert _is_positive_semidefinite_by_eigenvalues(m.T + m) - assert _is_positive_semidefinite_by_minors(m.T + m) + assert _is_positive_semidefinite_cholesky(m.T + m) + assert _is_positive_semidefinite_evals(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -541,9 +541,9 @@ def test_definite(): m = Matrix([[2, 3], [4, 8]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) - assert _is_positive_semidefinite_by_eigenvalues(m.T + m) - assert _is_positive_semidefinite_by_minors(m.T + m) + assert _is_positive_semidefinite_cholesky(m.T + m) + assert _is_positive_semidefinite_evals(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -552,9 +552,9 @@ def test_definite(): m = Matrix([[1, 2*I], [-I, 4]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.H + m) - assert _is_positive_semidefinite_by_eigenvalues(m.H + m) - assert _is_positive_semidefinite_by_minors(m.H + m) + assert _is_positive_semidefinite_cholesky(m.H + m) + assert _is_positive_semidefinite_evals(m.H + m) + assert _is_positive_semidefinite_minors(m.H + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -565,9 +565,9 @@ def test_definite(): m = Matrix([[a, 0, 0], [0, a, 0], [0, 0, a]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert _is_positive_semidefinite_by_eigenvalues(m) - assert _is_positive_semidefinite_by_minors(m) + assert _is_positive_semidefinite_cholesky(m) + assert _is_positive_semidefinite_evals(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -575,9 +575,9 @@ def test_definite(): m = Matrix([[b, 0, 0], [0, b, 0], [0, 0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False - assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert not _is_positive_semidefinite_by_eigenvalues(m) - assert not _is_positive_semidefinite_by_minors(m) + assert not _is_positive_semidefinite_cholesky(m) + assert not _is_positive_semidefinite_evals(m) + assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == True assert m.is_negative_semidefinite == True assert m.is_indefinite == False @@ -585,9 +585,9 @@ def test_definite(): m = Matrix([[a, 0], [0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False - assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert not _is_positive_semidefinite_by_eigenvalues(m) - assert not _is_positive_semidefinite_by_minors(m) + assert not _is_positive_semidefinite_cholesky(m) + assert not _is_positive_semidefinite_evals(m) + assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == True @@ -604,9 +604,9 @@ def test_definite(): ]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m.T + m) - assert _is_positive_semidefinite_by_eigenvalues(m.T + m) - assert _is_positive_semidefinite_by_minors(m.T + m) + assert _is_positive_semidefinite_cholesky(m.T + m) + assert _is_positive_semidefinite_evals(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_indefinite == False # test for issue 19547: https://github.com/sympy/sympy/issues/19547 @@ -617,7 +617,7 @@ def test_definite(): ]) assert not m.is_positive_definite assert not m.is_positive_semidefinite - assert not _is_positive_semidefinite_by_cholesky_factorization_with_pivots(m) - assert not _is_positive_semidefinite_by_eigenvalues(m) - assert not _is_positive_semidefinite_by_minors(m) + assert not _is_positive_semidefinite_cholesky(m) + assert not _is_positive_semidefinite_evals(m) + assert not _is_positive_semidefinite_minors(m) From 7bb906e6afea33585ff307bdfc093c0c9af40f57 Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Wed, 17 Jun 2020 08:13:20 -0600 Subject: [PATCH 7/9] Remove trailing whitespace and update flake8 config --- setup.cfg | 2 +- sympy/matrices/tests/test_eigen.py | 23 +++++++++++------------ 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/setup.cfg b/setup.cfg index 112c744bb447..1093aea0d65d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -5,7 +5,7 @@ long_description_content_type = text/markdown [flake8] doctests = True ignore = F403 -select = F,E722 +select = F,E722,W exclude = sympy/parsing/latex/_antlr/*, sympy/parsing/autolev/_antlr/*, sympy/parsing/autolev/test-examples/*, diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py index d8521115bd52..412b69422bcb 100644 --- a/sympy/matrices/tests/test_eigen.py +++ b/sympy/matrices/tests/test_eigen.py @@ -500,7 +500,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m) assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -510,7 +510,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m) assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -521,7 +521,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m) assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -531,7 +531,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m.T + m) assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -543,7 +543,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m.T + m) assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -554,7 +554,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m.H + m) assert _is_positive_semidefinite_evals(m.H + m) - assert _is_positive_semidefinite_minors(m.H + m) + assert _is_positive_semidefinite_minors(m.H + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -567,7 +567,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m) assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) + assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -577,7 +577,7 @@ def test_definite(): assert m.is_positive_semidefinite == False assert not _is_positive_semidefinite_cholesky(m) assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m) + assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == True assert m.is_negative_semidefinite == True assert m.is_indefinite == False @@ -587,7 +587,7 @@ def test_definite(): assert m.is_positive_semidefinite == False assert not _is_positive_semidefinite_cholesky(m) assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m) + assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == True @@ -606,7 +606,7 @@ def test_definite(): assert m.is_positive_semidefinite == True assert _is_positive_semidefinite_cholesky(m.T + m) assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) + assert _is_positive_semidefinite_minors(m.T + m) assert m.is_indefinite == False # test for issue 19547: https://github.com/sympy/sympy/issues/19547 @@ -619,5 +619,4 @@ def test_definite(): assert not m.is_positive_semidefinite assert not _is_positive_semidefinite_cholesky(m) assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m) - + assert not _is_positive_semidefinite_minors(m) From b15fcf096de9d4ccbf2df0d8e43aaa6bcc22d70c Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Wed, 17 Jun 2020 08:57:33 -0600 Subject: [PATCH 8/9] Revert setup.cfg --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 1093aea0d65d..112c744bb447 100644 --- a/setup.cfg +++ b/setup.cfg @@ -5,7 +5,7 @@ long_description_content_type = text/markdown [flake8] doctests = True ignore = F403 -select = F,E722,W +select = F,E722 exclude = sympy/parsing/latex/_antlr/*, sympy/parsing/autolev/_antlr/*, sympy/parsing/autolev/test-examples/*, From f2537e648efa90af3c12b1ed6f1b0ee543a288dd Mon Sep 17 00:00:00 2001 From: Wes Galbraith Date: Wed, 17 Jun 2020 12:15:46 -0600 Subject: [PATCH 9/9] Remove extra checks for semidefinitness and consolidate methods --- sympy/matrices/eigen.py | 79 ++++++------------------------ sympy/matrices/tests/test_eigen.py | 38 -------------- 2 files changed, 15 insertions(+), 102 deletions(-) diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 78c1dbf61cd8..e9af03276c56 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -1,11 +1,9 @@ from types import FunctionType from collections import Counter -import itertools from mpmath import mp, workprec from mpmath.libmp.libmpf import prec_to_dps -from sympy import Symbol from sympy.core.compatibility import default_sort_key from sympy.core.evalf import DEFAULT_MAXPREC, PrecisionExhausted from sympy.core.logic import fuzzy_and, fuzzy_or @@ -701,18 +699,21 @@ def _is_positive_semidefinite(M): if nonnegative_diagonals and M.is_weakly_diagonally_dominant: return True - try: - return _is_positive_semidefinite_cholesky(M) - except Exception: - pass - - if M.rows < 5 or all(a.func != Symbol for a in M.atoms()): - try: - return _is_positive_semidefinite_evals(M) - except Exception: - pass - - return _is_positive_semidefinite_minors(M) + # uses Cholesky factorization with complete pivoting + # see http://eprints.ma.man.ac.uk/1199/1/covered/MIMS_ep2008_116.pdf + M = M.copy() + for k in range(M.rows): + pivot = max(range(k, M.rows), key=lambda i: M[i, i]) + if pivot > k: + M.col_swap(k, pivot) + M.row_swap(k, pivot) + if M[k, k].is_negative: + return False + M[k, k] = sqrt(M[k, k]) + M[k, (k+1):] /= M[k, k] + for j in range(k+1, M.rows): + M[(k+1):(j+1), j] -= M[k, (k+1):(j+1)].T * M[k, j] + return M[-1, -1].is_nonnegative def _is_negative_definite(M): @@ -754,56 +755,6 @@ def _is_positive_definite_GE(M): return True -def _is_positive_semidefinite_minors(M): - """A method to evaluate all principal minors for testing - positive-semidefiniteness by using Sylvestre's crieterion.""" - return all( - M[idx, idx].det(method='berkowitz').is_nonnegative - for minor_size in range(1, M.rows+1) - for idx in itertools.combinations(range(M.rows), minor_size) - ) - - -def _is_positive_semidefinite_evals(M): - """Determines if a matrix is positive semidefinite by checking - if all the eigenvalues are nonnegative""" - return all(eigenvalue.is_nonnegative for eigenvalue in M.eigenvals()) - - -def _is_positive_semidefinite_cholesky(M): - """Determines if a matrix is positive semidefinite by computing its - Cholesky decomposition with complete pivoting - - Any symmetric positive semidefinite matrix has a factorization - P.T * A * P = R.T * R, where - - R = [R_11 R_22] - [0 0 ], - - and R_11 is rxr upper triangular with positive diagonal elements, - rank(A) = r, and P is a permutation matrix. - - An algorithm for computing R is described in [1]. If at any point in - the algorithm a negative diagonal term is obtained, then A is not - positive semidefinite, so the algorithm exits early and returns False. - - [1] http://eprints.ma.man.ac.uk/1199/1/covered/MIMS_ep2008_116.pdf - """ - M = M.copy() - for k in range(M.rows): - pivot = max(range(k, M.rows), key=lambda i: M[i, i]) - if pivot > k: - M.col_swap(k, pivot) - M.row_swap(k, pivot) - if M[k, k].is_negative: - return False - M[k, k] = sqrt(M[k, k]) - M[k, (k+1):] /= M[k, k] - for j in range(k+1, M.rows): - M[(k+1):(j+1), j] -= M[k, (k+1):(j+1)].T * M[k, j] - return M[-1, -1].is_nonnegative - - _doc_positive_definite = \ r"""Finds out the definiteness of a matrix. diff --git a/sympy/matrices/tests/test_eigen.py b/sympy/matrices/tests/test_eigen.py index 412b69422bcb..b8a9cc6b309c 100644 --- a/sympy/matrices/tests/test_eigen.py +++ b/sympy/matrices/tests/test_eigen.py @@ -4,11 +4,6 @@ from sympy.matrices import eye, Matrix from sympy.core.singleton import S from sympy.testing.pytest import raises, XFAIL -from sympy.matrices.eigen import ( - _is_positive_semidefinite_cholesky, - _is_positive_semidefinite_evals, - _is_positive_semidefinite_minors -) from sympy.matrices.matrices import NonSquareMatrixError, MatrixError from sympy.simplify.simplify import simplify from sympy.matrices.immutable import ImmutableMatrix @@ -498,9 +493,6 @@ def test_definite(): m = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m) - assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -508,9 +500,6 @@ def test_definite(): m = Matrix([[5, 4], [4, 5]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m) - assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -519,9 +508,6 @@ def test_definite(): m = Matrix([[2, -1, -1], [-1, 2, -1], [-1, -1, 2]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m) - assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -529,9 +515,6 @@ def test_definite(): m = Matrix([[1, 2], [2, 4]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m.T + m) - assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -541,9 +524,6 @@ def test_definite(): m = Matrix([[2, 3], [4, 8]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m.T + m) - assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -552,9 +532,6 @@ def test_definite(): m = Matrix([[1, 2*I], [-I, 4]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m.H + m) - assert _is_positive_semidefinite_evals(m.H + m) - assert _is_positive_semidefinite_minors(m.H + m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -565,9 +542,6 @@ def test_definite(): m = Matrix([[a, 0, 0], [0, a, 0], [0, 0, a]]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m) - assert _is_positive_semidefinite_evals(m) - assert _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == False @@ -575,9 +549,6 @@ def test_definite(): m = Matrix([[b, 0, 0], [0, b, 0], [0, 0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False - assert not _is_positive_semidefinite_cholesky(m) - assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == True assert m.is_negative_semidefinite == True assert m.is_indefinite == False @@ -585,9 +556,6 @@ def test_definite(): m = Matrix([[a, 0], [0, b]]) assert m.is_positive_definite == False assert m.is_positive_semidefinite == False - assert not _is_positive_semidefinite_cholesky(m) - assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m) assert m.is_negative_definite == False assert m.is_negative_semidefinite == False assert m.is_indefinite == True @@ -604,9 +572,6 @@ def test_definite(): ]) assert m.is_positive_definite == True assert m.is_positive_semidefinite == True - assert _is_positive_semidefinite_cholesky(m.T + m) - assert _is_positive_semidefinite_evals(m.T + m) - assert _is_positive_semidefinite_minors(m.T + m) assert m.is_indefinite == False # test for issue 19547: https://github.com/sympy/sympy/issues/19547 @@ -617,6 +582,3 @@ def test_definite(): ]) assert not m.is_positive_definite assert not m.is_positive_semidefinite - assert not _is_positive_semidefinite_cholesky(m) - assert not _is_positive_semidefinite_evals(m) - assert not _is_positive_semidefinite_minors(m)