From fe85934cbfb26ff0dc7e5ef81411fe73eed739ff Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 5 Oct 2020 15:23:00 -0400 Subject: [PATCH] mjo/ldlt.py: add fast block-LDLT-based is_positive_semidefinite(). --- mjo/ldlt.py | 47 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/mjo/ldlt.py b/mjo/ldlt.py index a86dbbe..1c2663e 100644 --- a/mjo/ldlt.py +++ b/mjo/ldlt.py @@ -24,3 +24,50 @@ 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 is_positive_semidefinite(A): + r""" + A fast positive-semidefinite check based on the block-LDLT + factorization. + + SETUP:: + + sage: from mjo.ldlt import (is_positive_semidefinite, + ....: is_positive_semidefinite_naive) + + TESTS: + + Check that the naive and fast answers are the same, in general:: + + sage: set_random_seed() + sage: F = NumberField(x^2 + 1, 'I') + sage: from sage.misc.prandom import choice + sage: ring = choice([ZZ,QQ,F]) + sage: A = matrix.random(ring, 10) + sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A) + True + + Check that the naive and fast answers are the same for a Hermitian + matrix:: + + sage: set_random_seed() + sage: F = NumberField(x^2 + 1, 'I') + sage: from sage.misc.prandom import choice + sage: ring = choice([ZZ,QQ,F]) + sage: A = matrix.random(ring, 10); A = A + A.conjugate_transpose() + sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A) + True + + """ + if not A.is_hermitian(): + return False + _,_,d = A._block_ldlt() + for d_i in d: + if d_i.nrows() == 1: + if d_i < 0: + return False + else: + # A 2x2 block indicates that it's indefinite + return False + return True -- 2.44.2