Refactor symmetric_pds/doubly_nonnegative.

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

diff --git a/mjo/cone/completely_positive.py b/mjo/cone/completely_positive.py
new file mode 100644 (file)
index 0000000..3353c24
--- /dev/null
+++ b/mjo/cone/completely_positive.py
@@ -0,0 +1,68 @@
+"""
+The completely positive cone `$\mathcal{K}$` over `\mathbb{R}^{n}$` is
+the set of all matrices `$A$`of the form `$\sum uu^{T}$` for `$u \in
+\mathbb{R}^{n}_{+}$`. Equivalently, `$A = XX{T}$` where all entries of
+`$X$` are nonnegative.
+"""
+
+# 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.
+from os.path import abspath
+from site import addsitedir
+addsitedir(abspath('../../'))
+from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
+
+
+def is_completely_positive(A):
+    """
+    Determine whether or not the matrix ``A`` is completely positive.
+
+    INPUT:
+
+      - ``A`` - The matrix in question
+
+    OUTPUT:
+
+    Either ``True`` if ``A`` is completely positive, or ``False``
+    otherwise.
+
+    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_completely_positive(A)
+        True
+
+    The following matrix isn't positive semidefinite, so it can't be
+    completely positive::
+
+        sage: A = matrix(QQ, [[1, 2], [2, 1]])
+        sage: is_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_completely_positive(A)
+        False
+
+    """
+
+    if A.base_ring() == SR:
+        msg = 'The matrix ``A`` cannot be the 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
+