mjo/ldlt.py

   1 from sage.all import *
   2
   3 def is_positive_semidefinite_naive(A):
   4     r"""
   5     A naive positive-semidefinite check that tests the eigenvalues for
   6     nonnegativity.  We follow the sage convention that positive
   7     (semi)definite matrices must be symmetric or Hermitian.
   8
   9     SETUP::
  10
  11         sage: from mjo.ldlt import is_positive_semidefinite_naive
  12
  13     TESTS:
  14
  15     The trivial matrix is vaciously positive-semidefinite::
  16
  17         sage: A = matrix(QQ, 0)
  18         sage: A
  19         []
  20         sage: is_positive_semidefinite_naive(A)
  21         True
  22
  23     """
  24     if A.nrows() == 0:
  25         return True # vacuously
  26     return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() )
  27
  28
  29 def is_positive_semidefinite(A):
  30     r"""
  31     A fast positive-semidefinite check based on the block-LDLT
  32     factorization.
  33
  34     SETUP::
  35
  36         sage: from mjo.ldlt import (is_positive_semidefinite,
  37         ....:                       is_positive_semidefinite_naive)
  38
  39     TESTS:
  40
  41     Check that the naive and fast answers are the same, in general::
  42
  43         sage: set_random_seed()
  44         sage: F = NumberField(x^2 + 1, 'I')
  45         sage: from sage.misc.prandom import choice
  46         sage: ring = choice([ZZ,QQ,F])
  47         sage: A = matrix.random(ring, 10)
  48         sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A)
  49         True
  50
  51     Check that the naive and fast answers are the same for a Hermitian
  52     matrix::
  53
  54         sage: set_random_seed()
  55         sage: F = NumberField(x^2 + 1, 'I')
  56         sage: from sage.misc.prandom import choice
  57         sage: ring = choice([ZZ,QQ,F])
  58         sage: A = matrix.random(ring, 10); A = A + A.conjugate_transpose()
  59         sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A)
  60         True
  61
  62     """
  63     if not A.is_hermitian():
  64         return False
  65     _,_,d = A._block_ldlt()
  66     for d_i in d:
  67         if d_i.nrows() == 1:
  68             if d_i < 0:
  69                 return False
  70         else:
  71             # A 2x2 block indicates that it's indefinite
  72             return False
  73     return True