--- /dev/null
+"""
+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
+SEMI-definite cone.
+"""
+
+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
+ vector space ``V``. That is, the returned matrix will be a linear
+ transformation on ``V``, with the same base ring as ``V``.
+
+ (We take a very loose interpretation of "random," here.)
+
+ INPUT:
+
+ - ``V`` - The vector space on which the returned matrix will act.
+
+ OUTPUT:
+
+ A random symmetric positive-definite matrix, i.e. a linear
+ transformation from ``V`` to itself.
+
+ ALGORITHM:
+
+ We request a full-rank positive-semidefinite matrix from the
+ :func:`mjo.cone.symmetric_psd.random_symmetric_psd` function.
+
+ SETUP::
+
+ sage: from mjo.cone.symmetric_pd import random_symmetric_pd
+
+ EXAMPLES:
+
+ Well, it doesn't crash at least::
+
+ sage: set_random_seed()
+ sage: V = VectorSpace(QQ, 2)
+ sage: A = random_symmetric_pd(V)
+ sage: A.matrix_space()
+ Full MatrixSpace of 2 by 2 dense matrices over Rational Field
+
+ """
+ # We accept_zero because the trivial matrix is (trivially)
+ # symmetric and positive-definite, but it's also "zero." The rank
+ # condition ensures that we don't get the zero matrix in a
+ # nontrivial space.
+ return random_symmetric_psd(V, rank=V.dimension())