From fdafcf90579fd42c869c669f0575dc5812485bc6 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 27 Sep 2020 21:40:13 -0400
Subject: [PATCH] mjo/ldlt.py: add naive, pivoted LDLT matrix factorization.

---
 mjo/all.py  |   1 +
 mjo/ldlt.py | 125 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 126 insertions(+)
 create mode 100644 mjo/ldlt.py

diff --git a/mjo/all.py b/mjo/all.py
index de2212d..3854ac3 100644
--- a/mjo/all.py
+++ b/mjo/all.py
@@ -6,6 +6,7 @@ in his script. Instead, he can just `from mjo.all import *`.
 from mjo.basis_repr import *
 from mjo.cone.all import *
 from mjo.eja.all import *
+from mjo.ldlt import *
 from mjo.interpolation import *
 from mjo.misc import *
 from mjo.orthogonal_polynomials import *
diff --git a/mjo/ldlt.py b/mjo/ldlt.py
new file mode 100644
index 0000000..8e5f693
--- /dev/null
+++ b/mjo/ldlt.py
@@ -0,0 +1,125 @@
+from sage.all import *
+
+def is_positive_semidefinite_naive(A):
+    r"""
+    A naive positive-semidefinite check that tests the eigenvalues for
+    nonnegativity.  We follow the sage convention that positive
+    (semi)definite matrices must be symmetric or Hermitian.
+
+    SETUP::
+
+        sage: from mjo.ldlt import is_positive_semidefinite_naive
+
+    TESTS:
+
+    The trivial matrix is vaciously positive-semidefinite::
+
+        sage: A = matrix(QQ, 0)
+        sage: A
+        []
+        sage: is_positive_semidefinite_naive(A)
+        True
+
+    """
+    if A.nrows() == 0:
+        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
+    positive-semidefinite matrix `A`.
+
+    This is a naive, recursive implementation that is inefficient due
+    to Python's lack of tail-call optimization. The pivot strategy is
+    to choose the largest diagonal entry of the matrix at each step,
+    and to permute it into the top-left position. Ultimately this
+    results in a factorization `A = PLDL^{T}P^{T}`, where `P` is a
+    permutation matrix, `L` is unit-lower-triangular, and `D` is
+    diagonal decreasing from top-left to bottom-right.
+
+    ALGORITHM:
+
+    The algorithm is based on the discussion in Golub and Van Loan, but with
+    some "typos" fixed.
+
+    OUTPUT:
+
+    A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
+
+      * `P` is a permutaiton matrix
+      * `L` is unit lower-triangular
+      * `D` is a diagonal matrix whose entries are decreasing from top-left
+        to bottom-right
+
+    SETUP::
+
+        sage: from mjo.ldlt import ldlt_naive, is_positive_semidefinite_naive
+
+    EXAMPLES:
+
+    All three factors should be the identity when the original matrix is::
+
+        sage: I = matrix.identity(QQ,4)
+        sage: P,L,D = ldlt_naive(I)
+        sage: P == I and L == I and D == I
+        True
+
+    TESTS:
+
+    Ensure that a "random" positive-semidefinite matrix is factored correctly::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: A = matrix.random(QQ, n)
+        sage: A = A*A.transpose()
+        sage: is_positive_semidefinite_naive(A)
+        True
+        sage: P,L,D = ldlt_naive(A)
+        sage: A == P*L*D*L.transpose()*P.transpose()
+        True
+
+    """
+    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)
+
+    A1 = A.change_ring(ring)
+    diags = A1.diagonal()
+    s = diags.index(max(diags))
+    P1 = copy(A1.matrix_space().identity_matrix())
+    A1 = P1.T * A1 * P1
+    alpha1 = A1[0,0]
+
+    # Golub and Van Loan mention in passing what to do here. This is
+    # only sensible if the matrix is positive-semidefinite, because we
+    # are assuming that we can set everything else to zero as soon as
+    # we hit the first on-diagonal zero.
+    if alpha1 == 0:
+        P = A1.matrix_space().identity_matrix()
+        L = P
+        D = A1.matrix_space().zero()
+        return (P,L,D)
+
+    v1 = A1[1:n,0]
+    A2 = A1[1:,1:]
+
+    P2, L2, D2 = ldlt_naive(A2 - (v1*v1.transpose())/alpha1)
+
+    P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
+                               [0*v1,  P2]])
+    L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
+                            [P2.transpose()*v1/alpha1, L2]])
+    D1 = block_matrix(2,2, [[alpha1, ZZ(0)],
+                            [0*v1,   D2]])
+
+    return (P1,L1,D1)
-- 
2.44.2