]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
Add more positive operator examples from the paper.
[sage.d.git] / mjo / cone / cone.py
1 from sage.all import *
2
3 def is_lyapunov_like(L,K):
4 r"""
5 Determine whether or not ``L`` is Lyapunov-like on ``K``.
6
7 We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
8 L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
9 `\left\langle x,s \right\rangle` in the complementarity set of
10 ``K``. It is known [Orlitzky]_ that this property need only be
11 checked for generators of ``K`` and its dual.
12
13 INPUT:
14
15 - ``L`` -- A linear transformation or matrix.
16
17 - ``K`` -- A polyhedral closed convex cone.
18
19 OUTPUT:
20
21 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
22 and ``False`` otherwise.
23
24 .. WARNING::
25
26 If this function returns ``True``, then ``L`` is Lyapunov-like
27 on ``K``. However, if ``False`` is returned, that could mean one
28 of two things. The first is that ``L`` is definitely not
29 Lyapunov-like on ``K``. The second is more of an "I don't know"
30 answer, returned (for example) if we cannot prove that an inner
31 product is zero.
32
33 REFERENCES:
34
35 M. Orlitzky. The Lyapunov rank of an improper cone.
36 http://www.optimization-online.org/DB_HTML/2015/10/5135.html
37
38 EXAMPLES:
39
40 The identity is always Lyapunov-like in a nontrivial space::
41
42 sage: set_random_seed()
43 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
44 sage: L = identity_matrix(K.lattice_dim())
45 sage: is_lyapunov_like(L,K)
46 True
47
48 As is the "zero" transformation::
49
50 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
51 sage: R = K.lattice().vector_space().base_ring()
52 sage: L = zero_matrix(R, K.lattice_dim())
53 sage: is_lyapunov_like(L,K)
54 True
55
56 Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
57 on ``K``::
58
59 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
60 sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
61 True
62
63 """
64 return all([(L*x).inner_product(s) == 0
65 for (x,s) in K.discrete_complementarity_set()])
66
67
68 def motzkin_decomposition(K):
69 r"""
70 Return the pair of components in the Motzkin decomposition of this cone.
71
72 Every convex cone is the direct sum of a strictly convex cone and a
73 linear subspace [Stoer-Witzgall]_. Return a pair ``(P,S)`` of cones
74 such that ``P`` is strictly convex, ``S`` is a subspace, and ``K``
75 is the direct sum of ``P`` and ``S``.
76
77 OUTPUT:
78
79 An ordered pair ``(P,S)`` of closed convex polyhedral cones where
80 ``P`` is strictly convex, ``S`` is a subspace, and ``K`` is the
81 direct sum of ``P`` and ``S``.
82
83 REFERENCES:
84
85 .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and
86 Optimization in Finite Dimensions I. Springer-Verlag, New
87 York, 1970.
88
89 EXAMPLES:
90
91 The nonnegative orthant is strictly convex, so it is its own
92 strictly convex component and its subspace component is trivial::
93
94 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
95 sage: (P,S) = motzkin_decomposition(K)
96 sage: K.is_equivalent(P)
97 True
98 sage: S.is_trivial()
99 True
100
101 Likewise, full spaces are their own subspace components::
102
103 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
104 sage: K.is_full_space()
105 True
106 sage: (P,S) = motzkin_decomposition(K)
107 sage: K.is_equivalent(S)
108 True
109 sage: P.is_trivial()
110 True
111
112 TESTS:
113
114 A random point in the cone should belong to either the strictly
115 convex component or the subspace component. If the point is nonzero,
116 it cannot be in both::
117
118 sage: set_random_seed()
119 sage: K = random_cone(max_ambient_dim=8)
120 sage: (P,S) = motzkin_decomposition(K)
121 sage: x = K.random_element()
122 sage: P.contains(x) or S.contains(x)
123 True
124 sage: x.is_zero() or (P.contains(x) != S.contains(x))
125 True
126
127 The strictly convex component should always be strictly convex, and
128 the subspace component should always be a subspace::
129
130 sage: set_random_seed()
131 sage: K = random_cone(max_ambient_dim=8)
132 sage: (P,S) = motzkin_decomposition(K)
133 sage: P.is_strictly_convex()
134 True
135 sage: S.lineality() == S.dim()
136 True
137
138 The generators of the components are obtained from orthogonal
139 projections of the original generators [Stoer-Witzgall]_::
140
141 sage: set_random_seed()
142 sage: K = random_cone(max_ambient_dim=8)
143 sage: (P,S) = motzkin_decomposition(K)
144 sage: A = S.linear_subspace().complement().matrix()
145 sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A
146 sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice())
147 sage: P.is_equivalent(expected_P)
148 True
149 sage: A = S.linear_subspace().matrix()
150 sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A
151 sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice())
152 sage: S.is_equivalent(expected_S)
153 True
154 """
155 # The lines() method only returns one generator per line. For a true
156 # line, we also need a generator pointing in the opposite direction.
157 S_gens = [ direction*gen for direction in [1,-1] for gen in K.lines() ]
158 S = Cone(S_gens, K.lattice())
159
160 # Since ``S`` is a subspace, the rays of its dual generate its
161 # orthogonal complement.
162 S_perp = Cone(S.dual(), K.lattice())
163 P = K.intersection(S_perp)
164
165 return (P,S)
166
167
168 def positive_operator_gens(K):
169 r"""
170 Compute generators of the cone of positive operators on this cone.
171
172 OUTPUT:
173
174 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
175 Each matrix ``P`` in the list should have the property that ``P*x``
176 is an element of ``K`` whenever ``x`` is an element of
177 ``K``. Moreover, any nonnegative linear combination of these
178 matrices shares the same property.
179
180 EXAMPLES:
181
182 Positive operators on the nonnegative orthant are nonnegative matrices::
183
184 sage: K = Cone([(1,)])
185 sage: positive_operator_gens(K)
186 [[1]]
187
188 sage: K = Cone([(1,0),(0,1)])
189 sage: positive_operator_gens(K)
190 [
191 [1 0] [0 1] [0 0] [0 0]
192 [0 0], [0 0], [1 0], [0 1]
193 ]
194
195 The trivial cone in a trivial space has no positive operators::
196
197 sage: K = Cone([], ToricLattice(0))
198 sage: positive_operator_gens(K)
199 []
200
201 Every operator is positive on the trivial cone::
202
203 sage: K = Cone([(0,)])
204 sage: positive_operator_gens(K)
205 [[1], [-1]]
206
207 sage: K = Cone([(0,0)])
208 sage: K.is_trivial()
209 True
210 sage: positive_operator_gens(K)
211 [
212 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
213 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
214 ]
215
216 Every operator is positive on the ambient vector space::
217
218 sage: K = Cone([(1,),(-1,)])
219 sage: K.is_full_space()
220 True
221 sage: positive_operator_gens(K)
222 [[1], [-1]]
223
224 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
225 sage: K.is_full_space()
226 True
227 sage: positive_operator_gens(K)
228 [
229 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
230 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
231 ]
232
233 A non-obvious application is to find the positive operators on the
234 right half-plane::
235
236 sage: K = Cone([(1,0),(0,1),(0,-1)])
237 sage: positive_operator_gens(K)
238 [
239 [1 0] [0 0] [ 0 0] [0 0] [ 0 0]
240 [0 0], [1 0], [-1 0], [0 1], [ 0 -1]
241 ]
242
243 TESTS:
244
245 Each positive operator generator should send the generators of the
246 cone into the cone::
247
248 sage: set_random_seed()
249 sage: K = random_cone(max_ambient_dim=5)
250 sage: pi_of_K = positive_operator_gens(K)
251 sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
252 True
253
254 Each positive operator generator should send a random element of the
255 cone into the cone::
256
257 sage: set_random_seed()
258 sage: K = random_cone(max_ambient_dim=5)
259 sage: pi_of_K = positive_operator_gens(K)
260 sage: all([ K.contains(P*K.random_element()) for P in pi_of_K ])
261 True
262
263 A random element of the positive operator cone should send the
264 generators of the cone into the cone::
265
266 sage: set_random_seed()
267 sage: K = random_cone(max_ambient_dim=5)
268 sage: pi_of_K = positive_operator_gens(K)
269 sage: L = ToricLattice(K.lattice_dim()**2)
270 sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
271 sage: P = matrix(K.lattice_dim(), pi_cone.random_element().list())
272 sage: all([ K.contains(P*x) for x in K ])
273 True
274
275 A random element of the positive operator cone should send a random
276 element of the cone into the cone::
277
278 sage: set_random_seed()
279 sage: K = random_cone(max_ambient_dim=5)
280 sage: pi_of_K = positive_operator_gens(K)
281 sage: L = ToricLattice(K.lattice_dim()**2)
282 sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
283 sage: P = matrix(K.lattice_dim(), pi_cone.random_element().list())
284 sage: K.contains(P*K.random_element())
285 True
286
287 The dimension of the cone of positive operators is given by the
288 corollary in my paper::
289
290 sage: set_random_seed()
291 sage: K = random_cone(max_ambient_dim=5)
292 sage: n = K.lattice_dim()
293 sage: m = K.dim()
294 sage: l = K.lineality()
295 sage: pi_of_K = positive_operator_gens(K)
296 sage: L = ToricLattice(n**2)
297 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
298 sage: expected = n**2 - l*(m - l) - (n - m)*m
299 sage: actual == expected
300 True
301
302 The lineality of the cone of positive operators is given by the
303 corollary in my paper::
304
305 sage: set_random_seed()
306 sage: K = random_cone(max_ambient_dim=5)
307 sage: n = K.lattice_dim()
308 sage: pi_of_K = positive_operator_gens(K)
309 sage: L = ToricLattice(n**2)
310 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
311 sage: expected = n**2 - K.dim()*K.dual().dim()
312 sage: actual == expected
313 True
314
315 The cone ``K`` is proper if and only if the cone of positive
316 operators on ``K`` is proper::
317
318 sage: set_random_seed()
319 sage: K = random_cone(max_ambient_dim=5)
320 sage: pi_of_K = positive_operator_gens(K)
321 sage: L = ToricLattice(K.lattice_dim()**2)
322 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
323 sage: K.is_proper() == pi_cone.is_proper()
324 True
325 """
326 # Matrices are not vectors in Sage, so we have to convert them
327 # to vectors explicitly before we can find a basis. We need these
328 # two values to construct the appropriate "long vector" space.
329 F = K.lattice().base_field()
330 n = K.lattice_dim()
331
332 tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
333
334 # Convert those tensor products to long vectors.
335 W = VectorSpace(F, n**2)
336 vectors = [ W(tp.list()) for tp in tensor_products ]
337
338 # Create the *dual* cone of the positive operators, expressed as
339 # long vectors..
340 pi_dual = Cone(vectors, ToricLattice(W.dimension()))
341
342 # Now compute the desired cone from its dual...
343 pi_cone = pi_dual.dual()
344
345 # And finally convert its rays back to matrix representations.
346 M = MatrixSpace(F, n)
347 return [ M(v.list()) for v in pi_cone.rays() ]
348
349
350 def Z_transformation_gens(K):
351 r"""
352 Compute generators of the cone of Z-transformations on this cone.
353
354 OUTPUT:
355
356 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
357 Each matrix ``L`` in the list should have the property that
358 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
359 discrete complementarity set of ``K``. Moreover, any nonnegative
360 linear combination of these matrices shares the same property.
361
362 EXAMPLES:
363
364 Z-transformations on the nonnegative orthant are just Z-matrices.
365 That is, matrices whose off-diagonal elements are nonnegative::
366
367 sage: K = Cone([(1,0),(0,1)])
368 sage: Z_transformation_gens(K)
369 [
370 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
371 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
372 ]
373 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
374 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
375 ....: for i in range(z.nrows())
376 ....: for j in range(z.ncols())
377 ....: if i != j ])
378 True
379
380 The trivial cone in a trivial space has no Z-transformations::
381
382 sage: K = Cone([], ToricLattice(0))
383 sage: Z_transformation_gens(K)
384 []
385
386 Z-transformations on a subspace are Lyapunov-like and vice-versa::
387
388 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
389 sage: K.is_full_space()
390 True
391 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
392 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
393 sage: zs == lls
394 True
395
396 TESTS:
397
398 The Z-property is possessed by every Z-transformation::
399
400 sage: set_random_seed()
401 sage: K = random_cone(max_ambient_dim=6)
402 sage: Z_of_K = Z_transformation_gens(K)
403 sage: dcs = K.discrete_complementarity_set()
404 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
405 ....: for (x,s) in dcs])
406 True
407
408 The lineality space of Z is LL::
409
410 sage: set_random_seed()
411 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
412 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
413 sage: z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ])
414 sage: z_cone.linear_subspace() == lls
415 True
416
417 And thus, the lineality of Z is the Lyapunov rank::
418
419 sage: set_random_seed()
420 sage: K = random_cone(max_ambient_dim=6)
421 sage: Z_of_K = Z_transformation_gens(K)
422 sage: L = ToricLattice(K.lattice_dim()**2)
423 sage: z_cone = Cone([ z.list() for z in Z_of_K ], lattice=L)
424 sage: z_cone.lineality() == K.lyapunov_rank()
425 True
426
427 The lineality spaces of pi-star and Z-star are equal:
428
429 sage: set_random_seed()
430 sage: K = random_cone(max_ambient_dim=5)
431 sage: pi_of_K = positive_operator_gens(K)
432 sage: Z_of_K = Z_transformation_gens(K)
433 sage: L = ToricLattice(K.lattice_dim()**2)
434 sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
435 sage: z_star = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
436 sage: pi_star.linear_subspace() == z_star.linear_subspace()
437 True
438 """
439 # Matrices are not vectors in Sage, so we have to convert them
440 # to vectors explicitly before we can find a basis. We need these
441 # two values to construct the appropriate "long vector" space.
442 F = K.lattice().base_field()
443 n = K.lattice_dim()
444
445 # These tensor products contain generators for the dual cone of
446 # the cross-positive transformations.
447 tensor_products = [ s.tensor_product(x)
448 for (x,s) in K.discrete_complementarity_set() ]
449
450 # Turn our matrices into long vectors...
451 W = VectorSpace(F, n**2)
452 vectors = [ W(m.list()) for m in tensor_products ]
453
454 # Create the *dual* cone of the cross-positive operators,
455 # expressed as long vectors..
456 Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
457
458 # Now compute the desired cone from its dual...
459 Sigma_cone = Sigma_dual.dual()
460
461 # And finally convert its rays back to matrix representations.
462 # But first, make them negative, so we get Z-transformations and
463 # not cross-positive ones.
464 M = MatrixSpace(F, n)
465 return [ -M(v.list()) for v in Sigma_cone.rays() ]
466
467
468 def Z_cone(K):
469 gens = Z_transformation_gens(K)
470 L = None
471 if len(gens) == 0:
472 L = ToricLattice(0)
473 return Cone([ g.list() for g in gens ], lattice=L)
474
475 def pi_cone(K):
476 gens = positive_operator_gens(K)
477 L = None
478 if len(gens) == 0:
479 L = ToricLattice(0)
480 return Cone([ g.list() for g in gens ], lattice=L)