]>
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 Every completely positive matrix is symmetric
26 positive-semidefinite::
28 sage: v = vector(map(abs, random_vector(ZZ, 10)))
29 sage: A = v.column() * v.row()
30 sage: is_symmetric_psd(A)
33 The following matrix is symmetric but not positive semidefinite::
35 sage: A = matrix(ZZ, [[1, 2], [2, 1]])
36 sage: is_symmetric_psd(A)
39 This matrix isn't even symmetric::
41 sage: A = matrix(ZZ, [[1, 2], [3, 4]])
42 sage: is_symmetric_psd(A)
47 if A
.base_ring() == SR
:
48 msg
= 'The matrix ``A`` cannot be symbolic.'
49 raise ValueError.new(msg
)
51 # First make sure that ``A`` is symmetric.
52 if not A
.is_symmetric():
55 # If ``A`` is symmetric, we only need to check that it is positive
56 # semidefinite. For that we can consult its minimum eigenvalue,
57 # which should be zero or greater. Since ``A`` is symmetric, its
58 # eigenvalues are guaranteed to be real.
59 return min(A
.eigenvalues()) >= 0
62 def unit_eigenvectors(A
):
64 Return the unit eigenvectors of a symmetric positive-definite matrix.
68 - ``A`` - The matrix whose eigenvectors we want to compute.
72 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
73 associated with its paired eigenvalue of ``A`` and has norm `1`.
77 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
78 sage: unit_evs = unit_eigenvectors(A)
79 sage: bool(unit_evs[0][1].norm() == 1)
81 sage: bool(unit_evs[1][1].norm() == 1)
83 sage: bool(unit_evs[2][1].norm() == 1)
87 # This will give us a list of lists whose elements are the
88 # eigenvectors we want.
89 ev_lists
= [ (val
,vecs
) for (val
,vecs
,multiplicity
)
90 in A
.eigenvectors_right() ]
92 # Pair each eigenvector with its eigenvalue and normalize it.
93 evs
= [ [(l
, vec
/vec
.norm()) for vec
in vecs
] for (l
,vecs
) in ev_lists
]
95 # Flatten the list, abusing the fact that "+" is overloaded on lists.
103 Factor a symmetric positive-semidefinite matrix ``A`` into
108 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
109 exact or the symbolic ring (to compute eigenvalues), and it
110 must be a field so that we can take its algebraic closure
111 (necessary to e.g. take square roots).
115 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
116 be the algebraic closure of the base field of ``A``.
120 Since ``A`` is symmetric and positive-semidefinite, we can
121 diagonalize it by some matrix `$Q$` whose columns are orthogonal
122 eigenvectors of ``A``. Then,
126 From this representation we can take the square root of `$D$`
127 (since all eigenvalues of ``A`` are nonnegative). If we then let
128 `$X = Q*sqrt(D)*Q^{T}$`, we have,
130 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
134 In principle, this is the algorithm used, although we ignore the
135 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
136 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
137 `$D$` will have dimension `$k \times k$`. In the end everything
142 Create a symmetric positive-semidefinite matrix over the symbolic
145 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
146 sage: X = factor_psd(A)
147 sage: A2 = (X*X.transpose()).simplify_full()
151 Attempt to factor the same matrix over ``RR`` which won't work
152 because ``RR`` isn't exact::
154 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
156 Traceback (most recent call last):
158 ValueError: The base ring of ``A`` must be either exact or symbolic.
160 Attempt to factor the same matrix over ``ZZ`` which won't work
161 because ``ZZ`` isn't a field::
163 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
165 Traceback (most recent call last):
167 ValueError: The base ring of ``A`` must be a field.
171 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
172 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
173 raise ValueError(msg
)
175 if not A
.base_ring().is_field():
176 raise ValueError('The base ring of ``A`` must be a field.')
178 if not A
.base_ring() is SR
:
179 # Change the base field of ``A`` so that we are sure we can take
180 # roots. The symbolic ring has no algebraic_closure method.
181 A
= A
.change_ring(A
.base_ring().algebraic_closure())
184 # Get the eigenvectors, and filter out the ones that correspond to
185 # the eigenvalue zero.
186 all_evs
= unit_eigenvectors(A
)
187 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
189 d
= [ sqrt(val
) for (val
,vec
) in evs
]
190 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
192 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
194 return Q
*root_D
*Q
.transpose()
197 def random_psd(V
, accept_zero
=True, rank
=None):
199 Generate a random symmetric positive-semidefinite matrix over the
200 vector space ``V``. That is, the returned matrix will be a linear
201 transformation on ``V``, with the same base ring as ``V``.
203 We take a very loose interpretation of "random," here. Otherwise we
204 would never (for example) choose a matrix on the boundary of the
205 cone (with a zero eigenvalue).
209 - ``V`` - The vector space on which the returned matrix will act.
211 - ``accept_zero`` - Do you want to accept the zero matrix (which
212 is symmetric PSD? Defaults to ``True``.
214 - ``rank`` - Require the returned matrix to have the given rank
219 A random symmetric positive semidefinite matrix, i.e. a linear
220 transformation from ``V`` to itself.
224 The matrix is constructed from some number of spectral projectors,
225 which in turn are created at "random" from the underlying vector
228 If no particular ``rank`` is desired, we choose the number of
229 projectors at random. Otherwise, we keep adding new projectors until
230 the desired rank is achieved.
232 Finally, before returning, we check if the matrix is zero. If
233 ``accept_zero`` is ``False``, we restart the process from the
238 Well, it doesn't crash at least::
240 sage: V = VectorSpace(QQ, 2)
241 sage: A = random_psd(V)
242 sage: A.matrix_space()
243 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
244 sage: is_symmetric_psd(A)
247 A matrix with the desired rank is returned::
249 sage: V = VectorSpace(QQ, 5)
250 sage: A = random_psd(V,False,1)
253 sage: A = random_psd(V,False,2)
256 sage: A = random_psd(V,False,3)
259 sage: A = random_psd(V,False,4)
262 sage: A = random_psd(V,False,5)
266 If the user asks for a rank that's too high, we fail::
268 sage: V = VectorSpace(QQ, 2)
269 sage: A = random_psd(V,False,3)
270 Traceback (most recent call last):
272 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
276 # We construct the matrix from its spectral projectors. Since
277 # there can be at most ``n`` of them, where ``n`` is the dimension
278 # of our vector space, we want to choose a random integer between
279 # ``0`` and ``n`` and then construct that many random elements of
285 # Choose one randomly
286 rank_A
= ZZ
.random_element(n
+1)
287 elif (rank
< 0) or (rank
> n
):
288 # The rank of ``A`` can be at most ``n``.
289 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
290 raise ValueError(msg
)
292 # Use the one the user gave us.
295 # Begin with the zero matrix, and add projectors to it if we have any.
296 A
= V
.zero().column()*V
.zero().row()
298 # Careful, begin at idx=1 so that we only generate a projector
299 # when rank_A is greater than zero.
300 while A
.rank() < rank_A
:
301 v
= V
.random_element()
302 A
+= v
.column()*v
.row()
304 if accept_zero
or not A
.is_zero():
305 # We either don't care what ``A`` is, or it's non-zero, so
309 # Uh oh, we need to generate a new one.
310 return random_psd(V
, accept_zero
, rank
)