]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/symmetric_psd.py
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}$`
9 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
10 # have to explicitly mangle our sitedir here so that "mjo.symbolic"
12 from os
.path
import abspath
13 from site
import addsitedir
14 addsitedir(abspath('../../'))
17 def is_symmetric_psd(A
):
19 Determine whether or not the matrix ``A`` is symmetric
20 positive-semidefinite.
24 - ``A`` - The matrix in question
28 Either ``True`` if ``A`` is symmetric positive-semidefinite, or
33 Every completely positive matrix is symmetric
34 positive-semidefinite::
36 sage: v = vector(map(abs, random_vector(ZZ, 10)))
37 sage: A = v.column() * v.row()
38 sage: is_symmetric_psd(A)
41 The following matrix is symmetric but not positive semidefinite::
43 sage: A = matrix(ZZ, [[1, 2], [2, 1]])
44 sage: is_symmetric_psd(A)
47 This matrix isn't even symmetric::
49 sage: A = matrix(ZZ, [[1, 2], [3, 4]])
50 sage: is_symmetric_psd(A)
55 if A
.base_ring() == SR
:
56 msg
= 'The matrix ``A`` cannot be symbolic.'
57 raise ValueError.new(msg
)
59 # First make sure that ``A`` is symmetric.
60 if not A
.is_symmetric():
63 # If ``A`` is symmetric, we only need to check that it is positive
64 # semidefinite. For that we can consult its minimum eigenvalue,
65 # which should be zero or greater. Since ``A`` is symmetric, its
66 # eigenvalues are guaranteed to be real.
67 return min(A
.eigenvalues()) >= 0
70 def unit_eigenvectors(A
):
72 Return the unit eigenvectors of a symmetric positive-definite matrix.
76 - ``A`` - The matrix whose eigenvectors we want to compute.
80 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
81 associated with its paired eigenvalue of ``A`` and has norm `1`.
85 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
86 sage: unit_evs = unit_eigenvectors(A)
87 sage: bool(unit_evs[0][1].norm() == 1)
89 sage: bool(unit_evs[1][1].norm() == 1)
91 sage: bool(unit_evs[2][1].norm() == 1)
95 # This will give us a list of lists whose elements are the
96 # eigenvectors we want.
97 ev_lists
= [ (val
,vecs
) for (val
,vecs
,multiplicity
)
98 in A
.eigenvectors_right() ]
100 # Pair each eigenvector with its eigenvalue and normalize it.
101 evs
= [ [(l
, vec
/vec
.norm()) for vec
in vecs
] for (l
,vecs
) in ev_lists
]
103 # Flatten the list, abusing the fact that "+" is overloaded on lists.
111 Factor a symmetric positive-semidefinite matrix ``A`` into
116 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
117 exact or the symbolic ring (to compute eigenvalues), and it
118 must be a field so that we can take its algebraic closure
119 (necessary to e.g. take square roots).
123 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
124 be the algebraic closure of the base field of ``A``.
128 Since ``A`` is symmetric and positive-semidefinite, we can
129 diagonalize it by some matrix `$Q$` whose columns are orthogonal
130 eigenvectors of ``A``. Then,
134 From this representation we can take the square root of `$D$`
135 (since all eigenvalues of ``A`` are nonnegative). If we then let
136 `$X = Q*sqrt(D)*Q^{T}$`, we have,
138 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
142 In principle, this is the algorithm used, although we ignore the
143 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
144 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
145 `$D$` will have dimension `$k \times k$`. In the end everything
150 Create a symmetric positive-semidefinite matrix over the symbolic
153 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
154 sage: X = factor_psd(A)
155 sage: A2 = (X*X.transpose()).simplify_full()
159 Attempt to factor the same matrix over ``RR`` which won't work
160 because ``RR`` isn't exact::
162 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
164 Traceback (most recent call last):
166 ValueError: The base ring of ``A`` must be either exact or symbolic.
168 Attempt to factor the same matrix over ``ZZ`` which won't work
169 because ``ZZ`` isn't a field::
171 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
173 Traceback (most recent call last):
175 ValueError: The base ring of ``A`` must be a field.
179 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
180 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
181 raise ValueError(msg
)
183 if not A
.base_ring().is_field():
184 raise ValueError('The base ring of ``A`` must be a field.')
186 if not A
.base_ring() is SR
:
187 # Change the base field of ``A`` so that we are sure we can take
188 # roots. The symbolic ring has no algebraic_closure method.
189 A
= A
.change_ring(A
.base_ring().algebraic_closure())
192 # Get the eigenvectors, and filter out the ones that correspond to
193 # the eigenvalue zero.
194 all_evs
= unit_eigenvectors(A
)
195 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
197 d
= [ sqrt(val
) for (val
,vec
) in evs
]
198 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
200 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
202 return Q
*root_D
*Q
.transpose()
205 def random_psd(V
, accept_zero
=True, rank
=None):
207 Generate a random symmetric positive-semidefinite matrix over the
208 vector space ``V``. That is, the returned matrix will be a linear
209 transformation on ``V``, with the same base ring as ``V``.
211 We take a very loose interpretation of "random," here. Otherwise we
212 would never (for example) choose a matrix on the boundary of the
213 cone (with a zero eigenvalue).
217 - ``V`` - The vector space on which the returned matrix will act.
219 - ``accept_zero`` - Do you want to accept the zero matrix (which
220 is symmetric PSD? Defaults to ``True``.
222 - ``rank`` - Require the returned matrix to have the given rank
227 A random symmetric positive semidefinite matrix, i.e. a linear
228 transformation from ``V`` to itself.
232 The matrix is constructed from some number of spectral projectors,
233 which in turn are created at "random" from the underlying vector
236 If no particular ``rank`` is desired, we choose the number of
237 projectors at random. Otherwise, we keep adding new projectors until
238 the desired rank is achieved.
240 Finally, before returning, we check if the matrix is zero. If
241 ``accept_zero`` is ``False``, we restart the process from the
246 Well, it doesn't crash at least::
248 sage: V = VectorSpace(QQ, 2)
249 sage: A = random_psd(V)
250 sage: A.matrix_space()
251 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
252 sage: is_symmetric_psd(A)
255 A matrix with the desired rank is returned::
257 sage: V = VectorSpace(QQ, 5)
258 sage: A = random_psd(V,False,1)
261 sage: A = random_psd(V,False,2)
264 sage: A = random_psd(V,False,3)
267 sage: A = random_psd(V,False,4)
270 sage: A = random_psd(V,False,5)
274 If the user asks for a rank that's too high, we fail::
276 sage: V = VectorSpace(QQ, 2)
277 sage: A = random_psd(V,False,3)
278 Traceback (most recent call last):
280 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
284 # We construct the matrix from its spectral projectors. Since
285 # there can be at most ``n`` of them, where ``n`` is the dimension
286 # of our vector space, we want to choose a random integer between
287 # ``0`` and ``n`` and then construct that many random elements of
293 # Choose one randomly
294 rank_A
= ZZ
.random_element(n
+1)
295 elif (rank
< 0) or (rank
> n
):
296 # The rank of ``A`` can be at most ``n``.
297 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
298 raise ValueError(msg
)
300 # Use the one the user gave us.
303 # Begin with the zero matrix, and add projectors to it if we have
305 A
= V
.zero_element().column()*V
.zero_element().row()
307 # Careful, begin at idx=1 so that we only generate a projector
308 # when rank_A is greater than zero.
309 while A
.rank() < rank_A
:
310 v
= V
.random_element()
311 A
+= v
.column()*v
.row()
313 if accept_zero
or not A
.is_zero():
314 # We either don't care what ``A`` is, or it's non-zero, so
318 # Uh oh, we need to generate a new one.
319 return random_psd(V
, accept_zero
, rank
)