mjo/cone/completely_positive.py

   1 r"""
   2 The completely positive cone `$\mathcal{K}$` over `\mathbb{R}^{n}$` is
   3 the set of all matrices `$A$`of the form `$\sum uu^{T}$` for `$u \in
   4 \mathbb{R}^{n}_{+}$`. Equivalently, `$A = XX{T}$` where all entries of
   5 `$X$` are nonnegative.
   6 """
   7
   8 from sage.all import *
   9 from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
  10 from mjo.cone.doubly_nonnegative import (is_doubly_nonnegative,
  11                                          is_extreme_doubly_nonnegative)
  12
  13 def is_completely_positive(A):
  14     """
  15     Determine whether or not the matrix ``A`` is completely
  16     positive. This is a known-hard problem, and we'll just give up and
  17     throw a ``ValueError`` if we can't come to a decision one way or
  18     another.
  19
  20     INPUT:
  21
  22     - ``A`` - The matrix in question
  23
  24     OUTPUT:
  25
  26     Either ``True`` if ``A`` is completely positive, or ``False``
  27     otherwise.
  28
  29     SETUP::
  30
  31         sage: from mjo.cone.completely_positive import is_completely_positive
  32
  33     EXAMPLES:
  34
  35     Generate an extreme completely positive matrix and check that we
  36     identify it correctly::
  37
  38         sage: v = vector([1,2,3,4])
  39         sage: A = v.column() * v.row()
  40         sage: A = A.change_ring(QQ)
  41         sage: is_completely_positive(A)
  42         True
  43
  44     Generate an intractible matrix (which does happen to be completely
  45     positive) and explode::
  46
  47         sage: v = vector(map(abs, random_vector(ZZ, 10)))
  48         sage: A = v.column() * v.row()
  49         sage: A = A.change_ring(QQ)
  50         sage: is_completely_positive(A)
  51         Traceback (most recent call last):
  52         ...
  53         ValueError: Unable to determine extremity of ``A``.
  54
  55     Generate a *non*-extreme completely positive matrix and check we
  56     we identify it correctly. This doesn't work so well if we generate
  57     random vectors because our algorithm can give up::
  58
  59         sage: v1 = vector([1,2,3])
  60         sage: v2 = vector([4,5,6])
  61         sage: A = v1.column()*v1.row() + v2.column()*v2.row()
  62         sage: A = A.change_ring(QQ)
  63         sage: is_completely_positive(A)
  64         True
  65
  66     The following matrix isn't positive semidefinite, so it can't be
  67     completely positive::
  68
  69         sage: A = matrix(QQ, [[1, 2], [2, 1]])
  70         sage: is_completely_positive(A)
  71         False
  72
  73     This matrix isn't even symmetric, so it can't be completely
  74     positive::
  75
  76         sage: A = matrix(QQ, [[1, 2], [3, 4]])
  77         sage: is_completely_positive(A)
  78         False
  79
  80     This one is symmetric but has full rank::
  81
  82         sage: A = matrix(QQ, [[1, 0], [0, 4]])
  83         sage: is_completely_positive(A)
  84         True
  85
  86     """
  87
  88     if A.base_ring() == SR:
  89         msg = 'The matrix ``A`` cannot be symbolic.'
  90         raise ValueError.new(msg)
  91
  92     if not is_symmetric_psd(A):
  93         return False
  94
  95     n = A.nrows() # Makes sense since ``A`` is symmetric.
  96
  97     if n <= 4:
  98         # The two sets are the same for n <= 4.
  99         return is_doubly_nonnegative(A)
 100
 101     # No luck.
 102     raise ValueError('Unable to determine extremity of ``A``.')
 103
 104
 105
 106 def is_extreme_completely_positive(A):
 107     """
 108     Determine whether or not the matrix ``A`` is an extreme vector of
 109     the completely positive cone. This is in general a known-hard
 110     problem, so our algorithm will give up and throw a ``ValueError`` if
 111     a decision cannot be made one way or another.
 112
 113     INPUT:
 114
 115     - ``A`` - The matrix whose extremity we want to discover
 116
 117     OUTPUT:
 118
 119     Either ``True`` if ``A`` is an extreme vector of the completely
 120     positive cone, or ``False`` otherwise.
 121
 122     REFERENCES:
 123
 124     1. Berman, Abraham and Shaked-Monderer, Naomi. Completely Positive
 125        Matrices. World Scientific, 2003.
 126
 127     SETUP::
 128
 129         sage: from mjo.cone.completely_positive import is_extreme_completely_positive
 130
 131     EXAMPLES:
 132
 133     Generate an extreme completely positive matrix and check that we
 134     identify it correctly::
 135
 136         sage: v = vector(map(abs, random_vector(ZZ, 10)))
 137         sage: A = v.column() * v.row()
 138         sage: A = A.change_ring(QQ)
 139         sage: is_extreme_completely_positive(A)
 140         True
 141
 142     Generate a *non*-extreme completely positive matrix and check we
 143     we identify it correctly. This doesn't work so well if we generate
 144     random vectors because our algorithm can give up::
 145
 146         sage: v1 = vector([1,2,3])
 147         sage: v2 = vector([4,5,6])
 148         sage: A = v1.column()*v1.row() + v2.column()*v2.row()
 149         sage: A = A.change_ring(QQ)
 150         sage: is_extreme_completely_positive(A)
 151         False
 152
 153     The following matrix isn't positive semidefinite, so it can't be
 154     completely positive::
 155
 156         sage: A = matrix(QQ, [[1, 2], [2, 1]])
 157         sage: is_extreme_completely_positive(A)
 158         False
 159
 160     This matrix isn't even symmetric, so it can't be completely
 161     positive::
 162
 163         sage: A = matrix(QQ, [[1, 2], [3, 4]])
 164         sage: is_extreme_completely_positive(A)
 165         False
 166
 167     This one is symmetric but has full rank::
 168
 169         sage: A = matrix(QQ, [[1, 0], [0, 4]])
 170         sage: is_extreme_completely_positive(A)
 171         False
 172
 173     """
 174
 175     if A.base_ring() == SR:
 176         msg = 'The matrix ``A`` cannot be symbolic.'
 177         raise ValueError(msg)
 178
 179     if not is_symmetric_psd(A):
 180         return False
 181
 182     n = A.nrows() # Makes sense since ``A`` is symmetric.
 183
 184     # Easy case, see the reference (Remark 2.4).
 185     if n <= 4:
 186         return is_extreme_doubly_nonnegative(A)
 187
 188     # If n > 5, we don't really know. But we might get lucky? See the
 189     # reference again for this claim.
 190     if A.rank() == 1 and is_doubly_nonnegative(A):
 191         return True
 192
 193     # We didn't get lucky. We have a characterization of the extreme
 194     # vectors, but it isn't useful if we start with ``A`` because the
 195     # factorization into `$XX^{T}$` may not be unique!
 196     raise ValueError('Unable to determine extremity of ``A``.')
 197