Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix 19547 is positive semidefinite incorrect #19573

Merged
merged 11 commits into from Jun 18, 2020
Merged

Fix 19547 is positive semidefinite incorrect #19573

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
bf4fcef
Add check on eigenvalues to determine semidefiniteness
galbwe Jun 16, 2020
1bf89e8
Got cholesky algorithm to pass existing tests
galbwe Jun 17, 2020
74e6807
Merge branch 'master' of https://github.com/sympy/sympy into fix-1954…
galbwe Jun 17, 2020
50fbb84
vectorize inner loops
galbwe Jun 17, 2020
b07e43f
Add reference
galbwe Jun 17, 2020
9e80ddf
Add individual tests for subroutines
galbwe Jun 17, 2020
1682ac3
Use shorter function names
galbwe Jun 17, 2020
7612cb8
Merge branch 'master' of https://github.com/sympy/sympy into fix-1954…
galbwe Jun 17, 2020
7bb906e
Remove trailing whitespace and update flake8 config
galbwe Jun 17, 2020
b15fcf0
Revert setup.cfg
galbwe Jun 17, 2020
f2537e6
Remove extra checks for semidefinitness and consolidate methods
galbwe Jun 17, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
52 changes: 52 additions & 0 deletions sympy/matrices/eigen.py
View file
Open in desktop
Expand Up @@ -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
Expand Down Expand Up @@ -700,6 +701,17 @@ 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)
except Exception:
pass

return _is_positive_semidefinite_by_minors(M)


Expand Down Expand Up @@ -752,6 +764,46 @@ 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())


def _is_positive_semidefinite_by_cholesky_factorization_with_pivots(M):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are too long names that can easily exceed 80 characters.
What about making them shorter like _is_positive_semidefinite_cholesky.
I don't think that it should need much mathematical meaning, especially if these functions are only used internaly.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree they are a bit verbose. I tend to err on the side of more specificity in my variable names and comments. I shortened them.

"""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.

Expand Down
40 changes: 40 additions & 0 deletions sympy/matrices/tests/test_eigen.py
View file
Open in desktop
Expand Up @@ -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
Expand Down Expand Up @@ -493,13 +498,19 @@ 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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a trailing whitespace at the end of line 503 which is causing quality tests to fail, please fix that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My bad there. I removed whitespace and updated the flake8 part of setup.cfg to catch it earlier next time.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess adding whitespace to setup.cfg was a bad idea, since it caused a whole bunch of linting errors in the pipeline. Removed in latest commit.

assert m.is_negative_definite == False
assert m.is_negative_semidefinite == False
assert m.is_indefinite == False

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
Expand All @@ -508,13 +519,19 @@ 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

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
Expand All @@ -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
Expand All @@ -541,20 +565,29 @@ 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

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

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
Expand All @@ -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
Expand All @@ -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)