]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/ldlt.py
8e5f6930da7630d11704458ca14613d9aab64718
3 def is_positive_semidefinite_naive(A
):
5 A naive positive-semidefinite check that tests the eigenvalues for
6 nonnegativity. We follow the sage convention that positive
7 (semi)definite matrices must be symmetric or Hermitian.
11 sage: from mjo.ldlt import is_positive_semidefinite_naive
15 The trivial matrix is vaciously positive-semidefinite::
17 sage: A = matrix(QQ, 0)
20 sage: is_positive_semidefinite_naive(A)
25 return True # vacuously
26 return A
.is_hermitian() and all( v
>= 0 for v
in A
.eigenvalues() )
30 Perform a pivoted `LDL^{T}` factorization of the Hermitian
31 positive-semidefinite matrix `A`.
33 This is a naive, recursive implementation that is inefficient due
34 to Python's lack of tail-call optimization. The pivot strategy is
35 to choose the largest diagonal entry of the matrix at each step,
36 and to permute it into the top-left position. Ultimately this
37 results in a factorization `A = PLDL^{T}P^{T}`, where `P` is a
38 permutation matrix, `L` is unit-lower-triangular, and `D` is
39 diagonal decreasing from top-left to bottom-right.
43 The algorithm is based on the discussion in Golub and Van Loan, but with
48 A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
50 * `P` is a permutaiton matrix
51 * `L` is unit lower-triangular
52 * `D` is a diagonal matrix whose entries are decreasing from top-left
57 sage: from mjo.ldlt import ldlt_naive, is_positive_semidefinite_naive
61 All three factors should be the identity when the original matrix is::
63 sage: I = matrix.identity(QQ,4)
64 sage: P,L,D = ldlt_naive(I)
65 sage: P == I and L == I and D == I
70 Ensure that a "random" positive-semidefinite matrix is factored correctly::
72 sage: set_random_seed()
73 sage: n = ZZ.random_element(5)
74 sage: A = matrix.random(QQ, n)
75 sage: A = A*A.transpose()
76 sage: is_positive_semidefinite_naive(A)
78 sage: P,L,D = ldlt_naive(A)
79 sage: A == P*L*D*L.transpose()*P.transpose()
85 # Use the fraction field of the given matrix so that division will work
86 # when (for example) our matrix consists of integer entries.
87 ring
= A
.base_ring().fraction_field()
90 # We can get n == 0 if someone feeds us a trivial matrix.
91 P
= matrix
.identity(ring
, n
)
92 L
= matrix
.identity(ring
, n
)
96 A1
= A
.change_ring(ring
)
98 s
= diags
.index(max(diags
))
99 P1
= copy(A1
.matrix_space().identity_matrix())
103 # Golub and Van Loan mention in passing what to do here. This is
104 # only sensible if the matrix is positive-semidefinite, because we
105 # are assuming that we can set everything else to zero as soon as
106 # we hit the first on-diagonal zero.
108 P
= A1
.matrix_space().identity_matrix()
110 D
= A1
.matrix_space().zero()
116 P2
, L2
, D2
= ldlt_naive(A2
- (v1
*v1
.transpose())/alpha1
)
118 P1
= P1
*block_matrix(2,2, [[ZZ(1), ZZ(0)],
120 L1
= block_matrix(2,2, [[ZZ(1), ZZ(0)],
121 [P2
.transpose()*v1
/alpha1
, L2
]])
122 D1
= block_matrix(2,2, [[alpha1
, ZZ(0)],