-"""
+r"""
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
-
+from sage.all import *
+from mjo.cone.symmetric_psd import factor_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:
Either ``True`` if ``A`` is completely positive, or ``False``
otherwise.
+ SETUP::
+
+ sage: from mjo.cone.completely_positive import is_completely_positive
+
EXAMPLES:
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
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):
+ if not A.is_positive_semidefinite():
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.
+
+ SETUP::
+
+ sage: from mjo.cone.completely_positive import is_extreme_completely_positive
+
+ 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 A.is_positive_semidefinite():
+ 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``.')
projects / sage.d.git / blobdiff