]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/ldlt.py
mjo/ldlt.py: various efficiency improvements to ldlt_fast().
[sage.d.git] / mjo / ldlt.py
1 from sage.all import *
2
3 def is_positive_semidefinite_naive(A):
4 r"""
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.
8
9 SETUP::
10
11 sage: from mjo.ldlt import is_positive_semidefinite_naive
12
13 TESTS:
14
15 The trivial matrix is vaciously positive-semidefinite::
16
17 sage: A = matrix(QQ, 0)
18 sage: A
19 []
20 sage: is_positive_semidefinite_naive(A)
21 True
22
23 """
24 if A.nrows() == 0:
25 return True # vacuously
26 return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() )
27
28
29 def ldlt_naive(A):
30 r"""
31 Perform a pivoted `LDL^{T}` factorization of the Hermitian
32 positive-semidefinite matrix `A`.
33
34 This is a naive, recursive implementation that is inefficient due
35 to Python's lack of tail-call optimization. The pivot strategy is
36 to choose the largest diagonal entry of the matrix at each step,
37 and to permute it into the top-left position. Ultimately this
38 results in a factorization `A = PLDL^{T}P^{T}`, where `P` is a
39 permutation matrix, `L` is unit-lower-triangular, and `D` is
40 diagonal decreasing from top-left to bottom-right.
41
42 ALGORITHM:
43
44 The algorithm is based on the discussion in Golub and Van Loan, but with
45 some "typos" fixed.
46
47 OUTPUT:
48
49 A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
50
51 * `P` is a permutaiton matrix
52 * `L` is unit lower-triangular
53 * `D` is a diagonal matrix whose entries are decreasing from top-left
54 to bottom-right
55
56 SETUP::
57
58 sage: from mjo.ldlt import ldlt_naive, is_positive_semidefinite_naive
59
60 EXAMPLES:
61
62 All three factors should be the identity when the original matrix is::
63
64 sage: I = matrix.identity(QQ,4)
65 sage: P,L,D = ldlt_naive(I)
66 sage: P == I and L == I and D == I
67 True
68
69 TESTS:
70
71 Ensure that a "random" positive-semidefinite matrix is factored correctly::
72
73 sage: set_random_seed()
74 sage: n = ZZ.random_element(5)
75 sage: A = matrix.random(QQ, n)
76 sage: A = A*A.transpose()
77 sage: is_positive_semidefinite_naive(A)
78 True
79 sage: P,L,D = ldlt_naive(A)
80 sage: A == P*L*D*L.transpose()*P.transpose()
81 True
82
83 """
84 n = A.nrows()
85
86 # Use the fraction field of the given matrix so that division will work
87 # when (for example) our matrix consists of integer entries.
88 ring = A.base_ring().fraction_field()
89
90 if n == 0 or n == 1:
91 # We can get n == 0 if someone feeds us a trivial matrix.
92 P = matrix.identity(ring, n)
93 L = matrix.identity(ring, n)
94 D = A
95 return (P,L,D)
96
97 A1 = A.change_ring(ring)
98 diags = A1.diagonal()
99 s = diags.index(max(diags))
100 P1 = copy(A1.matrix_space().identity_matrix())
101 P1.swap_rows(0,s)
102 A1 = P1.T * A1 * P1
103 alpha1 = A1[0,0]
104
105 # Golub and Van Loan mention in passing what to do here. This is
106 # only sensible if the matrix is positive-semidefinite, because we
107 # are assuming that we can set everything else to zero as soon as
108 # we hit the first on-diagonal zero.
109 if alpha1 == 0:
110 P = A1.matrix_space().identity_matrix()
111 L = P
112 D = A1.matrix_space().zero()
113 return (P,L,D)
114
115 v1 = A1[1:n,0]
116 A2 = A1[1:,1:]
117
118 P2, L2, D2 = ldlt_naive(A2 - (v1*v1.transpose())/alpha1)
119
120 P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
121 [0*v1, P2]])
122 L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
123 [P2.transpose()*v1/alpha1, L2]])
124 D1 = block_matrix(2,2, [[alpha1, ZZ(0)],
125 [0*v1, D2]])
126
127 return (P1,L1,D1)
128
129
130
131 def ldlt_fast(A):
132 r"""
133 Perform a fast, pivoted `LDL^{T}` factorization of the Hermitian
134 positive-semidefinite matrix `A`.
135
136 This function is much faster than ``ldlt_naive`` because the
137 tail-recursion has been unrolled into a loop.
138 """
139 ring = A.base_ring().fraction_field()
140 A = A.change_ring(ring)
141
142 # Keep track of the permutations in a vector rather than in a
143 # matrix, for efficiency.
144 n = A.nrows()
145 p = list(range(n))
146
147 for k in range(n):
148 # We need to loop once for every diagonal entry in the
149 # matrix. So, as many times as it has rows/columns. At each
150 # step, we obtain the permutation needed to put things in the
151 # right place, then the "next" entry (alpha) of D, and finally
152 # another column of L.
153 diags = A.diagonal()[k:n]
154 alpha = max(diags)
155
156 # We're working *within* the matrix ``A``, so every index is
157 # offset by k. For example: after the second step, we should
158 # only be looking at the lower 3-by-3 block of a 5-by-5 matrix.
159 s = k + diags.index(alpha)
160
161 # Move the largest diagonal element up into the top-left corner
162 # of the block we're working on (the one starting from index k,k).
163 # Presumably this is faster than hitting the thing with a
164 # permutation matrix.
165 #
166 # Since "L" is stored in the lower-left "half" of "A", it's a
167 # good thing that we need to permuts "L," too. This is due to
168 # how P2.T appears in the recursive algorithm applied to the
169 # "current" column of L There, P2.T is computed recusively, as
170 # 1 x P3.T, and P3.T = 1 x P4.T, etc, from the bottom up. All
171 # are eventually applied to "v" in order. Here we're working
172 # from the top down, and rather than keep track of what
173 # permutations we need to perform, we just perform them as we
174 # go along. No recursion needed.
175 A.swap_columns(k,s)
176 A.swap_rows(k,s)
177
178 # Update the permutation "matrix" with the swap we just did.
179 p_k = p[k]
180 p[k] = p[s]
181 p[s] = p_k
182
183 # Now the largest diagonal is in the top-left corner of the
184 # block below and to the right of index k,k. When alpha is
185 # zero, we can just leave the rest of the D/L entries
186 # zero... which is exactly how they start out.
187 if alpha != 0:
188 # Update the "next" block of A that we'll work on during
189 # the following iteration. I think it's faster to get the
190 # entries of a row than a column here?
191 for i in range(n-k-1):
192 for j in range(i+1):
193 A[k+1+j,k+1+i] = A[k+1+j,k+1+i] - A[k,k+1+j]*A[k,k+1+i]/alpha
194 A[k+1+i,k+1+j] = A[k+1+j,k+1+i] # keep it symmetric!
195
196 for i in range(n-k-1):
197 # Store the "new" (kth) column of L, being sure to set
198 # the lower-left "half" from the upper-right "half"
199 A[k+i+1,k] = A[k,k+1+i]/alpha
200
201 MS = A.matrix_space()
202 P = MS.matrix(lambda i,j: p[j] == i)
203 D = MS.diagonal_matrix(A.diagonal())
204
205 for i in range(n):
206 A[i,i] = 1
207 for j in range(i+1,n):
208 A[i,j] = 0
209
210 return P,A,D