INPUT:
- - ``V`` - The vector space on which the returned matrix will act.
+ - ``V`` -- The vector space on which the returned matrix will act.
OUTPUT:
# 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