]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/symmetric_pd.py
makefile: add --long to the test command.
[sage.d.git] / 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