From e4e4685ec815a523d9ee73a61a3de8b8b2799a1b Mon Sep 17 00:00:00 2001 From: Michael Orlitzky 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