mjo/ldlt.py: start a naive Bunch-Kaufman block LDLT.

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 696d78d2053793164a49cd1ac7f25b452b1e2d2c..5c1e65df20c35c44127463101b9753349f6a9ae9 100644 (file)
--- 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