""" 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. """ from sage.all import * 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. 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: - ``A`` - The matrix in question OUTPUT: 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 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 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 symbolic.' raise ValueError.new(msg) if not is_symmetric_psd(A): return False 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 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``.') def completely_positive_operators_gens(K): r""" Return a list of generators (matrices) for the completely-positive cone of ``K``. INPUT: - ``K`` -- a closed convex rational polyhedral cone. OUTPUT: A list of matrices, the conic hull of which is the completely-positive cone of ``K``. SETUP:: sage: from mjo.cone.completely_positive import ( ....: completely_positive_operators_gens, ....: is_completely_positive ) sage: from mjo.cone.nonnegative_orthant import nonnegative_orthant sage: from mjo.matrix_vector import isomorphism EXAMPLES:: sage: K = nonnegative_orthant(2) sage: completely_positive_operators_gens(K) [ [1 0] [0 0] [0 0], [0 1] ] sage: all( is_completely_positive(M) ....: for M in completely_positive_operators_gens(K) ) True TESTS: The completely-positive cone of ``K`` is subdual:: sage: K = random_cone(max_ambient_dim=8, max_rays=10) sage: cp_gens = completely_positive_operators_gens(K) sage: n = K.lattice_dim() sage: M = MatrixSpace(QQ, n, n) sage: (p, p_inv) = isomorphism(M) sage: L = ToricLattice(n**2) sage: cp_cone = Cone( (p(m) for m in cp_gens), lattice=L ) sage: copos_cone = Cone(cp_cone.dual().rays(), lattice=L ) sage: all( x in copos_cone for x in cp_cone ) True """ return [ x.tensor_product(x) for x in K ]