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: A = identity_matrix(ZZ, ZZ.random_element(10))
  35         sage: is_symmetric_pd(A)
  36         True
  37
  38     The following matrix is symmetric but not positive definite::
  39
  40         sage: A = matrix(ZZ, [[1, 2], [2, 1]])
  41         sage: is_symmetric_pd(A)
  42         False
  43
  44     This matrix isn't even symmetric::
  45
  46         sage: A = matrix(ZZ, [[1, 2], [3, 4]])
  47         sage: is_symmetric_pd(A)
  48         False
  49
  50     The trivial matrix in a trivial space is trivially symmetric and
  51     positive-definite::
  52
  53         sage: A = matrix(QQ, 0,0)
  54         sage: is_symmetric_pd(A)
  55         True
  56
  57     The implementation of "is_positive_definite" in Sage leaves a bit to
  58     be desired. This matrix is technically positive definite over the
  59     real numbers, but isn't symmetric::
  60
  61         sage: A = matrix(RR,[[1,1],[-1,1]])
  62         sage: A.is_positive_definite()
  63         Traceback (most recent call last):
  64         ...
  65         ValueError: Could not see Real Field with 53 bits of precision as a
  66         subring of the real or complex numbers
  67         sage: is_symmetric_pd(A)
  68         False
  69
  70     TESTS:
  71
  72     Every gram matrix is positive-definite, and we can sum two
  73     positive-definite matrices (``A`` and its transpose) to get a new
  74     positive-definite matrix that happens to be symmetric::
  75
  76         sage: set_random_seed()
  77         sage: V = VectorSpace(QQ, ZZ.random_element(5))
  78         sage: A = random_symmetric_pd(V)
  79         sage: is_symmetric_pd(A)
  80         True
  81
  82     """
  83
  84     if A.base_ring() == SR:
  85         msg = 'The matrix ``A`` cannot be symbolic.'
  86         raise ValueError.new(msg)
  87
  88     # First make sure that ``A`` is symmetric.
  89     if not A.is_symmetric():
  90         return False
  91
  92     # If ``A`` is symmetric, we only need to check that it is positive
  93     # definite. For that we can consult its minimum eigenvalue, which
  94     # should be greater than zero. Since ``A`` is symmetric, its
  95     # eigenvalues are guaranteed to be real.
  96     if A.is_zero():
  97         # A is trivial... so trivially positive-definite.
  98         return True
  99     else:
 100         return min(A.eigenvalues()) > 0
 101
 102
 103 def random_symmetric_pd(V):
 104     r"""
 105     Generate a random symmetric positive-definite matrix over the
 106     vector space ``V``. That is, the returned matrix will be a linear
 107     transformation on ``V``, with the same base ring as ``V``.
 108
 109     (We take a very loose interpretation of "random," here.)
 110
 111     INPUT:
 112
 113     - ``V`` -- The vector space on which the returned matrix will act.
 114
 115     OUTPUT:
 116
 117     A random symmetric positive-definite matrix, i.e. a linear
 118     transformation from ``V`` to itself.
 119
 120     ALGORITHM:
 121
 122     We request a full-rank positive-semidefinite matrix from the
 123     :func:`mjo.cone.symmetric_psd.random_symmetric_psd` function.
 124
 125     SETUP::
 126
 127         sage: from mjo.cone.symmetric_pd import random_symmetric_pd
 128
 129     EXAMPLES:
 130
 131     Well, it doesn't crash at least::
 132
 133         sage: set_random_seed()
 134         sage: V = VectorSpace(QQ, 2)
 135         sage: A = random_symmetric_pd(V)
 136         sage: A.matrix_space()
 137         Full MatrixSpace of 2 by 2 dense matrices over Rational Field
 138
 139     """
 140     # We accept_zero because the trivial matrix is (trivially)
 141     # symmetric and positive-definite, but it's also "zero." The rank
 142     # condition ensures that we don't get the zero matrix in a
 143     # nontrivial space.
 144     return random_symmetric_psd(V, rank=V.dimension())
 145
 146
 147 def inverse_symmetric_pd(M):
 148     r"""
 149     Compute the inverse of a symmetric positive-definite matrix.
 150
 151     The symmetric positive-definite matrix ``M`` has a Cholesky
 152     factorization, which can be used to invert ``M`` in a fast,
 153     numerically-stable way.
 154
 155     INPUT:
 156
 157     - ``M`` -- The matrix to invert.
 158
 159     OUTPUT:
 160
 161     The inverse of ``M``. The Cholesky factorization uses roots, so the
 162     returned matrix will have as its base ring the algebraic closure of
 163     the base ring of ``M``.
 164
 165     SETUP::
 166
 167         sage: from mjo.cone.symmetric_pd import inverse_symmetric_pd
 168
 169     EXAMPLES:
 170
 171     If we start with a rational matrix, its inverse will be over a field
 172     extension of that rationals that is necessarily contained in
 173     ``QQbar`` (which contails ALL roots)::
 174
 175         sage: M = matrix(QQ, [[3,1],[1,5]])
 176         sage: M_inv = inverse_symmetric_pd(M)
 177         sage: M_inv.base_ring().is_subring(QQbar)
 178         True
 179
 180     TESTS:
 181
 182     This "inverse" should actually act like an inverse::
 183
 184         sage: set_random_seed()
 185         sage: n = ZZ.random_element(5)
 186         sage: V = VectorSpace(QQ,n)
 187         sage: from mjo.cone.symmetric_pd import random_symmetric_pd
 188         sage: M = random_symmetric_pd(V)
 189         sage: M_inv = inverse_symmetric_pd(M)
 190         sage: I = identity_matrix(QQ,n)
 191         sage: M_inv * M == I
 192         True
 193         sage: M * M_inv == I
 194         True
 195
 196     """
 197     # M = L * L.transpose()
 198     # The default "echelonize" inverse() method works just fine for
 199     # triangular matrices.
 200     L_inverse = M.cholesky().inverse()
 201     return L_inverse.transpose()*L_inverse