]>
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 unit_eigenvectors(A
):
11 Return the unit eigenvectors of a symmetric positive-definite matrix.
15 - ``A`` -- The matrix whose unit eigenvectors we want to compute.
19 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
20 associated with its paired eigenvalue of ``A`` and has norm `1`. If
21 the base ring of ``A`` is not algebraically closed, then returned
22 eigenvectors may (necessarily) be over its algebraic closure and not
23 the base ring of ``A`` itself.
27 sage: from mjo.cone.symmetric_psd import unit_eigenvectors
31 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
32 sage: unit_evs = list(unit_eigenvectors(A))
33 sage: bool(unit_evs[0][1].norm() == 1)
35 sage: bool(unit_evs[1][1].norm() == 1)
37 sage: bool(unit_evs[2][1].norm() == 1)
41 return ( (val
,vec
.normalized())
42 for (val
,vecs
,multiplicity
) in A
.eigenvectors_right()
49 Factor a symmetric positive-semidefinite matrix ``A`` into
54 - ``A`` - The matrix to factor. The base ring of ``A`` must be either
55 exact or the symbolic ring (to compute eigenvalues), and it
56 must be a field so that we can take its algebraic closure
57 (necessary to e.g. take square roots).
61 A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
62 be the algebraic closure of the base field of ``A``.
66 Since ``A`` is symmetric and positive-semidefinite, we can
67 diagonalize it by some matrix `$Q$` whose columns are orthogonal
68 eigenvectors of ``A``. Then,
72 From this representation we can take the square root of `$D$`
73 (since all eigenvalues of ``A`` are nonnegative). If we then let
74 `$X = Q*sqrt(D)*Q^{T}$`, we have,
76 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
80 In principle, this is the algorithm used, although we ignore the
81 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
82 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
83 `$D$` will have dimension `$k \times k$`. In the end everything
88 sage: from mjo.cone.symmetric_psd import factor_psd
92 Create a symmetric positive-semidefinite matrix over the symbolic
95 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
96 sage: X = factor_psd(A)
97 sage: A2 = (X*X.transpose()).simplify_full()
101 Attempt to factor the same matrix over ``RR`` which won't work
102 because ``RR`` isn't exact::
104 sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
106 Traceback (most recent call last):
108 ValueError: The base ring of ``A`` must be either exact or symbolic.
110 Attempt to factor the same matrix over ``ZZ`` which won't work
111 because ``ZZ`` isn't a field::
113 sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
115 Traceback (most recent call last):
117 ValueError: The base ring of ``A`` must be a field.
121 if not A
.base_ring().is_exact() and not A
.base_ring() is SR
:
122 msg
= 'The base ring of ``A`` must be either exact or symbolic.'
123 raise ValueError(msg
)
125 if not A
.base_ring().is_field():
126 raise ValueError('The base ring of ``A`` must be a field.')
128 if not A
.base_ring() is SR
:
129 # Change the base field of ``A`` so that we are sure we can take
130 # roots. The symbolic ring has no algebraic_closure method.
131 A
= A
.change_ring(A
.base_ring().algebraic_closure())
134 # Get the eigenvectors, and filter out the ones that correspond to
135 # the eigenvalue zero.
136 all_evs
= unit_eigenvectors(A
)
137 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
139 d
= ( val
.sqrt() for (val
,vec
) in evs
)
140 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
142 Q
= matrix(A
.base_ring(), ( vec
for (val
,vec
) in evs
)).transpose()
144 return Q
*root_D
*Q
.transpose()
147 def random_symmetric_psd(V
, accept_zero
=True, rank
=None):
149 Generate a random symmetric positive-semidefinite matrix over the
150 vector space ``V``. That is, the returned matrix will be a linear
151 transformation on ``V``, with the same base ring as ``V``.
153 We take a very loose interpretation of "random," here. Otherwise we
154 would never (for example) choose a matrix on the boundary of the
155 cone (with a zero eigenvalue).
159 - ``V`` - The vector space on which the returned matrix will act.
161 - ``accept_zero`` - Do you want to accept the zero matrix (which
162 is symmetric PSD? Defaults to ``True``.
164 - ``rank`` - Require the returned matrix to have the given rank
169 A random symmetric positive semidefinite matrix, i.e. a linear
170 transformation from ``V`` to itself.
174 The matrix is constructed from some number of spectral projectors,
175 which in turn are created at "random" from the underlying vector
178 If no particular ``rank`` is desired, we choose the number of
179 projectors at random. Otherwise, we keep adding new projectors until
180 the desired rank is achieved.
182 Finally, before returning, we check if the matrix is zero. If
183 ``accept_zero`` is ``False``, we restart the process from the
188 sage: from mjo.cone.symmetric_psd import random_symmetric_psd
192 Well, it doesn't crash at least::
194 sage: V = VectorSpace(QQ, 2)
195 sage: A = random_symmetric_psd(V)
196 sage: A.matrix_space()
197 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
198 sage: A.is_positive_semidefinite()
201 A matrix with the desired rank is returned::
203 sage: V = VectorSpace(QQ, 5)
204 sage: A = random_symmetric_psd(V,False,1)
207 sage: A = random_symmetric_psd(V,False,2)
210 sage: A = random_symmetric_psd(V,False,3)
213 sage: A = random_symmetric_psd(V,False,4)
216 sage: A = random_symmetric_psd(V,False,5)
220 If the user asks for a rank that's too high, we fail::
222 sage: V = VectorSpace(QQ, 2)
223 sage: A = random_symmetric_psd(V,False,3)
224 Traceback (most recent call last):
226 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
230 # We construct the matrix from its spectral projectors. Since
231 # there can be at most ``n`` of them, where ``n`` is the dimension
232 # of our vector space, we want to choose a random integer between
233 # ``0`` and ``n`` and then construct that many random elements of
239 # Choose one randomly
240 rank_A
= ZZ
.random_element(n
+1)
241 elif (rank
< 0) or (rank
> n
):
242 # The rank of ``A`` can be at most ``n``.
243 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
244 raise ValueError(msg
)
246 # Use the one the user gave us.
249 if n
== 0 and not accept_zero
:
250 # We're gonna loop forever trying to satisfy this...
251 raise ValueError('You must have accept_zero=True when V is trivial')
253 # Loop until we find a suitable "A" that will then be returned.
255 # Begin with the zero matrix, and add projectors to it if we
257 A
= matrix
.zero(V
.base_ring(), n
, n
)
259 # Careful, begin at idx=1 so that we only generate a projector
260 # when rank_A is greater than zero.
261 while A
.rank() < rank_A
:
262 v
= V
.random_element()
263 A
+= v
.column()*v
.row()
265 if accept_zero
or not A
.is_zero():
266 # We either don't care what ``A`` is, or it's non-zero, so