From e4e4685ec815a523d9ee73a61a3de8b8b2799a1b Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 30 Nov 2018 15:17:16 -0500
Subject: [PATCH] mjo/cone/symmetric_pd: add inverse_symmetric_pd() function.

---
 mjo/cone/symmetric_pd.py | 59 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 58 insertions(+), 1 deletion(-)

diff --git a/mjo/cone/symmetric_pd.py b/mjo/cone/symmetric_pd.py
index b6c8620..9c0d6ba 100644
--- a/mjo/cone/symmetric_pd.py
+++ b/mjo/cone/symmetric_pd.py
@@ -110,7 +110,7 @@ def random_symmetric_pd(V):
 
     INPUT:
 
-    - ``V`` - The vector space on which the returned matrix will act.
+    - ``V`` -- The vector space on which the returned matrix will act.
 
     OUTPUT:
 
@@ -142,3 +142,60 @@ def random_symmetric_pd(V):
     # condition ensures that we don't get the zero matrix in a
     # nontrivial space.
     return random_symmetric_psd(V, rank=V.dimension())
+
+
+def inverse_symmetric_pd(M):
+    r"""
+    Compute the inverse of a symmetric positive-definite matrix.
+
+    The symmetric positive-definite matrix ``M`` has a Cholesky
+    factorization, which can be used to invert ``M`` in a fast,
+    numerically-stable way.
+
+    INPUT:
+
+    - ``M`` -- The matrix to invert.
+
+    OUTPUT:
+
+    The inverse of ``M``. The Cholesky factorization uses roots, so the
+    returned matrix will have as its base ring the algebraic closure of
+    the base ring of ``M``.
+
+    SETUP::
+
+        sage: from mjo.cone.symmetric_pd import inverse_symmetric_pd
+
+    EXAMPLES:
+
+    If we start with a rational matrix, its inverse will be over a field
+    extension of that rationals that is necessarily contained in
+    ``QQbar`` (which contails ALL roots)::
+
+        sage: M = matrix(QQ, [[3,1],[1,5]])
+        sage: M_inv = inverse_symmetric_pd(M)
+        sage: M_inv.base_ring().is_subring(QQbar)
+        True
+
+    TESTS:
+
+    This "inverse" should actually act like an inverse::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: V = VectorSpace(QQ,n)
+        sage: from mjo.cone.symmetric_pd import random_symmetric_pd
+        sage: M = random_symmetric_pd(V)
+        sage: M_inv = inverse_symmetric_pd(M)
+        sage: I = identity_matrix(QQ,n)
+        sage: M_inv * M == I
+        True
+        sage: M * M_inv == I
+        True
+
+    """
+    # M = L * L.transpose()
+    # The default "echelonize" inverse() method works just fine for
+    # triangular matrices.
+    L_inverse = M.cholesky().inverse()
+    return L_inverse.transpose()*L_inverse
-- 
2.44.2