Trac #33100: drop special RDF/CDF cholesky() and is_positive_definite().

These methods are buggy, and the superclass implementation should work, so the easy path to not returning incorrect results is to just get rid of the buggy overrides for now.
sagemath · Jan 3, 2022 · f6e4b64 · f6e4b64
1 parent 38c86e5
commit f6e4b64
Show file tree

Hide file tree

Showing 2 changed files with 112 additions and 311 deletions.
diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx
@@ -12583,6 +12583,89 @@ cdef class Matrix(Matrix1):
             ...
             ValueError: matrix is not positive definite
 
+        ::
+
+            sage: A = matrix(RDF, [[ 3,  -6,   9,   6,  -9],
+            ....:                  [-6,  11, -16, -11,  17],
+            ....:                  [ 9, -16,  28,  16, -40],
+            ....:                  [ 6, -11,  16,   9, -19],
+            ....:                  [-9,  17, -40, -19,  68]])
+            sage: A.is_symmetric()
+            True
+            sage: A.eigenvalues()
+            [108.07..., 13.02..., -0.02..., -0.70..., -1.37...]
+            sage: A.is_positive_definite()
+            False
+            sage: A.cholesky()
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix is not positive definite
+
+        ::
+
+            sage: B = matrix(CDF, [[      2, 4 - 2*I, 2 + 2*I],
+            ....:                  [4 + 2*I,       8,    10*I],
+            ....:                  [2 - 2*I,   -10*I,      -3]])
+            sage: B.is_hermitian()
+            True
+            sage: [ev.real() for ev in B.eigenvalues()]
+            [15.88..., 0.08..., -8.97...]
+            sage: B.is_positive_definite()
+            False
+            sage: B.cholesky()
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix is not positive definite
+
+        A real matrix that is symmetric, Hermitian, and positive definite::
+
+            sage: M = matrix(RDF,[[ 1,  1,    1,     1,     1],
+            ....:                 [ 1,  5,   31,   121,   341],
+            ....:                 [ 1, 31,  341,  1555,  4681],
+            ....:                 [ 1,121, 1555,  7381, 22621],
+            ....:                 [ 1,341, 4681, 22621, 69905]])
+            sage: M.is_symmetric()
+            True
+            sage: M.is_hermitian()
+            True
+            sage: M.is_positive_definite()
+            True
+            sage: L = M.cholesky()
+            sage: L.round(6).zero_at(10^-10)
+            [   1.0    0.0         0.0        0.0     0.0]
+            [   1.0    2.0         0.0        0.0     0.0]
+            [   1.0   15.0   10.723805        0.0     0.0]
+            [   1.0   60.0   60.985814   7.792973     0.0]
+            [   1.0  170.0  198.623524  39.366567  1.7231]
+            sage: (L*L.transpose()).round(6).zero_at(10^-10)
+            [ 1.0     1.0     1.0     1.0     1.0]
+            [ 1.0     5.0    31.0   121.0   341.0]
+            [ 1.0    31.0   341.0  1555.0  4681.0]
+            [ 1.0   121.0  1555.0  7381.0 22621.0]
+            [ 1.0   341.0  4681.0 22621.0 69905.0]
+
+        A complex matrix that is Hermitian and positive definite::
+
+            sage: A = matrix(CDF, [[        23,  17*I + 3,  24*I + 25,     21*I],
+            ....:                  [ -17*I + 3,        38, -69*I + 89, 7*I + 15],
+            ....:                  [-24*I + 25, 69*I + 89,        976, 24*I + 6],
+            ....:                  [     -21*I, -7*I + 15,  -24*I + 6,       28]])
+            sage: A.is_hermitian()
+            True
+            sage: A.is_positive_definite()
+            True
+            sage: L = A.cholesky()
+            sage: L.round(6).zero_at(10^-10)
+            [               4.795832                     0.0                    0.0       0.0]
+            [  0.625543 - 3.544745*I                5.004346                    0.0       0.0]
+            [   5.21286 - 5.004346*I 13.588189 + 10.721116*I              24.984023       0.0]
+            [            -4.378803*I  -0.104257 - 0.851434*I  -0.21486 + 0.371348*I  2.811799]
+            sage: (L*L.conjugate_transpose()).round(6).zero_at(10^-10)
+            [         23.0  3.0 + 17.0*I 25.0 + 24.0*I        21.0*I]
+            [ 3.0 - 17.0*I          38.0 89.0 - 69.0*I  15.0 + 7.0*I]
+            [25.0 - 24.0*I 89.0 + 69.0*I         976.0  6.0 + 24.0*I]
+            [      -21.0*I  15.0 - 7.0*I  6.0 - 24.0*I          28.0]
+
         TESTS:
 
         This verifies that :trac:`11274` is resolved::
@@ -14641,6 +14724,14 @@ cdef class Matrix(Matrix1):
             sage: [A[:i,:i].determinant() for i in range(1,A.nrows()+1)]
             [2, 2, 6]
 
+        A large random matrix that is guaranteed by theory to be
+        positive definite::
+
+            sage: R = matrix.random(CDF, 200)
+            sage: H = R.conjugate_transpose()*R
+            sage: H.is_positive_definite()
+            True
+
         TESTS:
 
         If the base ring does not make sense as a subfield of the real
@@ -14677,6 +14768,27 @@ cdef class Matrix(Matrix1):
             True
             sage: matrix.identity(SR,4).is_positive_definite()
             True
+
+        A trivially small case::
+
+            sage: S = matrix(CDF, [])
+            sage: S.nrows(), S.ncols()
+            (0, 0)
+            sage: S.is_positive_definite()
+            True
+
+        A rectangular matrix will never be positive definite::
+
+            sage: R = matrix(RDF, 2, 3, range(6))
+            sage: R.is_positive_definite()
+            False
+
+        A non-Hermitian matrix will never be positive definite::
+
+            sage: T = matrix(CDF, 8, 8, range(64))
+            sage: T.is_positive_definite()
+            False
+
         """
         result = self._is_positive_definite_or_semidefinite(False)
         if certificate: