]>
gitweb.michael.orlitzky.com - sage.d.git/blob - symmetric_psd.py
fae0ae2307bb4e8b0583cf5b2200951f66a7e256
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 unit_eigenvectors(A
):
20 Return the unit eigenvectors of a symmetric positive-definite matrix.
24 - ``A`` - The matrix whose eigenvectors we want to compute.
28 A list of (eigenvalue, eigenvector) pairs where each eigenvector is
29 associated with its paired eigenvalue of ``A`` and has norm `1`.
33 sage: A = matrix(QQ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
34 sage: unit_evs = unit_eigenvectors(A)
35 sage: bool(unit_evs[0][1].norm() == 1)
37 sage: bool(unit_evs[1][1].norm() == 1)
39 sage: bool(unit_evs[2][1].norm() == 1)
43 # This will give us a list of lists whose elements are the
44 # eigenvectors we want.
45 ev_lists
= [ (val
,vecs
) for (val
,vecs
,multiplicity
)
46 in A
.eigenvectors_right() ]
48 # Pair each eigenvector with its eigenvalue and normalize it.
49 evs
= [ [(l
, vec
/vec
.norm()) for vec
in vecs
] for (l
,vecs
) in ev_lists
]
51 # Flatten the list, abusing the fact that "+" is overloaded on lists.
59 Factor a symmetric positive-semidefinite matrix ``A`` into
64 - ``A`` - The matrix to factor.
68 A matrix ``X`` such that `A = XX^{T}`.
72 Since ``A`` is symmetric and positive-semidefinite, we can
73 diagonalize it by some matrix `$Q$` whose columns are orthogonal
74 eigenvectors of ``A``. Then,
78 From this representation we can take the square root of `$D$`
79 (since all eigenvalues of ``A`` are nonnegative). If we then let
80 `$X = Q*sqrt(D)*Q^{T}$`, we have,
82 `$XX^{T} = Q*sqrt(D)*Q^{T}Q*sqrt(D)*Q^{T} = Q*D*Q^{T} = A$`
86 In principle, this is the algorithm used, although we ignore the
87 eigenvectors corresponding to the eigenvalue zero. Thus if `$rank(A)
88 = k$`, the matrix `$Q$` will have dimention `$n \times k$`, and
89 `$D$` will have dimension `$k \times k$`. In the end everything
94 sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
95 sage: X = factor_psd(A)
96 sage: A2 = matrix_simplify_full(X*X.transpose())
102 # Get the eigenvectors, and filter out the ones that correspond to
103 # the eigenvalue zero.
104 all_evs
= unit_eigenvectors(A
)
105 evs
= [ (val
,vec
) for (val
,vec
) in all_evs
if not val
== 0 ]
107 d
= [ sqrt(val
) for (val
,vec
) in evs
]
108 root_D
= diagonal_matrix(d
).change_ring(A
.base_ring())
110 Q
= matrix(A
.base_ring(), [ vec
for (val
,vec
) in evs
]).transpose()
112 return Q
*root_D
*Q
.transpose()