]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/symmetric_psd.py
makefile: add --long to the test command.
[sage.d.git] / mjo / cone / symmetric_psd.py
1 r"""
2 The positive semidefinite cone `$S^{n}_{+}$` is the cone consisting of
3 all symmetric positive-semidefinite matrices (as a subset of
4 `$\mathbb{R}^{n \times n}$`
5 """
6
7 from sage.all import *
8
9 def is_symmetric_psd(A):
10 """
11 Determine whether or not the matrix ``A`` is symmetric
12 positive-semidefinite.
13
14 INPUT:
15
16 - ``A`` - The matrix in question
17
18 OUTPUT:
19
20 Either ``True`` if ``A`` is symmetric positive-semidefinite, or
21 ``False`` otherwise.
22
23 SETUP::
24
25 sage: from mjo.cone.symmetric_psd import is_symmetric_psd
26
27 EXAMPLES:
28
29 Every completely positive matrix is symmetric
30 positive-semidefinite::
31
32 sage: set_random_seed()
33 sage: v = vector(map(abs, random_vector(ZZ, 10)))
34 sage: A = v.column() * v.row()
35 sage: is_symmetric_psd(A)
36 True
37
38 The following matrix is symmetric but not positive semidefinite::
39
40 sage: A = matrix(ZZ, [[1, 2], [2, 1]])
41 sage: is_symmetric_psd(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_psd(A)
48 False
49
50 The trivial matrix in a trivial space is trivially symmetric and
51 positive-semidefinite::
52
53 sage: A = matrix(QQ, 0,0)
54 sage: is_symmetric_psd(A)
55 True
56
57 """
58
59 if A.base_ring() == SR:
60 msg = 'The matrix ``A`` cannot be symbolic.'
61 raise ValueError.new(msg)
62
63 # First make sure that ``A`` is symmetric.
64 if not A.is_symmetric():
65 return False
66
67 # If ``A`` is symmetric, we only need to check that it is positive
68 # semidefinite. For that we can consult its minimum eigenvalue,
69 # which should be zero or greater. Since ``A`` is symmetric, its
70 # eigenvalues are guaranteed to be real.
71 if A.is_zero():
72 # A is trivial... so trivially positive-semudefinite.
73 return True
74 else:
75 return min(A.eigenvalues()) >= 0
76
77
78 def unit_eigenvectors(A):
79 """
80 Return the unit eigenvectors of a symmetric positive-definite matrix.
81
82 INPUT:
83
84 - ``A`` -- The matrix whose unit eigenvectors we want to compute.
85
86 OUTPUT:
87
88 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
89 associated with its paired eigenvalue of ``A`` and has norm `1`. If
90 the base ring of ``A`` is not algebraically closed, then returned
91 eigenvectors may (necessarily) be over its algebraic closure and not
92 the base ring of ``A`` itself.
93
94 SETUP::
95
96 sage: from mjo.cone.symmetric_psd import unit_eigenvectors
97
98 EXAMPLES::
99
100 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
101 sage: unit_evs = list(unit_eigenvectors(A))
102 sage: bool(unit_evs[0][1].norm() == 1)
103 True
104 sage: bool(unit_evs[1][1].norm() == 1)
105 True
106 sage: bool(unit_evs[2][1].norm() == 1)
107 True
108
109 """
110 return ( (val,vec.normalized())
111 for (val,vecs,multiplicity) in A.eigenvectors_right()
112 for vec in vecs )
113
114
115
116 def factor_psd(A):
117 """
118 Factor a symmetric positive-semidefinite matrix ``A`` into
119 `XX^{T}`.
120
121 INPUT:
122
123 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
124 exact or the symbolic ring (to compute eigenvalues), and it
125 must be a field so that we can take its algebraic closure
126 (necessary to e.g. take square roots).
127
128 OUTPUT:
129
130 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
131 be the algebraic closure of the base field of ``A``.
132
133 ALGORITHM:
134
135 Since ``A`` is symmetric and positive-semidefinite, we can
136 diagonalize it by some matrix `$Q$` whose columns are orthogonal
137 eigenvectors of ``A``. Then,
138
139 `$A = QDQ^{T}$`
140
141 From this representation we can take the square root of `$D$`
142 (since all eigenvalues of ``A`` are nonnegative). If we then let
143 `$X = Q*sqrt(D)*Q^{T}$`, we have,
144
145 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
146
147 as desired.
148
149 In principle, this is the algorithm used, although we ignore the
150 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
151 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
152 `$D$` will have dimension `$k \times k$`. In the end everything
153 works out the same.
154
155 SETUP::
156
157 sage: from mjo.cone.symmetric_psd import factor_psd
158
159 EXAMPLES:
160
161 Create a symmetric positive-semidefinite matrix over the symbolic
162 ring and factor it::
163
164 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
165 sage: X = factor_psd(A)
166 sage: A2 = (X*X.transpose()).simplify_full()
167 sage: A == A2
168 True
169
170 Attempt to factor the same matrix over ``RR`` which won't work
171 because ``RR`` isn't exact::
172
173 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
174 sage: factor_psd(A)
175 Traceback (most recent call last):
176 ...
177 ValueError: The base ring of ``A`` must be either exact or symbolic.
178
179 Attempt to factor the same matrix over ``ZZ`` which won't work
180 because ``ZZ`` isn't a field::
181
182 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
183 sage: factor_psd(A)
184 Traceback (most recent call last):
185 ...
186 ValueError: The base ring of ``A`` must be a field.
187
188 """
189
190 if not A.base_ring().is_exact() and not A.base_ring() is SR:
191 msg = 'The base ring of ``A`` must be either exact or symbolic.'
192 raise ValueError(msg)
193
194 if not A.base_ring().is_field():
195 raise ValueError('The base ring of ``A`` must be a field.')
196
197 if not A.base_ring() is SR:
198 # Change the base field of ``A`` so that we are sure we can take
199 # roots. The symbolic ring has no algebraic_closure method.
200 A = A.change_ring(A.base_ring().algebraic_closure())
201
202
203 # Get the eigenvectors, and filter out the ones that correspond to
204 # the eigenvalue zero.
205 all_evs = unit_eigenvectors(A)
206 evs = [ (val,vec) for (val,vec) in all_evs if not val == 0 ]
207
208 d = ( val.sqrt() for (val,vec) in evs )
209 root_D = diagonal_matrix(d).change_ring(A.base_ring())
210
211 Q = matrix(A.base_ring(), ( vec for (val,vec) in evs )).transpose()
212
213 return Q*root_D*Q.transpose()
214
215
216 def random_symmetric_psd(V, accept_zero=True, rank=None):
217 """
218 Generate a random symmetric positive-semidefinite matrix over the
219 vector space ``V``. That is, the returned matrix will be a linear
220 transformation on ``V``, with the same base ring as ``V``.
221
222 We take a very loose interpretation of "random," here. Otherwise we
223 would never (for example) choose a matrix on the boundary of the
224 cone (with a zero eigenvalue).
225
226 INPUT:
227
228 - ``V`` - The vector space on which the returned matrix will act.
229
230 - ``accept_zero`` - Do you want to accept the zero matrix (which
231 is symmetric PSD? Defaults to ``True``.
232
233 - ``rank`` - Require the returned matrix to have the given rank
234 (optional).
235
236 OUTPUT:
237
238 A random symmetric positive semidefinite matrix, i.e. a linear
239 transformation from ``V`` to itself.
240
241 ALGORITHM:
242
243 The matrix is constructed from some number of spectral projectors,
244 which in turn are created at "random" from the underlying vector
245 space ``V``.
246
247 If no particular ``rank`` is desired, we choose the number of
248 projectors at random. Otherwise, we keep adding new projectors until
249 the desired rank is achieved.
250
251 Finally, before returning, we check if the matrix is zero. If
252 ``accept_zero`` is ``False``, we restart the process from the
253 beginning.
254
255 SETUP::
256
257 sage: from mjo.cone.symmetric_psd import (is_symmetric_psd,
258 ....: random_symmetric_psd)
259
260 EXAMPLES:
261
262 Well, it doesn't crash at least::
263
264 sage: set_random_seed()
265 sage: V = VectorSpace(QQ, 2)
266 sage: A = random_symmetric_psd(V)
267 sage: A.matrix_space()
268 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
269 sage: is_symmetric_psd(A)
270 True
271
272 A matrix with the desired rank is returned::
273
274 sage: set_random_seed()
275 sage: V = VectorSpace(QQ, 5)
276 sage: A = random_symmetric_psd(V,False,1)
277 sage: A.rank()
278 1
279 sage: A = random_symmetric_psd(V,False,2)
280 sage: A.rank()
281 2
282 sage: A = random_symmetric_psd(V,False,3)
283 sage: A.rank()
284 3
285 sage: A = random_symmetric_psd(V,False,4)
286 sage: A.rank()
287 4
288 sage: A = random_symmetric_psd(V,False,5)
289 sage: A.rank()
290 5
291
292 If the user asks for a rank that's too high, we fail::
293
294 sage: set_random_seed()
295 sage: V = VectorSpace(QQ, 2)
296 sage: A = random_symmetric_psd(V,False,3)
297 Traceback (most recent call last):
298 ...
299 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
300
301 """
302
303 # We construct the matrix from its spectral projectors. Since
304 # there can be at most ``n`` of them, where ``n`` is the dimension
305 # of our vector space, we want to choose a random integer between
306 # ``0`` and ``n`` and then construct that many random elements of
307 # ``V``.
308 n = V.dimension()
309
310 rank_A = 0
311 if rank is None:
312 # Choose one randomly
313 rank_A = ZZ.random_element(n+1)
314 elif (rank < 0) or (rank > n):
315 # The rank of ``A`` can be at most ``n``.
316 msg = 'The ``rank`` must be between 0 and the dimension of ``V``.'
317 raise ValueError(msg)
318 else:
319 # Use the one the user gave us.
320 rank_A = rank
321
322 if n == 0 and not accept_zero:
323 # We're gonna loop forever trying to satisfy this...
324 raise ValueError('You must have accept_zero=True when V is trivial')
325
326 # Loop until we find a suitable "A" that will then be returned.
327 while True:
328 # Begin with the zero matrix, and add projectors to it if we
329 # have any.
330 A = matrix.zero(V.base_ring(), n, n)
331
332 # Careful, begin at idx=1 so that we only generate a projector
333 # when rank_A is greater than zero.
334 while A.rank() < rank_A:
335 v = V.random_element()
336 A += v.column()*v.row()
337
338 if accept_zero or not A.is_zero():
339 # We either don't care what ``A`` is, or it's non-zero, so
340 # just return it.
341 return A