X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fldlt.py;h=5c1e65df20c35c44127463101b9753349f6a9ae9;hb=e060818e0ceaf3bc2bf1dedfed4addd2ab26d036;hp=8e5f6930da7630d11704458ca14613d9aab64718;hpb=fdafcf90579fd42c869c669f0575dc5812485bc6;p=sage.d.git diff --git a/mjo/ldlt.py b/mjo/ldlt.py index 8e5f693..5c1e65d 100644 --- a/mjo/ldlt.py +++ b/mjo/ldlt.py @@ -25,6 +25,7 @@ def is_positive_semidefinite_naive(A): return True # vacuously return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() ) + def ldlt_naive(A): r""" Perform a pivoted `LDL^{T}` factorization of the Hermitian @@ -97,6 +98,7 @@ def ldlt_naive(A): diags = A1.diagonal() s = diags.index(max(diags)) P1 = copy(A1.matrix_space().identity_matrix()) + P1.swap_rows(0,s) A1 = P1.T * A1 * P1 alpha1 = A1[0,0] @@ -123,3 +125,215 @@ def ldlt_naive(A): [0*v1, D2]]) return (P1,L1,D1) + + + +def ldlt_fast(A): + r""" + Perform a fast, pivoted `LDL^{T}` factorization of the Hermitian + positive-semidefinite matrix `A`. + + This function is much faster than ``ldlt_naive`` because the + tail-recursion has been unrolled into a loop. + """ + ring = A.base_ring().fraction_field() + A = A.change_ring(ring) + + # Keep track of the permutations in a vector rather than in a + # matrix, for efficiency. + n = A.nrows() + p = list(range(n)) + + for k in range(n): + # We need to loop once for every diagonal entry in the + # matrix. So, as many times as it has rows/columns. At each + # step, we obtain the permutation needed to put things in the + # right place, then the "next" entry (alpha) of D, and finally + # another column of L. + diags = A.diagonal()[k:n] + alpha = max(diags) + + # We're working *within* the matrix ``A``, so every index is + # offset by k. For example: after the second step, we should + # only be looking at the lower 3-by-3 block of a 5-by-5 matrix. + s = k + diags.index(alpha) + + # Move the largest diagonal element up into the top-left corner + # of the block we're working on (the one starting from index k,k). + # Presumably this is faster than hitting the thing with a + # permutation matrix. + # + # Since "L" is stored in the lower-left "half" of "A", it's a + # good thing that we need to permuts "L," too. This is due to + # how P2.T appears in the recursive algorithm applied to the + # "current" column of L There, P2.T is computed recusively, as + # 1 x P3.T, and P3.T = 1 x P4.T, etc, from the bottom up. All + # are eventually applied to "v" in order. Here we're working + # from the top down, and rather than keep track of what + # permutations we need to perform, we just perform them as we + # go along. No recursion needed. + A.swap_columns(k,s) + A.swap_rows(k,s) + + # Update the permutation "matrix" with the swap we just did. + p_k = p[k] + p[k] = p[s] + p[s] = p_k + + # Now the largest diagonal is in the top-left corner of the + # block below and to the right of index k,k. When alpha is + # zero, we can just leave the rest of the D/L entries + # zero... which is exactly how they start out. + if alpha != 0: + # Update the "next" block of A that we'll work on during + # the following iteration. I think it's faster to get the + # entries of a row than a column here? + for i in range(n-k-1): + for j in range(i+1): + 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 + A[k+1+i,k+1+j] = A[k+1+j,k+1+i] # keep it symmetric! + + for i in range(n-k-1): + # Store the "new" (kth) column of L, being sure to set + # the lower-left "half" from the upper-right "half" + A[k+i+1,k] = A[k,k+1+i]/alpha + + MS = A.matrix_space() + P = MS.matrix(lambda i,j: p[j] == i) + D = MS.diagonal_matrix(A.diagonal()) + + for i in range(n): + A[i,i] = 1 + for j in range(i+1,n): + A[i,j] = 0 + + return P,A,D + + +def block_ldlt_naive(A, check_hermitian=False): + r""" + Perform a block-`LDL^{T}` factorization of the Hermitian + matrix `A`. + + This is a naive, recursive implementation akin to + ``ldlt_naive()``, where the pivots (and resulting diagonals) are + either `1 \times 1` or `2 \times 2` blocks. The pivots are chosen + using the Bunch-Kaufmann scheme that is both fast and numerically + stable. + + OUTPUT: + + A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where, + + * `P` is a permutation matrix + * `L` is unit lower-triangular + * `D` is a block-diagonal matrix whose entries are decreasing + from top-left to bottom-right and whose blocks are of size + one or two. + """ + n = A.nrows() + + # Use the fraction field of the given matrix so that division will work + # when (for example) our matrix consists of integer entries. + ring = A.base_ring().fraction_field() + + if n == 0 or n == 1: + # We can get n == 0 if someone feeds us a trivial matrix. + P = matrix.identity(ring, n) + L = matrix.identity(ring, n) + D = A + return (P,L,D) + + alpha = (1 + ZZ(17).sqrt()) * ~ZZ(8) + A1 = A.change_ring(ring) + + # Bunch-Kaufmann step 1, Higham step "zero." We use Higham's + # "omega" notation instead of Bunch-Kaufman's "lamda" because + # lambda means other things in the same context. + column_1_subdiag = A1[1:,0].list() + omega_1 = max([ a_i1.abs() for a_i1 in column_1_subdiag ]) + + if omega_1 == 0: + # "There's nothing to do at this step of the algorithm," + # which means that our matrix looks like, + # + # [ 1 0 ] + # [ 0 B ] + # + B = A1[1:,1:] + P2, L2, D2 = ldlt_naive(B) + P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [ZZ(0), P2]]) + L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [ZZ(0), L2]]) + D1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [ZZ(0), D2]]) + return (P1,L1,D1) + + if A1[0,0].abs() > alpha*omega_1: + # Higham step (1) + # The top-left entry is our 1x1 pivot. + C = A1[1:n,0] + B = A1[1:,1:] + + P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0]) + + P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [0*C, P2]]) + L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [P2.transpose()*C/A1[0,0], L2]]) + D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)], + [0*C, D2]]) + + return (P1,L1,D1) + + + r = 1 + column_1_subdiag.index(omega_1) + + # If the matrix is Hermitian, we need only look at the above- + # diagonal entries to find the off-diagonal of maximal magnitude. + omega_r = max( a_rj.abs() for a_rj in A1[:r,r].list() ) + + if A1[0,0].abs()*omega_r >= alpha*(omega_1^2): + # Higham step (2) + # The top-left entry is our 1x1 pivot. + C = A1[1:n,0] + B = A1[1:,1:] + + P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0]) + + P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [0*C, P2]]) + L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [P2.transpose()*C/A1[0,0], L2]]) + D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)], + [0*C, D2]]) + + return (P1,L1,D1) + + + if A1[r,r].abs() > alpha*omega_r: + # Higham step (3) + # Another 1x1 pivot, but this time swapping indices 0,r. + P1 = copy(A1.matrix_space().identity_matrix()) + P1.swap_rows(0,s) + A1 = P1.T * A1 * P1 + + # The top-left entry is now our 1x1 pivot. + C = A1[1:n,0] + B = A1[1:,1:] + + P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0]) + + P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)], + [0*C, P2]]) + L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)], + [P2.transpose()*C/A1[0,0], L2]]) + D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)], + [0*C, D2]]) + + return (P1,L1,D1) + + # Higham step (4) + # If we made it here, we have to do a 2x2 pivot. + return None