Implement the is_extreme_doubly_nonnegative() function, first shot.

diff --git a/mjo/cone/doubly_nonnegative.py b/mjo/cone/doubly_nonnegative.py
index 03d23b4ddcbd785b48fe9dd447876cefb586a067..4f7f950bb2c897509c2cb4b41d5200bccfd67c49 100644 (file)
--- a/mjo/cone/doubly_nonnegative.py
+++ b/mjo/cone/doubly_nonnegative.py
@@ -14,13 +14,13 @@ It is represented typically by either `\mathcal{D}^{n}` or
 from sage.all import *
 
 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
-# have to explicitly mangle our sitedir here so that "mjo.cone"
-# resolves.
+# have to explicitly mangle our sitedir here so that our module names
+# resolve.
 from os.path import abspath
 from site import addsitedir
 addsitedir(abspath('../../'))
 from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
-
+from mjo.matrix_vector import isomorphism
 
 
 def is_doubly_nonnegative(A):
@@ -188,6 +188,16 @@ def is_extreme_doubly_nonnegative(A):
     Returns ``True`` if the given matrix is an extreme matrix of the
     doubly-nonnegative cone, and ``False`` otherwise.
 
+    REFERENCES:
+
+    1. 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.
+
+    2. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive
+       Matrices. World Scientific, 2003.
+
     EXAMPLES:
 
     The zero matrix is an extreme matrix::
@@ -196,15 +206,101 @@ def is_extreme_doubly_nonnegative(A):
         sage: is_extreme_doubly_nonnegative(A)
         True
 
+    Any extreme vector of the completely positive cone is an extreme
+    vector of the doubly-nonnegative cone::
+
+        sage: v = vector([1,2,3,4,5,6])
+        sage: A = v.column() * v.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_extreme_doubly_nonnegative(A)
+        True
+
+    We should be able to generate the extreme completely positive
+    vectors randomly::
+
+        sage: v = vector(map(abs, random_vector(ZZ, 4)))
+        sage: A = v.column() * v.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_extreme_doubly_nonnegative(A)
+        True
+        sage: v = vector(map(abs, random_vector(ZZ, 10)))
+        sage: A = v.column() * v.row()
+        sage: A = A.change_ring(QQ)
+        sage: is_extreme_doubly_nonnegative(A)
+        True
+
+    The following matrix is completely positive but has rank 3, so by a
+    remark in reference #1 it is not extreme::
+
+        sage: A = matrix(QQ, [[1,2,1],[2,6,3],[1,3,5]])
+        sage: is_extreme_doubly_nonnegative(A)
+        False
+
+    The following matrix is completely positive (diagonal) with rank 2,
+    so it is also not extreme::
+
+        sage: A = matrix(QQ, [[1,0,0],[2,0,0],[0,0,0]])
+        sage: is_extreme_doubly_nonnegative(A)
+        False
+
     """
 
-    r = A.rank()
+    if not A.base_ring().is_exact() and not A.base_ring() is SR:
+        msg = 'The base ring of ``A`` must be either exact or symbolic.'
+        raise ValueError(msg)
 
-    if r == 0:
+    if not A.base_ring().is_field():
+        raise ValueError('The base ring of ``A`` must be a field.')
+
+    if not A.base_ring() is SR:
+        # Change the base field of ``A`` so that we are sure we can take
+        # roots. The symbolic ring has no algebraic_closure method.
+        A = A.change_ring(A.base_ring().algebraic_closure())
+
+    # Step 1 (see reference #1)
+    k = A.rank()
+
+    if k == 0:
         # Short circuit, we know the zero matrix is extreme.
         return True
 
-    if not is_admissible_extreme_rank(r):
+    if not is_symmetric_psd(A):
+        return False
+
+    # Step 1.5, appeal to Theorem 3.1 in reference #1 to short
+    # circuit.
+    if not has_admissible_extreme_rank(A):
         return False
 
-    raise NotImplementedError()
+    # Step 2
+    X = factor_psd(A)
+
+    # Step 3
+    #
+    # Begin with an empty spanning set, and add a new matrix to it
+    # whenever we come across an index pair `$(i,j)$` with
+    # `$A_{ij} = 0$`.
+    spanning_set = []
+    for j in range(0, A.ncols()):
+        for i in range(0,j):
+            if A[i,j] == 0:
+                M = A.matrix_space()
+                S = X.transpose() * (E(M,i,j) + E(M,j,i)) * X
+                spanning_set.append(S)
+
+    # The spanning set that we have at this point is of matrices.  We
+    # only care about the dimension of the spanned space, and Sage
+    # can't compute the dimension of a set of matrices anyway, so we
+    # convert them all to vectors and just ask for the dimension of the
+    # resulting vector space.
+    (phi, phi_inverse) = isomorphism(A.matrix_space())
+    vectors = map(phi,spanning_set)
+
+    V = span(vectors, A.base_ring())
+    d = V.dimension()
+
+    # Needed to safely divide by two here (we don't want integer
+    # division). We ensured that the base ring of ``A`` is a field
+    # earlier.
+    two = A.base_ring()(2)
+    return d == (k*(k + 1)/two - 1)