X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcompletely_positive.py;h=3dc66b454dfb8a85a4aa2ba4a4be03978a5b3690;hb=7af2b9d146a6bf2fb8acc3c342983de577b417ce;hp=0f94ca50e9d42d3540a7bb34e84c8ac2b9f57b55;hpb=3a800989b6e935fdab52418ac644d926997b83d0;p=sage.d.git diff --git a/mjo/cone/completely_positive.py b/mjo/cone/completely_positive.py index 0f94ca5..3dc66b4 100644 --- a/mjo/cone/completely_positive.py +++ b/mjo/cone/completely_positive.py @@ -1,22 +1,21 @@ -""" +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: @@ -27,15 +26,41 @@ def is_completely_positive(A): 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 @@ -52,17 +77,121 @@ 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): + 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``.')