]>
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('../../'))
15 from mjo
.symbolic
import matrix_simplify_full
18 def is_symmetric_psd(A
):
20 Determine whether or not the matrix ``A`` is symmetric
21 positive-semidefinite.
25 - ``A`` - The matrix in question
29 Either ``True`` if ``A`` is symmetric positive-semidefinite, or
34 Every completely positive matrix is symmetric
35 positive-semidefinite::
37 sage: v = vector(map(abs, random_vector(ZZ, 10)))
38 sage: A = v.column() * v.row()
39 sage: is_symmetric_psd(A)
42 The following matrix is symmetric but not positive semidefinite::
44 sage: A = matrix(ZZ, [[1, 2], [2, 1]])
45 sage: is_symmetric_psd(A)
48 This matrix isn't even symmetric::
50 sage: A = matrix(ZZ, [[1, 2], [3, 4]])
51 sage: is_symmetric_psd(A)
56 if A
.base_ring() == SR
:
57 msg
= 'The matrix ``A`` cannot be symbolic.'
58 raise ValueError.new(msg
)
60 # First make sure that ``A`` is symmetric.
61 if not A
.is_symmetric():
64 # If ``A`` is symmetric, we only need to check that it is positive
65 # semidefinite. For that we can consult its minimum eigenvalue,
66 # which should be zero or greater. Since ``A`` is symmetric, its
67 # eigenvalues are guaranteed to be real.
68 return min(A
.eigenvalues()) >= 0
71 def unit_eigenvectors(A
):
73 Return the unit eigenvectors of a symmetric positive-definite matrix.
77 - ``A`` - The matrix whose eigenvectors we want to compute.
81 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
82 associated with its paired eigenvalue of ``A`` and has norm `1`.
86 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
87 sage: unit_evs = unit_eigenvectors(A)
88 sage: bool(unit_evs[0][1].norm() == 1)
90 sage: bool(unit_evs[1][1].norm() == 1)
92 sage: bool(unit_evs[2][1].norm() == 1)
96 # This will give us a list of lists whose elements are the
97 # eigenvectors we want.
98 ev_lists
= [ (val
,vecs
) for (val
,vecs
,multiplicity
)
99 in A
.eigenvectors_right() ]
101 # Pair each eigenvector with its eigenvalue and normalize it.
102 evs
= [ [(l
, vec
/vec
.norm()) for vec
in vecs
] for (l
,vecs
) in ev_lists
]
104 # Flatten the list, abusing the fact that "+" is overloaded on lists.
112 Factor a symmetric positive-semidefinite matrix ``A`` into
117 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
118 exact or the symbolic ring (to compute eigenvalues), and it
119 must be a field so that we can take its algebraic closure
120 (necessary to e.g. take square roots).
124 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
125 be the algebraic closure of the base field of ``A``.
129 Since ``A`` is symmetric and positive-semidefinite, we can
130 diagonalize it by some matrix `$Q$` whose columns are orthogonal
131 eigenvectors of ``A``. Then,
135 From this representation we can take the square root of `$D$`
136 (since all eigenvalues of ``A`` are nonnegative). If we then let
137 `$X = Q*sqrt(D)*Q^{T}$`, we have,
139 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
143 In principle, this is the algorithm used, although we ignore the
144 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
145 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
146 `$D$` will have dimension `$k \times k$`. In the end everything
151 Create a symmetric positive-semidefinite matrix over the symbolic
154 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
155 sage: X = factor_psd(A)
156 sage: A2 = matrix_simplify_full(X*X.transpose())
160 Attempt to factor the same matrix over ``RR`` which won't work
161 because ``RR`` isn't exact::
163 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
165 Traceback (most recent call last):
167 ValueError: The base ring of ``A`` must be either exact or symbolic.
169 Attempt to factor the same matrix over ``ZZ`` which won't work
170 because ``ZZ`` isn't a field::
172 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
174 Traceback (most recent call last):
176 ValueError: The base ring of ``A`` must be a field.
180 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
181 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
182 raise ValueError(msg
)
184 if not A
.base_ring().is_field():
185 raise ValueError('The base ring of ``A`` must be a field.')
187 if not A
.base_ring() is SR
:
188 # Change the base field of ``A`` so that we are sure we can take
189 # roots. The symbolic ring has no algebraic_closure method.
190 A
= A
.change_ring(A
.base_ring().algebraic_closure())
193 # Get the eigenvectors, and filter out the ones that correspond to
194 # the eigenvalue zero.
195 all_evs
= unit_eigenvectors(A
)
196 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
198 d
= [ sqrt(val
) for (val
,vec
) in evs
]
199 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
201 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
203 return Q
*root_D
*Q
.transpose()
206 def random_psd(V
, accept_zero
=True, rank
=None):
208 Generate a random symmetric positive-semidefinite matrix over the
209 vector space ``V``. That is, the returned matrix will be a linear
210 transformation on ``V``, with the same base ring as ``V``.
212 We take a very loose interpretation of "random," here. Otherwise we
213 would never (for example) choose a matrix on the boundary of the
214 cone (with a zero eigenvalue).
218 - ``V`` - The vector space on which the returned matrix will act.
220 - ``accept_zero`` - Do you want to accept the zero matrix (which
221 is symmetric PSD? Defaults to ``True``.
223 - ``rank`` - Require the returned matrix to have the given rank
228 A random symmetric positive semidefinite matrix, i.e. a linear
229 transformation from ``V`` to itself.
233 The matrix is constructed from some number of spectral projectors,
234 which in turn are created at "random" from the underlying vector
237 If no particular ``rank`` is desired, we choose the number of
238 projectors at random. Otherwise, we keep adding new projectors until
239 the desired rank is achieved.
241 Finally, before returning, we check if the matrix is zero. If
242 ``accept_zero`` is ``False``, we restart the process from the
247 Well, it doesn't crash at least::
249 sage: V = VectorSpace(QQ, 2)
250 sage: A = random_psd(V)
251 sage: A.matrix_space()
252 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
253 sage: is_symmetric_psd(A)
256 A matrix with the desired rank is returned::
258 sage: V = VectorSpace(QQ, 5)
259 sage: A = random_psd(V,False,1)
262 sage: A = random_psd(V,False,2)
265 sage: A = random_psd(V,False,3)
268 sage: A = random_psd(V,False,4)
271 sage: A = random_psd(V,False,5)
275 If the user asks for a rank that's too high, we fail::
277 sage: V = VectorSpace(QQ, 2)
278 sage: A = random_psd(V,False,3)
279 Traceback (most recent call last):
281 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
285 # We construct the matrix from its spectral projectors. Since
286 # there can be at most ``n`` of them, where ``n`` is the dimension
287 # of our vector space, we want to choose a random integer between
288 # ``0`` and ``n`` and then construct that many random elements of
294 # Choose one randomly
295 rank_A
= ZZ
.random_element(n
+1)
296 elif (rank
< 0) or (rank
> n
):
297 # The rank of ``A`` can be at most ``n``.
298 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
299 raise ValueError(msg
)
301 # Use the one the user gave us.
304 # Begin with the zero matrix, and add projectors to it if we have
306 A
= V
.zero_element().column()*V
.zero_element().row()
308 # Careful, begin at idx=1 so that we only generate a projector
309 # when rank_A is greater than zero.
310 while A
.rank() < rank_A
:
311 v
= V
.random_element()
312 A
+= v
.column()*v
.row()
314 if accept_zero
or not A
.is_zero():
315 # We either don't care what ``A`` is, or it's non-zero, so
319 # Uh oh, we need to generate a new one.
320 return random_psd(V
, accept_zero
, rank
)