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