]>
gitweb.michael.orlitzky.com - sage.d.git/blob - 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 raise ValueError('The base ring of ``A`` must be either exact or symbolic.')
183 if not A
.base_ring().is_field():
184 raise ValueError('The base ring of ``A`` must be a field.')
186 if not A
.base_ring() is SR
:
187 # Change the base field of ``A`` so that we are sure we can take
188 # roots. The symbolic ring has no algebraic closure.
189 A
= A
.change_ring(A
.base_ring().algebraic_closure())
192 # Get the eigenvectors, and filter out the ones that correspond to
193 # the eigenvalue zero.
194 all_evs
= unit_eigenvectors(A
)
195 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
197 d
= [ sqrt(val
) for (val
,vec
) in evs
]
198 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
200 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
202 return Q
*root_D
*Q
.transpose()