X-Git-Url: http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=blobdiff_plain;f=mjo%2Fcone%2Fsymmetric_pd.py;h=5dcacfb8ea82b658a0c6d9bb46aed9a3034e2512;hp=9c0d6bae7573a5c085a96c797b9b991b29ec0e88;hb=c0a2a45d526131c9218f6b63b1728308ac963535;hpb=e4e4685ec815a523d9ee73a61a3de8b8b2799a1b diff --git a/mjo/cone/symmetric_pd.py b/mjo/cone/symmetric_pd.py index 9c0d6ba..5dcacfb 100644 --- a/mjo/cone/symmetric_pd.py +++ b/mjo/cone/symmetric_pd.py @@ -1,4 +1,4 @@ -""" +r""" The symmetric positive definite cone `$S^{n}_{++}$` is the cone consisting of all symmetric positive-definite matrices (as a subset of $\mathbb{R}^{n \times n}$`. It is the interior of the symmetric positive @@ -142,60 +142,3 @@ 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