From 46c389b0eafe43e54afc9f81f855bfb5f889d25d Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 3 Nov 2014 09:49:02 -0500 Subject: [PATCH] Implement the is(_extreme)_completely_positive functions. --- mjo/cone/completely_positive.py | 139 ++++++++++++++++++++++++++++++-- 1 file changed, 132 insertions(+), 7 deletions(-) diff --git a/mjo/cone/completely_positive.py b/mjo/cone/completely_positive.py index 0f94ca5..ac5f144 100644 --- a/mjo/cone/completely_positive.py +++ b/mjo/cone/completely_positive.py @@ -12,11 +12,14 @@ from os.path import abspath from site import addsitedir addsitedir(abspath('../../')) from mjo.cone.symmetric_psd import factor_psd, is_symmetric_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: @@ -32,10 +35,32 @@ def is_completely_positive(A): 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,117 @@ 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): 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. + + 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 is_symmetric_psd(A): + 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``.') -- 2.43.2