From f471b4e3c164cbc6ed79594a30540f4ef07911bc Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 28 Sep 2020 15:27:37 -0400 Subject: [PATCH] mjo/ldlt.py: various efficiency improvements to ldlt_fast(). --- mjo/ldlt.py | 61 ++++++++++++++++++++++++++--------------------------- 1 file changed, 30 insertions(+), 31 deletions(-) diff --git a/mjo/ldlt.py b/mjo/ldlt.py index 4561528..696d78d 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 @@ -135,18 +136,12 @@ def ldlt_fast(A): This function is much faster than ``ldlt_naive`` because the tail-recursion has been unrolled into a loop. """ - n = A.nrows() ring = A.base_ring().fraction_field() - A = A.change_ring(ring) - # Don't try to store the results in the lower-left-hand corner of - # "A" itself; there lies madness. - L = copy(A.matrix_space().identity_matrix()) - D = copy(A.matrix_space().zero()) - # 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): @@ -167,29 +162,27 @@ def ldlt_fast(A): # 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) - # Have to do L, too, to keep track of the "P2.T" (which is 1 x - # P3.T which is 1 x P4 T)... in the recursive - # algorithm. There, we compute P2^T from the bottom up. Here, - # we apply the permutations one at a time, essentially - # building them top-down (but just applying them instead of - # building them. - L.swap_columns(k,s) - L.swap_rows(k,s) - - # Update the permutation "matrix" with the next swap. + # 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.... - # Note: same as ``pivot``. - D[k,k] = alpha - - # When alpha is zero, we can just leave the rest of the D/L entries + # 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 @@ -197,15 +190,21 @@ def ldlt_fast(A): # entries of a row than a column here? for i in range(n-k-1): for j in range(i+1): - A[k+1+i,k+1+j] = A[k+1+i,k+1+j] - A[k,k+1+i]*A[k,k+1+j]/alpha - A[k+1+j,k+1+i] = A[k+1+i,k+1+j] # keep it symmetric! + 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! - # Store the "new" (kth) column of L. for i in range(n-k-1): - # Set the lower-left "half" from the upper-right "half"... - L[k+i+1,k] = A[k,k+1+i]/alpha + # 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()) - I = A.matrix_space().identity_matrix() - P = matrix.column( I.row(p[j]) for j in range(n) ) + for i in range(n): + A[i,i] = 1 + for j in range(i+1,n): + A[i,j] = 0 - return P,L,D + return P,A,D -- 2.44.2