From fdafcf90579fd42c869c669f0575dc5812485bc6 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky 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