-"""
+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
from sage.all import *
from mjo.cone.symmetric_psd import random_symmetric_psd
-def is_symmetric_pd(A):
- """
- Determine whether or not the matrix ``A`` is symmetric
- positive-definite.
-
- INPUT:
-
- - ``A`` - The matrix in question.
-
- OUTPUT:
-
- Either ``True`` if ``A`` is symmetric positive-definite, or
- ``False`` otherwise.
-
- SETUP::
-
- sage: from mjo.cone.symmetric_pd import (is_symmetric_pd,
- ....: random_symmetric_pd)
-
- EXAMPLES:
-
- The identity matrix is obviously symmetric and positive-definite::
-
- sage: A = identity_matrix(ZZ, ZZ.random_element(10))
- sage: is_symmetric_pd(A)
- True
-
- The following matrix is symmetric but not positive definite::
-
- sage: A = matrix(ZZ, [[1, 2], [2, 1]])
- sage: is_symmetric_pd(A)
- False
-
- This matrix isn't even symmetric::
-
- sage: A = matrix(ZZ, [[1, 2], [3, 4]])
- sage: is_symmetric_pd(A)
- False
-
- The trivial matrix in a trivial space is trivially symmetric and
- positive-definite::
-
- sage: A = matrix(QQ, 0,0)
- sage: is_symmetric_pd(A)
- True
-
- The implementation of "is_positive_definite" in Sage leaves a bit to
- be desired. This matrix is technically positive definite over the
- real numbers, but isn't symmetric::
-
- sage: A = matrix(RR,[[1,1],[-1,1]])
- sage: A.is_positive_definite()
- Traceback (most recent call last):
- ...
- ValueError: Could not see Real Field with 53 bits of precision as a
- subring of the real or complex numbers
- sage: is_symmetric_pd(A)
- False
-
- TESTS:
-
- Every gram matrix is positive-definite, and we can sum two
- positive-definite matrices (``A`` and its transpose) to get a new
- positive-definite matrix that happens to be symmetric::
-
- sage: set_random_seed()
- sage: V = VectorSpace(QQ, ZZ.random_element(5))
- sage: A = random_symmetric_pd(V)
- sage: is_symmetric_pd(A)
- True
-
- """
-
- if A.base_ring() == SR:
- msg = 'The matrix ``A`` cannot be symbolic.'
- raise ValueError.new(msg)
-
- # First make sure that ``A`` is symmetric.
- if not A.is_symmetric():
- return False
-
- # If ``A`` is symmetric, we only need to check that it is positive
- # definite. For that we can consult its minimum eigenvalue, which
- # should be greater than zero. Since ``A`` is symmetric, its
- # eigenvalues are guaranteed to be real.
- if A.is_zero():
- # A is trivial... so trivially positive-definite.
- return True
- else:
- return min(A.eigenvalues()) > 0
-
-
def random_symmetric_pd(V):
r"""
Generate a random symmetric positive-definite matrix over the
# 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
projects / sage.d.git / blobdiff