]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/symmetric_psd.py
68be78fbdd53630303e5b409727fa686d8ad7e32
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 unit 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`. If
78 the base ring of ``A`` is not algebraically closed, then returned
79 eigenvectors may (necessarily) be over its algebraic closure and not
80 the base ring of ``A`` itself.
84 sage: from mjo.cone.symmetric_psd import unit_eigenvectors
88 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
89 sage: unit_evs = list(unit_eigenvectors(A))
90 sage: bool(unit_evs[0][1].norm() == 1)
92 sage: bool(unit_evs[1][1].norm() == 1)
94 sage: bool(unit_evs[2][1].norm() == 1)
98 return ( (val
,vec
.normalized())
99 for (val
,vecs
,multiplicity
) in A
.eigenvectors_right()
106 Factor a symmetric positive-semidefinite matrix ``A`` into
111 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
112 exact or the symbolic ring (to compute eigenvalues), and it
113 must be a field so that we can take its algebraic closure
114 (necessary to e.g. take square roots).
118 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
119 be the algebraic closure of the base field of ``A``.
123 Since ``A`` is symmetric and positive-semidefinite, we can
124 diagonalize it by some matrix `$Q$` whose columns are orthogonal
125 eigenvectors of ``A``. Then,
129 From this representation we can take the square root of `$D$`
130 (since all eigenvalues of ``A`` are nonnegative). If we then let
131 `$X = Q*sqrt(D)*Q^{T}$`, we have,
133 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
137 In principle, this is the algorithm used, although we ignore the
138 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
139 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
140 `$D$` will have dimension `$k \times k$`. In the end everything
145 sage: from mjo.cone.symmetric_psd import factor_psd
149 Create a symmetric positive-semidefinite matrix over the symbolic
152 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
153 sage: X = factor_psd(A)
154 sage: A2 = (X*X.transpose()).simplify_full()
158 Attempt to factor the same matrix over ``RR`` which won't work
159 because ``RR`` isn't exact::
161 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
163 Traceback (most recent call last):
165 ValueError: The base ring of ``A`` must be either exact or symbolic.
167 Attempt to factor the same matrix over ``ZZ`` which won't work
168 because ``ZZ`` isn't a field::
170 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
172 Traceback (most recent call last):
174 ValueError: The base ring of ``A`` must be a field.
178 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
179 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
180 raise ValueError(msg
)
182 if not A
.base_ring().is_field():
183 raise ValueError('The base ring of ``A`` must be a field.')
185 if not A
.base_ring() is SR
:
186 # Change the base field of ``A`` so that we are sure we can take
187 # roots. The symbolic ring has no algebraic_closure method.
188 A
= A
.change_ring(A
.base_ring().algebraic_closure())
191 # Get the eigenvectors, and filter out the ones that correspond to
192 # the eigenvalue zero.
193 all_evs
= unit_eigenvectors(A
)
194 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
196 d
= [ sqrt(val
) for (val
,vec
) in evs
]
197 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
199 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
201 return Q
*root_D
*Q
.transpose()
204 def random_psd(V
, accept_zero
=True, rank
=None):
206 Generate a random symmetric positive-semidefinite matrix over the
207 vector space ``V``. That is, the returned matrix will be a linear
208 transformation on ``V``, with the same base ring as ``V``.
210 We take a very loose interpretation of "random," here. Otherwise we
211 would never (for example) choose a matrix on the boundary of the
212 cone (with a zero eigenvalue).
216 - ``V`` - The vector space on which the returned matrix will act.
218 - ``accept_zero`` - Do you want to accept the zero matrix (which
219 is symmetric PSD? Defaults to ``True``.
221 - ``rank`` - Require the returned matrix to have the given rank
226 A random symmetric positive semidefinite matrix, i.e. a linear
227 transformation from ``V`` to itself.
231 The matrix is constructed from some number of spectral projectors,
232 which in turn are created at "random" from the underlying vector
235 If no particular ``rank`` is desired, we choose the number of
236 projectors at random. Otherwise, we keep adding new projectors until
237 the desired rank is achieved.
239 Finally, before returning, we check if the matrix is zero. If
240 ``accept_zero`` is ``False``, we restart the process from the
245 sage: from mjo.cone.symmetric_psd import is_symmetric_psd, random_psd
249 Well, it doesn't crash at least::
251 sage: V = VectorSpace(QQ, 2)
252 sage: A = random_psd(V)
253 sage: A.matrix_space()
254 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
255 sage: is_symmetric_psd(A)
258 A matrix with the desired rank is returned::
260 sage: V = VectorSpace(QQ, 5)
261 sage: A = random_psd(V,False,1)
264 sage: A = random_psd(V,False,2)
267 sage: A = random_psd(V,False,3)
270 sage: A = random_psd(V,False,4)
273 sage: A = random_psd(V,False,5)
277 If the user asks for a rank that's too high, we fail::
279 sage: V = VectorSpace(QQ, 2)
280 sage: A = random_psd(V,False,3)
281 Traceback (most recent call last):
283 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
287 # We construct the matrix from its spectral projectors. Since
288 # there can be at most ``n`` of them, where ``n`` is the dimension
289 # of our vector space, we want to choose a random integer between
290 # ``0`` and ``n`` and then construct that many random elements of
296 # Choose one randomly
297 rank_A
= ZZ
.random_element(n
+1)
298 elif (rank
< 0) or (rank
> n
):
299 # The rank of ``A`` can be at most ``n``.
300 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
301 raise ValueError(msg
)
303 # Use the one the user gave us.
306 # Begin with the zero matrix, and add projectors to it if we have any.
307 A
= V
.zero().column()*V
.zero().row()
309 # Careful, begin at idx=1 so that we only generate a projector
310 # when rank_A is greater than zero.
311 while A
.rank() < rank_A
:
312 v
= V
.random_element()
313 A
+= v
.column()*v
.row()
315 if accept_zero
or not A
.is_zero():
316 # We either don't care what ``A`` is, or it's non-zero, so
320 # Uh oh, we need to generate a new one.
321 return random_psd(V
, accept_zero
, rank
)