]>
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 def is_symmetric_psd(A
):
11 Determine whether or not the matrix ``A`` is symmetric
12 positive-semidefinite.
16 - ``A`` - The matrix in question
20 Either ``True`` if ``A`` is symmetric positive-semidefinite, or
25 sage: from mjo.cone.symmetric_psd import is_symmetric_psd
29 Every completely positive matrix is symmetric
30 positive-semidefinite::
32 sage: v = vector(map(abs, random_vector(ZZ, 10)))
33 sage: A = v.column() * v.row()
34 sage: is_symmetric_psd(A)
37 The following matrix is symmetric but not positive semidefinite::
39 sage: A = matrix(ZZ, [[1, 2], [2, 1]])
40 sage: is_symmetric_psd(A)
43 This matrix isn't even symmetric::
45 sage: A = matrix(ZZ, [[1, 2], [3, 4]])
46 sage: is_symmetric_psd(A)
51 if A
.base_ring() == SR
:
52 msg
= 'The matrix ``A`` cannot be symbolic.'
53 raise ValueError.new(msg
)
55 # First make sure that ``A`` is symmetric.
56 if not A
.is_symmetric():
59 # If ``A`` is symmetric, we only need to check that it is positive
60 # semidefinite. For that we can consult its minimum eigenvalue,
61 # which should be zero or greater. Since ``A`` is symmetric, its
62 # eigenvalues are guaranteed to be real.
63 return min(A
.eigenvalues()) >= 0
66 def unit_eigenvectors(A
):
68 Return the unit eigenvectors of a symmetric positive-definite matrix.
72 - ``A`` - The matrix whose eigenvectors we want to compute.
76 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
77 associated with its paired eigenvalue of ``A`` and has norm `1`.
81 sage: from mjo.cone.symmetric_psd import unit_eigenvectors
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 sage: from mjo.cone.symmetric_psd import factor_psd
154 Create a symmetric positive-semidefinite matrix over the symbolic
157 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
158 sage: X = factor_psd(A)
159 sage: A2 = (X*X.transpose()).simplify_full()
163 Attempt to factor the same matrix over ``RR`` which won't work
164 because ``RR`` isn't exact::
166 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
168 Traceback (most recent call last):
170 ValueError: The base ring of ``A`` must be either exact or symbolic.
172 Attempt to factor the same matrix over ``ZZ`` which won't work
173 because ``ZZ`` isn't a field::
175 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
177 Traceback (most recent call last):
179 ValueError: The base ring of ``A`` must be a field.
183 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
184 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
185 raise ValueError(msg
)
187 if not A
.base_ring().is_field():
188 raise ValueError('The base ring of ``A`` must be a field.')
190 if not A
.base_ring() is SR
:
191 # Change the base field of ``A`` so that we are sure we can take
192 # roots. The symbolic ring has no algebraic_closure method.
193 A
= A
.change_ring(A
.base_ring().algebraic_closure())
196 # Get the eigenvectors, and filter out the ones that correspond to
197 # the eigenvalue zero.
198 all_evs
= unit_eigenvectors(A
)
199 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
201 d
= [ sqrt(val
) for (val
,vec
) in evs
]
202 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
204 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
206 return Q
*root_D
*Q
.transpose()
209 def random_psd(V
, accept_zero
=True, rank
=None):
211 Generate a random symmetric positive-semidefinite matrix over the
212 vector space ``V``. That is, the returned matrix will be a linear
213 transformation on ``V``, with the same base ring as ``V``.
215 We take a very loose interpretation of "random," here. Otherwise we
216 would never (for example) choose a matrix on the boundary of the
217 cone (with a zero eigenvalue).
221 - ``V`` - The vector space on which the returned matrix will act.
223 - ``accept_zero`` - Do you want to accept the zero matrix (which
224 is symmetric PSD? Defaults to ``True``.
226 - ``rank`` - Require the returned matrix to have the given rank
231 A random symmetric positive semidefinite matrix, i.e. a linear
232 transformation from ``V`` to itself.
236 The matrix is constructed from some number of spectral projectors,
237 which in turn are created at "random" from the underlying vector
240 If no particular ``rank`` is desired, we choose the number of
241 projectors at random. Otherwise, we keep adding new projectors until
242 the desired rank is achieved.
244 Finally, before returning, we check if the matrix is zero. If
245 ``accept_zero`` is ``False``, we restart the process from the
250 sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd
254 Well, it doesn't crash at least::
256 sage: V = VectorSpace(QQ, 2)
257 sage: A = random_psd(V)
258 sage: A.matrix_space()
259 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
260 sage: is_symmetric_psd(A)
263 A matrix with the desired rank is returned::
265 sage: V = VectorSpace(QQ, 5)
266 sage: A = random_psd(V,False,1)
269 sage: A = random_psd(V,False,2)
272 sage: A = random_psd(V,False,3)
275 sage: A = random_psd(V,False,4)
278 sage: A = random_psd(V,False,5)
282 If the user asks for a rank that's too high, we fail::
284 sage: V = VectorSpace(QQ, 2)
285 sage: A = random_psd(V,False,3)
286 Traceback (most recent call last):
288 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
292 # We construct the matrix from its spectral projectors. Since
293 # there can be at most ``n`` of them, where ``n`` is the dimension
294 # of our vector space, we want to choose a random integer between
295 # ``0`` and ``n`` and then construct that many random elements of
301 # Choose one randomly
302 rank_A
= ZZ
.random_element(n
+1)
303 elif (rank
< 0) or (rank
> n
):
304 # The rank of ``A`` can be at most ``n``.
305 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
306 raise ValueError(msg
)
308 # Use the one the user gave us.
311 # Begin with the zero matrix, and add projectors to it if we have any.
312 A
= V
.zero().column()*V
.zero().row()
314 # Careful, begin at idx=1 so that we only generate a projector
315 # when rank_A is greater than zero.
316 while A
.rank() < rank_A
:
317 v
= V
.random_element()
318 A
+= v
.column()*v
.row()
320 if accept_zero
or not A
.is_zero():
321 # We either don't care what ``A`` is, or it's non-zero, so
325 # Uh oh, we need to generate a new one.
326 return random_psd(V
, accept_zero
, rank
)