# 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