mjo/cone/symmetric_pd.py

   1 r"""
   2 The symmetric positive definite cone `$S^{n}_{++}$` is the cone
   3 consisting of all symmetric positive-definite matrices (as a subset of
   4 $\mathbb{R}^{n \times n}$`. It is the interior of the symmetric positive
   5 SEMI-definite cone.
   6 """
   7
   8 from sage.all import *
   9 from mjo.cone.symmetric_psd import random_symmetric_psd
  10
  11 def is_symmetric_pd(A):
  12     """
  13     Determine whether or not the matrix ``A`` is symmetric
  14     positive-definite.
  15
  16     INPUT:
  17
  18     - ``A`` - The matrix in question.
  19
  20     OUTPUT:
  21
  22     Either ``True`` if ``A`` is symmetric positive-definite, or
  23     ``False`` otherwise.
  24
  25     SETUP::
  26
  27         sage: from mjo.cone.symmetric_pd import (is_symmetric_pd,
  28         ....:                                    random_symmetric_pd)
  29
  30     EXAMPLES:
  31
  32     The identity matrix is obviously symmetric and positive-definite::
  33
  34         sage: set_random_seed()
  35         sage: A = identity_matrix(ZZ, ZZ.random_element(10))
  36         sage: is_symmetric_pd(A)
  37         True
  38
  39     The following matrix is symmetric but not positive definite::
  40
  41         sage: A = matrix(ZZ, [[1, 2], [2, 1]])
  42         sage: is_symmetric_pd(A)
  43         False
  44
  45     This matrix isn't even symmetric::
  46
  47         sage: A = matrix(ZZ, [[1, 2], [3, 4]])
  48         sage: is_symmetric_pd(A)
  49         False
  50
  51     The trivial matrix in a trivial space is trivially symmetric and
  52     positive-definite::
  53
  54         sage: A = matrix(QQ, 0,0)
  55         sage: is_symmetric_pd(A)
  56         True
  57
  58     The implementation of "is_positive_definite" in Sage leaves a bit to
  59     be desired. This matrix is technically positive definite over the
  60     real numbers, but isn't symmetric::
  61
  62         sage: A = matrix(RR,[[1,1],[-1,1]])
  63         sage: A.is_positive_definite()
  64         Traceback (most recent call last):
  65         ...
  66         ValueError: Could not see Real Field with 53 bits of precision as a
  67         subring of the real or complex numbers
  68         sage: is_symmetric_pd(A)
  69         False
  70
  71     TESTS:
  72
  73     Every gram matrix is positive-definite, and we can sum two
  74     positive-definite matrices (``A`` and its transpose) to get a new
  75     positive-definite matrix that happens to be symmetric::
  76
  77         sage: set_random_seed()
  78         sage: V = VectorSpace(QQ, ZZ.random_element(5))
  79         sage: A = random_symmetric_pd(V)
  80         sage: is_symmetric_pd(A)
  81         True
  82
  83     """
  84
  85     if A.base_ring() == SR:
  86         msg = 'The matrix ``A`` cannot be symbolic.'
  87         raise ValueError.new(msg)
  88
  89     # First make sure that ``A`` is symmetric.
  90     if not A.is_symmetric():
  91         return False
  92
  93     # If ``A`` is symmetric, we only need to check that it is positive
  94     # definite. For that we can consult its minimum eigenvalue, which
  95     # should be greater than zero. Since ``A`` is symmetric, its
  96     # eigenvalues are guaranteed to be real.
  97     if A.is_zero():
  98         # A is trivial... so trivially positive-definite.
  99         return True
 100     else:
 101         return min(A.eigenvalues()) > 0
 102
 103
 104 def random_symmetric_pd(V):
 105     r"""
 106     Generate a random symmetric positive-definite matrix over the
 107     vector space ``V``. That is, the returned matrix will be a linear
 108     transformation on ``V``, with the same base ring as ``V``.
 109
 110     (We take a very loose interpretation of "random," here.)
 111
 112     INPUT:
 113
 114     - ``V`` -- The vector space on which the returned matrix will act.
 115
 116     OUTPUT:
 117
 118     A random symmetric positive-definite matrix, i.e. a linear
 119     transformation from ``V`` to itself.
 120
 121     ALGORITHM:
 122
 123     We request a full-rank positive-semidefinite matrix from the
 124     :func:`mjo.cone.symmetric_psd.random_symmetric_psd` function.
 125
 126     SETUP::
 127
 128         sage: from mjo.cone.symmetric_pd import random_symmetric_pd
 129
 130     EXAMPLES:
 131
 132     Well, it doesn't crash at least::
 133
 134         sage: set_random_seed()
 135         sage: V = VectorSpace(QQ, 2)
 136         sage: A = random_symmetric_pd(V)
 137         sage: A.matrix_space()
 138         Full MatrixSpace of 2 by 2 dense matrices over Rational Field
 139
 140     """
 141     # We accept_zero because the trivial matrix is (trivially)
 142     # symmetric and positive-definite, but it's also "zero." The rank
 143     # condition ensures that we don't get the zero matrix in a
 144     # nontrivial space.
 145     return random_symmetric_psd(V, rank=V.dimension())