mjo/ldlt.py: move implementation to SageMath proper.

[sage.d.git] / mjo / ldlt.py

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 8e5f6930da7630d11704458ca14613d9aab64718..a86dbbefe5f396a197a76af33dc60583cb7cb971 100644 (file)
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -24,102 +24,3 @@ def is_positive_semidefinite_naive(A):
     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)