From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 5 Oct 2020 19:23:00 +0000 (-0400)
Subject: mjo/ldlt.py: add fast block-LDLT-based is_positive_semidefinite().
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=fe85934cbfb26ff0dc7e5ef81411fe73eed739ff;p=sage.d.git

mjo/ldlt.py: add fast block-LDLT-based is_positive_semidefinite().
---

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