mjo/ldlt.py: add examples and tests for block_ldlt().

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 729e01152033772217910b1fc46adf50f5d54b8c..461dda3062f1de5e1f0d2acac07d9589d28f8ca5 100644 (file)
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -353,14 +353,144 @@ def block_ldlt(A):
     Perform a block-`LDL^{T}` factorization of the Hermitian
     matrix `A`.
 
+    The standard `LDL^{T}` factorization of a positive-definite matrix
+    `A` factors it as `A = LDL^{T}` where `L` is unit-lower-triangular
+    and `D` is diagonal. If one allows row/column swaps via a
+    permutation matrix `P`, then this factorization can be extended to
+    some positive-semidefinite matrices `A` via the factorization
+    `P^{T}AP = LDL^{T}` that places the zeros at the bottom of `D` to
+    avoid division by zero. These factorizations extend easily to
+    complex Hermitian matrices when one replaces the transpose by the
+    conjugate-transpose.
+
+    However, we can go one step further. If, in addition, we allow `D`
+    to potentially contain `2 \times 2` blocks on its diagonal, then
+    every real or complex Hermitian matrix `A` can be factored as `A =
+    PLDL^{*}P^{T}`. When the row/column swaps are made intelligently,
+    this process is numerically stable over inexact rings like ``RDF``.
+    Bunch and Kaufman describe such a "pivot" scheme that is suitable
+    for the solution of Hermitian systems, and that is how we choose
+    our row and column swaps.
+
     OUTPUT:
 
-    A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
+    If the input matrix is Hermitian, we return a triple `(P,L,D)`
+    such that `A = PLDL^{*}P^{T}` and
 
-      * `P` is a permutation matrix
-      * `L` is unit lower-triangular
+      * `P` is a permutation matrix,
+      * `L` is unit lower-triangular,
       * `D` is a block-diagonal matrix whose blocks are of size
         one or two.
+
+    If the input matrix is not Hermitian, the output from this function
+    is undefined.
+
+    EXAMPLES:
+
+    This three-by-three real symmetric matrix has one positive, one
+    negative, and one zero eigenvalue -- so it is not any flavor of
+    (semi)definite, yet we can still factor it::
+
+        sage: A =  matrix(QQ, [[0, 1, 0],
+        ....:                  [1, 1, 2],
+        ....:                  [0, 2, 0]])
+        sage: P,L,D = block_ldlt(A)
+        sage: P
+        [0 0 1]
+        [1 0 0]
+        [0 1 0]
+        sage: L
+        [  1   0   0]
+        [  2   1   0]
+        [  1 1/2   1]
+        sage: D
+        [ 1| 0| 0]
+        [--+--+--]
+        [ 0|-4| 0]
+        [--+--+--]
+        [ 0| 0| 0]
+        sage: P.T*A*P == L*D*L.T
+        True
+
+    This two-by-two matrix has no standard factorization, but it
+    constitutes its own block-factorization::
+
+        sage: A = matrix(QQ, [ [0,1],
+        ....:                  [1,0] ])
+        sage: block_ldlt(A)
+        (
+        [1 0]  [1 0]  [0 1]
+        [0 1], [0 1], [1 0]
+        )
+
+    The same is true of the following complex Hermitian matrix::
+
+        sage: A = matrix(QQbar, [ [ 0,I],
+        ....:                     [-I,0] ])
+        sage: block_ldlt(A)
+        (
+        [1 0]  [1 0]  [ 0  I]
+        [0 1], [0 1], [-I  0]
+        )
+
+    TESTS:
+
+    All three factors should be the identity when the original matrix is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(6)
+        sage: I = matrix.identity(QQ,n)
+        sage: P,L,D = block_ldlt(I)
+        sage: P == I and L == I and D == I
+        True
+
+    Ensure that a "random" real symmetric matrix is factored correctly::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(6)
+        sage: F = NumberField(x^2 +1, 'I')
+        sage: A = matrix.random(F, n)
+        sage: A = A + A.transpose()
+        sage: P,L,D = block_ldlt(A)
+        sage: A == P*L*D*L.transpose()*P.transpose()
+        True
+
+    Ensure that a "random" complex Hermitian matrix is factored correctly::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(6)
+        sage: F = NumberField(x^2 +1, 'I')
+        sage: A = matrix.random(F, n)
+        sage: A = A + A.conjugate_transpose()
+        sage: P,L,D = block_ldlt(A)
+        sage: A == P*L*D*L.transpose()*P.transpose()
+        True
+
+    Ensure that a "random" complex positive-semidefinite matrix is
+    factored correctly and that the resulting block-diagonal matrix is
+    in fact diagonal::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(6)
+        sage: F = NumberField(x^2 +1, 'I')
+        sage: A = matrix.random(F, n)
+        sage: A = A*A.conjugate_transpose()
+        sage: P,L,D = block_ldlt(A)
+        sage: A == P*L*D*L.transpose()*P.transpose()
+        True
+        sage: diagonal_matrix(D.diagonal()) == D
+        True
+
+    The factorization should be a no-op on diagonal matrices::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(6)
+        sage: A = matrix.diagonal(random_vector(QQ, n))
+        sage: I = matrix.identity(QQ,n)
+        sage: P,L,D = block_ldlt(A)
+        sage: P == I and L == I and A == D
+        True
+
     """
 
     # We have to make at least one copy of the input matrix so that we