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

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index a86dbbefe5f396a197a76af33dc60583cb7cb971..1c2663ea3d574360059b07d22f9af8163ebdc2e7 100644 (file)
--- 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