]>
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: set_random_seed()
195 sage: V = VectorSpace(QQ, 2)
196 sage: A = random_symmetric_psd(V)
197 sage: A.matrix_space()
198 Full MatrixSpace of 2 by 2 dense matrices over Rational Field
199 sage: A.is_positive_semidefinite()
202 A matrix with the desired rank is returned::
204 sage: set_random_seed()
205 sage: V = VectorSpace(QQ, 5)
206 sage: A = random_symmetric_psd(V,False,1)
209 sage: A = random_symmetric_psd(V,False,2)
212 sage: A = random_symmetric_psd(V,False,3)
215 sage: A = random_symmetric_psd(V,False,4)
218 sage: A = random_symmetric_psd(V,False,5)
222 If the user asks for a rank that's too high, we fail::
224 sage: set_random_seed()
225 sage: V = VectorSpace(QQ, 2)
226 sage: A = random_symmetric_psd(V,False,3)
227 Traceback (most recent call last):
229 ValueError: The ``rank`` must be between 0 and the dimension of ``V``.
233 # We construct the matrix from its spectral projectors. Since
234 # there can be at most ``n`` of them, where ``n`` is the dimension
235 # of our vector space, we want to choose a random integer between
236 # ``0`` and ``n`` and then construct that many random elements of
242 # Choose one randomly
243 rank_A
= ZZ
.random_element(n
+1)
244 elif (rank
< 0) or (rank
> n
):
245 # The rank of ``A`` can be at most ``n``.
246 msg
= 'The ``rank`` must be between 0 and the dimension of ``V``.'
247 raise ValueError(msg
)
249 # Use the one the user gave us.
252 if n
== 0 and not accept_zero
:
253 # We're gonna loop forever trying to satisfy this...
254 raise ValueError('You must have accept_zero=True when V is trivial')
256 # Loop until we find a suitable "A" that will then be returned.
258 # Begin with the zero matrix, and add projectors to it if we
260 A
= matrix
.zero(V
.base_ring(), n
, n
)
262 # Careful, begin at idx=1 so that we only generate a projector
263 # when rank_A is greater than zero.
264 while A
.rank() < rank_A
:
265 v
= V
.random_element()
266 A
+= v
.column()*v
.row()
268 if accept_zero
or not A
.is_zero():
269 # We either don't care what ``A`` is, or it's non-zero, so