From e060818e0ceaf3bc2bf1dedfed4addd2ab26d036 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 2 Oct 2020 07:37:35 -0400 Subject: [PATCH] mjo/ldlt.py: start a naive Bunch-Kaufman block LDLT. --- mjo/ldlt.py | 129 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 129 insertions(+) diff --git a/mjo/ldlt.py b/mjo/ldlt.py index 696d78d..5c1e65d 100644 --- a/mjo/ldlt.py +++ b/mjo/ldlt.py @@ -208,3 +208,132 @@ def ldlt_fast(A): 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 -- 2.44.2