]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/completely_positive.py
makefile: add --long to the test command.
[sage.d.git] / 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