Repository: sympy/sympy

Issue:
Performance regression in Cholesky solve
Offending PR: https://github.com/sympy/sympy/pull/14474

Benchmark:

http://www.moorepants.info/misc/sympy-asv/#solve.TimeMatrixSolve2.time_cholesky_solve

Discussion:

Patch:
diff --git a/data/TeXmacs/bin/tm_sympy b/data/TeXmacs/bin/tm_sympy
index 2233dc6ef6..5e1a0de073 100755
--- a/data/TeXmacs/bin/tm_sympy
+++ b/data/TeXmacs/bin/tm_sympy
@@ -65,7 +65,7 @@ f = Function('f')
 Gamma, Zeta = gamma, zeta
 
 _greek = 'alpha beta gamma delta epsilon zeta eta '  \
-         'theta iota kappa mu nu xi omicron rho ' \
+         'theta iota kappa mu nu xi omicron pi rho ' \
          'sigma tau upsilon phi chi psi omega'
 
 for _symbol in _greek.split(' '):
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index cee6f0e7ef..f2ac15f25e 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -129,22 +129,32 @@ def __getitem__(self, key):
     def __setitem__(self, key, value):
         raise NotImplementedError()
 
-    def _cholesky(self):
+    def _cholesky(self, hermitian=True):
         """Helper function of cholesky.
         Without the error checks.
         To be used privately.
         Implements the Cholesky-Banachiewicz algorithm.
+        Returns L such that L*L.H == self if hermitian flag is True,
+        or L*L.T == self if hermitian is False.
         """
         L = zeros(self.rows, self.rows)
-        for i in range(self.rows):
-            for j in range(i):
-                L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
-                    sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
-            Lii2 = expand_mul(self[i, i] -
-                sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
-            if Lii2.is_positive is False:
-                raise ValueError("Matrix must be positive-definite")
-            L[i, i] = sqrt(Lii2)
+        if hermitian:
+            for i in range(self.rows):
+                for j in range(i):
+                    L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
+                        sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
+                Lii2 = expand_mul(self[i, i] -
+                    sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
+                if Lii2.is_positive is False:
+                    raise ValueError("Matrix must be positive-definite")
+                L[i, i] = sqrt(Lii2)
+        else:
+            for i in range(self.rows):
+                for j in range(i):
+                    L[i, j] = (1 / L[j, j])*(self[i, j] -
+                        sum(L[i, k]*L[j, k] for k in range(j)))
+                L[i, i] = sqrt(self[i, i] -
+                    sum(L[i, k]**2 for k in range(i)))
         return self._new(L)
 
     def _diagonal_solve(self, rhs):
@@ -281,22 +291,31 @@ def _eval_tolist(self):
         cols = self.cols
         return [mat[i*cols:(i + 1)*cols] for i in range(self.rows)]
 
-    def _LDLdecomposition(self):
+    def _LDLdecomposition(self, hermitian=True):
         """Helper function of LDLdecomposition.
         Without the error checks.
         To be used privately.
+        Returns L and D such that L*D*L.H == self if hermitian flag is True,
+        or L*D*L.T == self if hermitian is False.
         """
         # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2
         D = zeros(self.rows, self.rows)
         L = eye(self.rows)
-        for i in range(self.rows):
-            for j in range(i):
-                L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
-                    L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
-            D[i, i] = expand_mul(self[i, i] -
-                sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
-            if D[i, i].is_positive is False:
-                raise ValueError("Matrix must be positive-definite")
+        if hermitian:
+            for i in range(self.rows):
+                for j in range(i):
+                    L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
+                        L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
+                D[i, i] = expand_mul(self[i, i] -
+                    sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
+                if D[i, i].is_positive is False:
+                    raise ValueError("Matrix must be positive-definite")
+        else:
+            for i in range(self.rows):
+                for j in range(i):
+                    L[i, j] = (1 / D[j, j])*(self[i, j] - sum(
+                        L[i, k]*L[j, k]*D[k, k] for k in range(j)))
+                D[i, i] = self[i, i] - sum(L[i, k]**2*D[k, k] for k in range(i))
         return self._new(L), self._new(D)
 
     def _lower_triangular_solve(self, rhs):
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index f5eab5615e..952d7cf062 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -2361,22 +2361,31 @@ def cholesky_solve(self, rhs):
         QRsolve
         pinv_solve
         """
-        if self.is_hermitian:
-            L = self._cholesky()
+        hermitian = True
+        if self.is_symmetric():
+            hermitian = False
+            L = self._cholesky(hermitian=hermitian)
+        elif self.is_hermitian:
+            L = self._cholesky(hermitian=hermitian)
         elif self.rows >= self.cols:
-            L = (self.H * self)._cholesky()
+            L = (self.H * self)._cholesky(hermitian=hermitian)
             rhs = self.H * rhs
         else:
             raise NotImplementedError('Under-determined System. '
                                       'Try M.gauss_jordan_solve(rhs)')
         Y = L._lower_triangular_solve(rhs)
-        return (L.H)._upper_triangular_solve(Y)
+        if hermitian:
+            return (L.H)._upper_triangular_solve(Y)
+        else:
+            return (L.T)._upper_triangular_solve(Y)
 
-    def cholesky(self):
-        """Returns the Cholesky decomposition L of a matrix A
-        such that L * L.H = A
+    def cholesky(self, hermitian=True):
+        """Returns the Cholesky-type decomposition L of a matrix A
+        such that L * L.H == A if hermitian flag is True,
+        or L * L.T == A if hermitian is False.
 
-        A must be a Hermitian positive-definite matrix.
+        A must be a Hermitian positive-definite matrix if hermitian is True,
+        or a symmetric matrix if it is False.
 
         Examples
         ========
@@ -2407,6 +2416,18 @@ def cholesky(self):
         [   9, 3*I],
         [-3*I,   5]])
 
+        Non-hermitian Cholesky-type decomposition may be useful when the
+        matrix is not positive-definite.
+
+        >>> A = Matrix([[1, 2], [2, 1]])
+        >>> L = A.cholesky(hermitian=False)
+        >>> L
+        Matrix([
+        [1,         0],
+        [2, sqrt(3)*I]])
+        >>> L*L.T == A
+        True
+
         See Also
         ========
 
@@ -2417,9 +2438,11 @@ def cholesky(self):
 
         if not self.is_square:
             raise NonSquareMatrixError("Matrix must be square.")
-        if not self.is_hermitian:
+        if hermitian and not self.is_hermitian:
             raise ValueError("Matrix must be Hermitian.")
-        return self._cholesky()
+        if not hermitian and not self.is_symmetric():
+            raise ValueError("Matrix must be symmetric.")
+        return self._cholesky(hermitian=hermitian)
 
     def condition_number(self):
         """Returns the condition number of a matrix.
@@ -3095,12 +3118,14 @@ def key2ij(self, key):
         else:
             return divmod(a2idx_(key, len(self)), self.cols)
 
-    def LDLdecomposition(self):
+    def LDLdecomposition(self, hermitian=True):
         """Returns the LDL Decomposition (L, D) of matrix A,
-        such that L * D * L.H == A
+        such that L * D * L.H == A if hermitian flag is True, or
+        L * D * L.T == A if hermitian is False.
         This method eliminates the use of square root.
         Further this ensures that all the diagonal entries of L are 1.
-        A must be a Hermitian positive-definite matrix.
+        A must be a Hermitian positive-definite matrix if hermitian is True,
+        or a symmetric matrix otherwise.
 
         Examples
         ========
@@ -3146,9 +3171,11 @@ def LDLdecomposition(self):
         """
         if not self.is_square:
             raise NonSquareMatrixError("Matrix must be square.")
-        if not self.is_hermitian:
+        if hermitian and not self.is_hermitian:
             raise ValueError("Matrix must be Hermitian.")
-        return self._LDLdecomposition()
+        if not hermitian and not self.is_symmetric():
+            raise ValueError("Matrix must be symmetric.")
+        return self._LDLdecomposition(hermitian=hermitian)
 
     def LDLsolve(self, rhs):
         """Solves Ax = B using LDL decomposition,
@@ -3179,17 +3206,24 @@ def LDLsolve(self, rhs):
         QRsolve
         pinv_solve
         """
-        if self.is_hermitian:
-            L, D = self.LDLdecomposition()
+        hermitian = True
+        if self.is_symmetric():
+            hermitian = False
+            L, D = self.LDLdecomposition(hermitian=hermitian)
+        elif self.is_hermitian:
+            L, D = self.LDLdecomposition(hermitian=hermitian)
         elif self.rows >= self.cols:
-            L, D = (self.H * self).LDLdecomposition()
+            L, D = (self.H * self).LDLdecomposition(hermitian=hermitian)
             rhs = self.H * rhs
         else:
             raise NotImplementedError('Under-determined System. '
                                       'Try M.gauss_jordan_solve(rhs)')
         Y = L._lower_triangular_solve(rhs)
         Z = D._diagonal_solve(Y)
-        return (L.H)._upper_triangular_solve(Z)
+        if hermitian:
+            return (L.H)._upper_triangular_solve(Z)
+        else:
+            return (L.T)._upper_triangular_solve(Z)
 
     def lower_triangular_solve(self, rhs):
         """Solves Ax = B, where A is a lower triangular matrix.


Test patch:
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index d82aeeb017..0a07df4541 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -2027,10 +2027,18 @@ def test_cholesky():
     raises(ValueError, lambda: Matrix(((1, 2), (3, 4))).cholesky())
     raises(ValueError, lambda: Matrix(((5 + I, 0), (0, 1))).cholesky())
     raises(ValueError, lambda: Matrix(((1, 5), (5, 1))).cholesky())
+    raises(ValueError, lambda: Matrix(((1, 2), (3, 4))).cholesky(hermitian=False))
+    assert Matrix(((5 + I, 0), (0, 1))).cholesky(hermitian=False) == Matrix([
+        [sqrt(5 + I), 0], [0, 1]])
+    A = Matrix(((1, 5), (5, 1)))
+    L = A.cholesky(hermitian=False)
+    assert L == Matrix([[1, 0], [5, 2*sqrt(6)*I]])
+    assert L*L.T == A
     A = Matrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))
-    assert A.cholesky() * A.cholesky().T == A
-    assert A.cholesky().is_lower
-    assert A.cholesky() == Matrix([[5, 0, 0], [3, 3, 0], [-1, 1, 3]])
+    L = A.cholesky()
+    assert L * L.T == A
+    assert L.is_lower
+    assert L == Matrix([[5, 0, 0], [3, 3, 0], [-1, 1, 3]])
     A = Matrix(((4, -2*I, 2 + 2*I), (2*I, 2, -1 + I), (2 - 2*I, -1 - I, 11)))
     assert A.cholesky() == Matrix(((2, 0, 0), (I, 1, 0), (1 - I, 0, 3)))
 
@@ -2040,6 +2048,10 @@ def test_LDLdecomposition():
     raises(ValueError, lambda: Matrix(((1, 2), (3, 4))).LDLdecomposition())
     raises(ValueError, lambda: Matrix(((5 + I, 0), (0, 1))).LDLdecomposition())
     raises(ValueError, lambda: Matrix(((1, 5), (5, 1))).LDLdecomposition())
+    raises(ValueError, lambda: Matrix(((1, 2), (3, 4))).LDLdecomposition(hermitian=False))
+    A = Matrix(((1, 5), (5, 1)))
+    L, D = A.LDLdecomposition(hermitian=False)
+    assert L * D * L.T == A
     A = Matrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))
     L, D = A.LDLdecomposition()
     assert L * D * L.T == A
@@ -2069,6 +2081,11 @@ def test_cholesky_solve():
     b = A*x
     soln = A.cholesky_solve(b)
     assert soln == x
+    A = Matrix(((1, 5), (5, 1)))
+    x = Matrix((4, -3))
+    b = A*x
+    soln = A.cholesky_solve(b)
+    assert soln == x
     A = Matrix(((9, 3*I), (-3*I, 5)))
     x = Matrix((-2, 1))
     b = A*x
@@ -2079,6 +2096,11 @@ def test_cholesky_solve():
     b = A*x
     soln = A.cholesky_solve(b)
     assert expand_mul(soln) == x
+    a00, a01, a11, b0, b1 = symbols('a00, a01, a11, b0, b1')
+    A = Matrix(((a00, a01), (a01, a11)))
+    b = Matrix((b0, b1))
+    x = A.cholesky_solve(b)
+    assert simplify(A*x) == b
 
 
 def test_LDLsolve():