mjo/cone/symmetric_pd: new module for symmetric positive-definite matrices.

diff --git a/mjo/cone/all.py b/mjo/cone/all.py
index d6a6dd6a802939e41a77ee4ab2f4449b0482d50f..7ce56c0bf41f6f18123155c6b232bd52594b5905 100644 (file)
--- a/mjo/cone/all.py
+++ b/mjo/cone/all.py
@@ -10,5 +10,6 @@ from mjo.cone.nonnegative_orthant import *
 from mjo.cone.permutation_invariant import *
 from mjo.cone.rearrangement import *
 from mjo.cone.schur import *
+from mjo.cone.symmetric_pd import *
 from mjo.cone.symmetric_psd import *
 from mjo.cone.trivial_cone import *

diff --git a/mjo/cone/symmetric_pd.py b/mjo/cone/symmetric_pd.py
new file mode 100644 (file)
index 0000000..b6c8620
--- /dev/null
+++ b/mjo/cone/symmetric_pd.py
@@ -0,0 +1,144 @@
+"""
+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())