doubly_nonnegative.py: Begin implementing the extreme vector algorithms.

[sage.d.git] / mjo / cone / doubly_nonnegative.py

diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py
index b071e41dc332a709a527588bd8f665e80eeb4c17..b43c974b2076a5dbda14ba43f679fd8caee30327 100644 (file)
--- a/mjo/cone/doubly_nonnegative.py
+++ b/mjo/cone/doubly_nonnegative.py
@@ -19,7 +19,7 @@ from sage.all import *
 from os.path import abspath
 from site import addsitedir
 addsitedir(abspath('../../'))
-from mjo.cone.symmetric_psd import factor_psd
+from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
 
 
 
@@ -29,7 +29,7 @@ def is_doubly_nonnegative(A):
 
     INPUT:
 
-      - ``A`` - The matrix in question
+    - ``A`` - The matrix in question
 
     OUTPUT:
 
@@ -54,28 +54,94 @@ def is_doubly_nonnegative(A):
     """
 
     if A.base_ring() == SR:
-        msg = 'The base ring of ``A`` cannot be the Symbolic Ring'
+        msg = 'The matrix ``A`` cannot be the symbolic.'
         raise ValueError.new(msg)
 
-    # First make sure that ``A`` is symmetric.
-    if not A.is_symmetric():
-        return False
-
     # Check that all of the entries of ``A`` are nonnegative.
     if not all([ a >= 0 for a in A.list() ]):
         return False
 
-    # If ``A`` is symmetric and non-negative, we only need to check
-    # that it is positive semidefinite. For that we can consult its
-    # minimum eigenvalue, which should be zero or greater. Since ``A``
-    # is symmetric, its eigenvalues are guaranteed to be real.
-    return min(A.eigenvalues()) >= 0
+    # It's nonnegative, so all we need to do is check that it's
+    # symmetric positive-semidefinite.
+    return is_symmetric_psd(A)
+
+
 
+def has_admissible_extreme_rank(A):
+    """
+    The extreme matrices of the doubly-nonnegative cone have some
+    restrictions on their ranks. This function checks to see whether or
+    not ``A`` could be extreme based on its rank.
+
+    INPUT:
+
+    - ``A`` - The matrix in question
+
+    OUTPUT:
+
+    ``False`` if the rank of ``A`` precludes it from being an extreme
+    matrix of the doubly-nonnegative cone, ``True`` otherwise.
+
+    REFERENCE:
+
+    Hamilton-Jester, Crista Lee; Li, Chi-Kwong. Extreme Vectors of
+    Doubly Nonnegative Matrices. Rocky Mountain Journal of Mathematics
+    26 (1996), no. 4, 1371--1383. doi:10.1216/rmjm/1181071993.
+    http://projecteuclid.org/euclid.rmjm/1181071993.
+
+    EXAMPLES:
+
+    The zero matrix has rank zero, which is admissible::
+
+    sage: A = zero_matrix(QQ, 5, 5)
+    sage: has_admissible_extreme_rank(A)
+    True
+
+    """
+    if not A.is_symmetric():
+        raise ValueError('The matrix ``A`` must be symmetric.')
+
+    r = rank(A)
+    n = A.nrows() # Columns would work, too, since ``A`` is symmetric.
+
+    if r == 0:
+        # Zero is in the doubly-nonnegative cone.
+        return True
+
+    # See Theorem 3.1 in the cited reference.
+    if r == 2:
+        return False
+
+    if n.mod(2) == 0:
+        # n is even
+        return r <= max(1, n-3)
+    else:
+        # n is odd
+        return r <= max(1, n-2)
 
 
 def is_extreme_doubly_nonnegative(A):
     """
     Returns ``True`` if the given matrix is an extreme matrix of the
     doubly-nonnegative cone, and ``False`` otherwise.
+
+    EXAMPLES:
+
+    The zero matrix is an extreme matrix::
+
+        sage: A = zero_matrix(QQ, 5, 5)
+        sage: is_extreme_doubly_nonnegative(A)
+        True
+
     """
+
+    r = A.rank()
+
+    if r == 0:
+        # Short circuit, we know the zero matrix is extreme.
+        return True
+
+    if not is_admissible_extreme_rank(r):
+        return False
+
     raise NotImplementedError()