--- /dev/null
+"""
+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
+
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
"""
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)
from mjo.symbolic import matrix_simplify_full
+def is_symmetric_psd(A):
+ """
+ Determine whether or not the matrix ``A`` is symmetric
+ positive-semidefinite.
+
+ INPUT:
+
+ - ``A`` - The matrix in question
+
+ OUTPUT:
+
+ Either ``True`` if ``A`` is symmetric positive-semidefinite, or
+ ``False`` otherwise.
+
+ EXAMPLES:
+
+ Every completely positive matrix is symmetric
+ positive-semidefinite::
+
+ sage: v = vector(map(abs, random_vector(ZZ, 10)))
+ sage: A = v.column() * v.row()
+ sage: is_symmetric_psd(A)
+ True
+
+ The following matrix is symmetric but not positive semidefinite::
+
+ sage: A = matrix(ZZ, [[1, 2], [2, 1]])
+ sage: is_symmetric_psd(A)
+ False
+
+ This matrix isn't even symmetric::
+
+ sage: A = matrix(ZZ, [[1, 2], [3, 4]])
+ sage: is_symmetric_psd(A)
+ False
+
+ """
+
+ if A.base_ring() == SR:
+ msg = 'The matrix ``A`` cannot be symbolic.'
+ raise ValueError.new(msg)
+
+ # First make sure that ``A`` is symmetric.
+ if not A.is_symmetric():
+ return False
+
+ # If ``A`` is symmetric, 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
+
+
def unit_eigenvectors(A):
"""
Return the unit eigenvectors of a symmetric positive-definite matrix.
INPUT:
- - ``A`` - The matrix to factor.
+ - ``A`` - The matrix to factor. The base ring of ``A`` must be either
+ exact or the symbolic ring (to compute eigenvalues), and it
+ must be a field so that we can take its algebraic closure
+ (necessary to e.g. take square roots).
OUTPUT:
- A matrix ``X`` such that `A = XX^{T}`.
+ A matrix ``X`` such that `A = XX^{T}`. The base field of ``X`` will
+ be the algebraic closure of the base field of ``A``.
ALGORITHM:
`$D$` will have dimension `$k \times k$`. In the end everything
works out the same.
- EXAMPLES::
+ EXAMPLES:
+
+ Create a symmetric positive-semidefinite matrix over the symbolic
+ ring and factor it::
sage: A = matrix(SR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
sage: X = factor_psd(A)
sage: A == A2
True
+ Attempt to factor the same matrix over ``RR`` which won't work
+ because ``RR`` isn't exact::
+
+ sage: A = matrix(RR, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
+ sage: factor_psd(A)
+ Traceback (most recent call last):
+ ...
+ ValueError: The base ring of ``A`` must be either exact or symbolic.
+
+ Attempt to factor the same matrix over ``ZZ`` which won't work
+ because ``ZZ`` isn't a field::
+
+ sage: A = matrix(ZZ, [[0, 2, 3], [2, 0, 0], [3, 0, 0]])
+ sage: factor_psd(A)
+ Traceback (most recent call last):
+ ...
+ ValueError: The base ring of ``A`` must be a field.
+
"""
+ if not A.base_ring().is_exact() and not A.base_ring() is SR:
+ raise ValueError('The base ring of ``A`` must be either exact or symbolic.')
+
+ 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.
+ A = A.change_ring(A.base_ring().algebraic_closure())
+
+
# Get the eigenvectors, and filter out the ones that correspond to
# the eigenvalue zero.
all_evs = unit_eigenvectors(A)
projects / sage.d.git / commitdiff