mjo/cone/completely_positive.py

   1 """
   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 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
   9 # have to explicitly mangle our sitedir here so that "mjo.cone"
  10 # resolves.
  11 from os.path import abspath
  12 from site import addsitedir
  13 addsitedir(abspath('../../'))
  14 from mjo.cone.symmetric_psd import factor_psd, is_symmetric_psd
  15 from mjo.cone.doubly_nonnegative import is_doubly_nonnegative, is_extreme_doubly_nonnegative
  16
  17 def is_completely_positive(A):
  18     """
  19     Determine whether or not the matrix ``A`` is completely
  20     positive. This is a known-hard problem, and we'll just give up and
  21     throw a ``ValueError`` if we can't come to a decision one way or
  22     another.
  23
  24     INPUT:
  25
  26     - ``A`` - The matrix in question
  27
  28     OUTPUT:
  29
  30     Either ``True`` if ``A`` is completely positive, or ``False``
  31     otherwise.
  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     EXAMPLES:
 128
 129     Generate an extreme completely positive matrix and check that we
 130     identify it correctly::
 131
 132         sage: v = vector(map(abs, random_vector(ZZ, 10)))
 133         sage: A = v.column() * v.row()
 134         sage: A = A.change_ring(QQ)
 135         sage: is_extreme_completely_positive(A)
 136         True
 137
 138     Generate a *non*-extreme completely positive matrix and check we
 139     we identify it correctly. This doesn't work so well if we generate
 140     random vectors because our algorithm can give up::
 141
 142         sage: v1 = vector([1,2,3])
 143         sage: v2 = vector([4,5,6])
 144         sage: A = v1.column()*v1.row() + v2.column()*v2.row()
 145         sage: A = A.change_ring(QQ)
 146         sage: is_extreme_completely_positive(A)
 147         False
 148
 149     The following matrix isn't positive semidefinite, so it can't be
 150     completely positive::
 151
 152         sage: A = matrix(QQ, [[1, 2], [2, 1]])
 153         sage: is_extreme_completely_positive(A)
 154         False
 155
 156     This matrix isn't even symmetric, so it can't be completely
 157     positive::
 158
 159         sage: A = matrix(QQ, [[1, 2], [3, 4]])
 160         sage: is_extreme_completely_positive(A)
 161         False
 162
 163     This one is symmetric but has full rank::
 164
 165         sage: A = matrix(QQ, [[1, 0], [0, 4]])
 166         sage: is_extreme_completely_positive(A)
 167         False
 168
 169     """
 170
 171     if A.base_ring() == SR:
 172         msg = 'The matrix ``A`` cannot be symbolic.'
 173         raise ValueError(msg)
 174
 175     if not is_symmetric_psd(A):
 176         return False
 177
 178     n = A.nrows() # Makes sense since ``A`` is symmetric.
 179
 180     # Easy case, see the reference (Remark 2.4).
 181     if n <= 4:
 182         return is_extreme_doubly_nonnegative(A)
 183
 184     # If n > 5, we don't really know. But we might get lucky? See the
 185     # reference again for this claim.
 186     if A.rank() == 1 and is_doubly_nonnegative(A):
 187         return True
 188
 189     # We didn't get lucky. We have a characterization of the extreme
 190     # vectors, but it isn't useful if we start with ``A`` because the
 191     # factorization into `$XX^{T}$` may not be unique!
 192     raise ValueError('Unable to determine extremity of ``A``.')
 193