Implement the is(_extreme)_completely_positive functions.

diff --git a/mjo/cone/completely_positive.py b/mjo/cone/completely_positive.py
index 0f94ca50e9d42d3540a7bb34e84c8ac2b9f57b55..ac5f14435713b4b524c53a2e1c73607ee004b323 100644 (file)
--- a/mjo/cone/completely_positive.py
+++ b/mjo/cone/completely_positive.py
@@ -12,11 +12,14 @@ from os.path import abspath
 from site import addsitedir
 addsitedir(abspath('../../'))
 from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
-
+from mjo.cone.doubly_nonnegative import is_doubly_nonnegative, is_extreme_doubly_nonnegative
 
 def is_completely_positive(A):
     """
-    Determine whether or not the matrix ``A`` is completely positive.
+    Determine whether or not the matrix ``A`` is completely
+    positive. This is a known-hard problem, and we'll just give up and
+    throw a ``ValueError`` if we can't come to a decision one way or
+    another.
 
     INPUT:
 
@@ -32,10 +35,32 @@ def is_completely_positive(A):
     Generate an extreme completely positive matrix and check that we
     identify it correctly::
 
+        sage: v = vector([1,2,3,4])
+        sage: A = v.column() * v.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_completely_positive(A)
+        True
+
+    Generate an intractible matrix (which does happen to be completely
+    positive) and explode::
+
         sage: v = vector(map(abs, random_vector(ZZ, 10)))
         sage: A = v.column() * v.row()
         sage: A = A.change_ring(QQ)
         sage: is_completely_positive(A)
+        Traceback (most recent call last):
+        ...
+        ValueError: Unable to determine extremity of ``A``.
+
+    Generate a *non*-extreme completely positive matrix and check we
+    we identify it correctly. This doesn't work so well if we generate
+    random vectors because our algorithm can give up::
+
+        sage: v1 = vector([1,2,3])
+        sage: v2 = vector([4,5,6])
+        sage: A = v1.column()*v1.row() + v2.column()*v2.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_completely_positive(A)
         True
 
     The following matrix isn't positive semidefinite, so it can't be
@@ -52,17 +77,117 @@ def is_completely_positive(A):
         sage: is_completely_positive(A)
         False
 
+    This one is symmetric but has full rank::
+
+        sage: A = matrix(QQ, [[1, 0], [0, 4]])
+        sage: is_completely_positive(A)
+        True
+
     """
 
     if A.base_ring() == SR:
-        msg = 'The matrix ``A`` cannot be the symbolic.'
+        msg = 'The matrix ``A`` cannot be symbolic.'
         raise ValueError.new(msg)
 
     if not is_symmetric_psd(A):
         return False
 
-    # Would crash if we hadn't ensured that ``A`` was symmetric
-    # positive-semidefinite.
-    X = factor_psd(A)
-    return X.rank() == 1
+    n = A.nrows() # Makes sense since ``A`` is symmetric.
+
+    if n <= 4:
+        # The two sets are the same for n <= 4.
+        return is_doubly_nonnegative(A)
+
+    # No luck.
+    raise ValueError('Unable to determine extremity of ``A``.')
+
+
+
+def is_extreme_completely_positive(A):
+    """
+    Determine whether or not the matrix ``A`` is an extreme vector of
+    the completely positive cone. This is in general a known-hard
+    problem, so our algorithm will give up and throw a ``ValueError`` if
+    a decision cannot be made one way or another.
+
+    INPUT:
+
+    - ``A`` - The matrix whose extremity we want to discover
+
+    OUTPUT:
+
+    Either ``True`` if ``A`` is an extreme vector of the completely
+    positive cone, or ``False`` otherwise.
+
+    REFERENCES:
+
+    1. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive
+       Matrices. World Scientific, 2003.
+
+    EXAMPLES:
+
+    Generate an extreme completely positive matrix and check that we
+    identify it correctly::
+
+        sage: v = vector(map(abs, random_vector(ZZ, 10)))
+        sage: A = v.column() * v.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_extreme_completely_positive(A)
+        True
+
+    Generate a *non*-extreme completely positive matrix and check we
+    we identify it correctly. This doesn't work so well if we generate
+    random vectors because our algorithm can give up::
+
+        sage: v1 = vector([1,2,3])
+        sage: v2 = vector([4,5,6])
+        sage: A = v1.column()*v1.row() + v2.column()*v2.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_extreme_completely_positive(A)
+        False
+
+    The following matrix isn't positive semidefinite, so it can't be
+    completely positive::
+
+        sage: A = matrix(QQ, [[1, 2], [2, 1]])
+        sage: is_extreme_completely_positive(A)
+        False
+
+    This matrix isn't even symmetric, so it can't be completely
+    positive::
+
+        sage: A = matrix(QQ, [[1, 2], [3, 4]])
+        sage: is_extreme_completely_positive(A)
+        False
+
+    This one is symmetric but has full rank::
+
+        sage: A = matrix(QQ, [[1, 0], [0, 4]])
+        sage: is_extreme_completely_positive(A)
+        False
+
+    """
+
+    if A.base_ring() == SR:
+        msg = 'The matrix ``A`` cannot be symbolic.'
+        raise ValueError(msg)
+
+    if not is_symmetric_psd(A):
+        return False
+
+    n = A.nrows() # Makes sense since ``A`` is symmetric.
+
+    # Easy case, see the reference (Remark 2.4).
+    if n <= 4:
+        return is_extreme_doubly_nonnegative(A)
+
+    # If n > 5, we don't really know. But we might get lucky? See the
+    # reference again for this claim.
+    if A.rank() == 1 and is_doubly_nonnegative(A):
+        return True
+
+    # We didn't get lucky. We have a characterization of the extreme
+    # vectors, but it isn't useful if we start with ``A`` because the
+    # factorization into `$XX^{T}$` may not be unique!
+    raise ValueError('Unable to determine extremity of ``A``.')