]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
adba809cf26722f52b0ec75432b2e6068683c3d9
3 def is_lyapunov_like(L
,K
):
5 Determine whether or not ``L`` is Lyapunov-like on ``K``.
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.
15 - ``L`` -- A linear transformation or matrix.
17 - ``K`` -- A polyhedral closed convex cone.
21 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
22 and ``False`` otherwise.
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
35 M. Orlitzky. The Lyapunov rank of an improper cone.
36 http://www.optimization-online.org/DB_HTML/2015/10/5135.html
40 The identity is always Lyapunov-like in a nontrivial space::
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)
48 As is the "zero" transformation::
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)
56 Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
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() ])
64 return all([(L
*x
).inner_product(s
) == 0
65 for (x
,s
) in K
.discrete_complementarity_set()])
68 def motzkin_decomposition(K
):
70 Return the pair of components in the Motzkin decomposition of this cone.
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``.
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``.
85 .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and
86 Optimization in Finite Dimensions I. Springer-Verlag, New
91 The nonnegative orthant is strictly convex, so it is its own
92 strictly convex component and its subspace component is trivial::
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)
101 Likewise, full spaces are their own subspace components::
103 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
104 sage: K.is_full_space()
106 sage: (P,S) = motzkin_decomposition(K)
107 sage: K.is_equivalent(S)
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::
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(ring=QQ)
122 sage: P.contains(x) or S.contains(x)
124 sage: x.is_zero() or (P.contains(x) != S.contains(x))
127 The strictly convex component should always be strictly convex, and
128 the subspace component should always be a subspace::
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()
135 sage: S.lineality() == S.dim()
138 The generators of the components are obtained from orthogonal
139 projections of the original generators [Stoer-Witzgall]_::
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)
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)
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())
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
)
168 def positive_operator_gens(K
):
170 Compute generators of the cone of positive operators on this cone.
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.
182 Positive operators on the nonnegative orthant are nonnegative matrices::
184 sage: K = Cone([(1,)])
185 sage: positive_operator_gens(K)
188 sage: K = Cone([(1,0),(0,1)])
189 sage: positive_operator_gens(K)
191 [1 0] [0 1] [0 0] [0 0]
192 [0 0], [0 0], [1 0], [0 1]
195 The trivial cone in a trivial space has no positive operators::
197 sage: K = Cone([], ToricLattice(0))
198 sage: positive_operator_gens(K)
201 Every operator is positive on the trivial cone::
203 sage: K = Cone([(0,)])
204 sage: positive_operator_gens(K)
207 sage: K = Cone([(0,0)])
210 sage: positive_operator_gens(K)
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]
216 Every operator is positive on the ambient vector space::
218 sage: K = Cone([(1,),(-1,)])
219 sage: K.is_full_space()
221 sage: positive_operator_gens(K)
224 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
225 sage: K.is_full_space()
227 sage: positive_operator_gens(K)
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]
233 A non-obvious application is to find the positive operators on the
236 sage: K = Cone([(1,0),(0,1),(0,-1)])
237 sage: positive_operator_gens(K)
239 [1 0] [0 0] [ 0 0] [0 0] [ 0 0]
240 [0 0], [1 0], [-1 0], [0 1], [ 0 -1]
245 Each positive operator generator should send the generators of the
248 sage: set_random_seed()
249 sage: K = random_cone(max_ambient_dim=4)
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 ])
254 Each positive operator generator should send a random element of the
257 sage: set_random_seed()
258 sage: K = random_cone(max_ambient_dim=4)
259 sage: pi_of_K = positive_operator_gens(K)
260 sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
263 A random element of the positive operator cone should send the
264 generators of the cone into the cone::
266 sage: set_random_seed()
267 sage: K = random_cone(max_ambient_dim=4)
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 ],
273 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
274 sage: all([ K.contains(P*x) for x in K ])
277 A random element of the positive operator cone should send a random
278 element of the cone into the cone::
280 sage: set_random_seed()
281 sage: K = random_cone(max_ambient_dim=4)
282 sage: pi_of_K = positive_operator_gens(K)
283 sage: L = ToricLattice(K.lattice_dim()**2)
284 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
287 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
288 sage: K.contains(P*K.random_element(ring=QQ))
291 The lineality space of the dual of the cone of positive operators
292 can be computed from the lineality spaces of the cone and its dual::
294 sage: set_random_seed()
295 sage: K = random_cone(max_ambient_dim=4)
296 sage: pi_of_K = positive_operator_gens(K)
297 sage: L = ToricLattice(K.lattice_dim()**2)
298 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
301 sage: actual = pi_cone.dual().linear_subspace()
302 sage: U1 = [ vector((s.tensor_product(x)).list())
303 ....: for x in K.lines()
304 ....: for s in K.dual() ]
305 sage: U2 = [ vector((s.tensor_product(x)).list())
307 ....: for s in K.dual().lines() ]
308 sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
309 sage: actual == expected
312 The lineality of the dual of the cone of positive operators
313 is known from its lineality space::
315 sage: set_random_seed()
316 sage: K = random_cone(max_ambient_dim=4)
317 sage: n = K.lattice_dim()
319 sage: l = K.lineality()
320 sage: pi_of_K = positive_operator_gens(K)
321 sage: L = ToricLattice(n**2)
322 sage: pi_cone = Cone([p.list() for p in pi_of_K],
325 sage: actual = pi_cone.dual().lineality()
326 sage: expected = l*(m - l) + m*(n - m)
327 sage: actual == expected
330 The dimension of the cone of positive operators is given by the
331 corollary in my paper::
333 sage: set_random_seed()
334 sage: K = random_cone(max_ambient_dim=4)
335 sage: n = K.lattice_dim()
337 sage: l = K.lineality()
338 sage: pi_of_K = positive_operator_gens(K)
339 sage: L = ToricLattice(n**2)
340 sage: pi_cone = Cone([p.list() for p in pi_of_K],
343 sage: actual = pi_cone.dim()
344 sage: expected = n**2 - l*(m - l) - (n - m)*m
345 sage: actual == expected
348 The trivial cone, full space, and half-plane all give rise to the
349 expected dimensions::
351 sage: n = ZZ.random_element().abs()
352 sage: K = Cone([[0] * n], ToricLattice(n))
355 sage: L = ToricLattice(n^2)
356 sage: pi_of_K = positive_operator_gens(K)
357 sage: pi_cone = Cone([p.list() for p in pi_of_K],
360 sage: actual = pi_cone.dim()
364 sage: K.is_full_space()
366 sage: pi_of_K = positive_operator_gens(K)
367 sage: pi_cone = Cone([p.list() for p in pi_of_K],
370 sage: actual = pi_cone.dim()
373 sage: K = Cone([(1,0),(0,1),(0,-1)])
374 sage: pi_of_K = positive_operator_gens(K)
375 sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim()
379 The lineality of the cone of positive operators follows from the
380 description of its generators::
382 sage: set_random_seed()
383 sage: K = random_cone(max_ambient_dim=4)
384 sage: n = K.lattice_dim()
385 sage: pi_of_K = positive_operator_gens(K)
386 sage: L = ToricLattice(n**2)
387 sage: pi_cone = Cone([p.list() for p in pi_of_K],
390 sage: actual = pi_cone.lineality()
391 sage: expected = n**2 - K.dim()*K.dual().dim()
392 sage: actual == expected
395 The trivial cone, full space, and half-plane all give rise to the
396 expected linealities::
398 sage: n = ZZ.random_element().abs()
399 sage: K = Cone([[0] * n], ToricLattice(n))
402 sage: L = ToricLattice(n^2)
403 sage: pi_of_K = positive_operator_gens(K)
404 sage: pi_cone = Cone([p.list() for p in pi_of_K],
407 sage: actual = pi_cone.lineality()
411 sage: K.is_full_space()
413 sage: pi_of_K = positive_operator_gens(K)
414 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
417 sage: K = Cone([(1,0),(0,1),(0,-1)])
418 sage: pi_of_K = positive_operator_gens(K)
419 sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False)
420 sage: actual = pi_cone.lineality()
424 A cone is proper if and only if its cone of positive operators
427 sage: set_random_seed()
428 sage: K = random_cone(max_ambient_dim=4)
429 sage: pi_of_K = positive_operator_gens(K)
430 sage: L = ToricLattice(K.lattice_dim()**2)
431 sage: pi_cone = Cone([p.list() for p in pi_of_K],
434 sage: K.is_proper() == pi_cone.is_proper()
437 # Matrices are not vectors in Sage, so we have to convert them
438 # to vectors explicitly before we can find a basis. We need these
439 # two values to construct the appropriate "long vector" space.
440 F
= K
.lattice().base_field()
443 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
445 # Convert those tensor products to long vectors.
446 W
= VectorSpace(F
, n
**2)
447 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
450 if K
.is_solid() or K
.is_strictly_convex():
451 # The lineality space of either ``K`` or ``K.dual()`` is
452 # trivial and it's easy to show that our generating set is
453 # minimal. I would love a proof that this works when ``K`` is
454 # neither pointed nor solid.
456 # Note that in that case we can get *duplicates*, since the
457 # tensor product of (x,s) is the same as that of (-x,-s).
460 # Create the dual cone of the positive operators, expressed as
462 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()), check
=check
)
464 # Now compute the desired cone from its dual...
465 pi_cone
= pi_dual
.dual()
467 # And finally convert its rays back to matrix representations.
468 M
= MatrixSpace(F
, n
)
469 return [ M(v
.list()) for v
in pi_cone
]
472 def Z_transformation_gens(K
):
474 Compute generators of the cone of Z-transformations on this cone.
478 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
479 Each matrix ``L`` in the list should have the property that
480 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
481 discrete complementarity set of ``K``. Moreover, any nonnegative
482 linear combination of these matrices shares the same property.
486 Z-transformations on the nonnegative orthant are just Z-matrices.
487 That is, matrices whose off-diagonal elements are nonnegative::
489 sage: K = Cone([(1,0),(0,1)])
490 sage: Z_transformation_gens(K)
492 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
493 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
495 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
496 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
497 ....: for i in range(z.nrows())
498 ....: for j in range(z.ncols())
502 The trivial cone in a trivial space has no Z-transformations::
504 sage: K = Cone([], ToricLattice(0))
505 sage: Z_transformation_gens(K)
508 Z-transformations on a subspace are Lyapunov-like and vice-versa::
510 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
511 sage: K.is_full_space()
513 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
514 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
520 The Z-property is possessed by every Z-transformation::
522 sage: set_random_seed()
523 sage: K = random_cone(max_ambient_dim=4)
524 sage: Z_of_K = Z_transformation_gens(K)
525 sage: dcs = K.discrete_complementarity_set()
526 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
527 ....: for (x,s) in dcs])
530 The lineality space of Z is LL::
532 sage: set_random_seed()
533 sage: K = random_cone(max_ambient_dim=4)
534 sage: L = ToricLattice(K.lattice_dim()**2)
535 sage: z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ],
538 sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
539 sage: lls = L.vector_space().span(ll_basis)
540 sage: z_cone.linear_subspace() == lls
543 And thus, the lineality of Z is the Lyapunov rank::
545 sage: set_random_seed()
546 sage: K = random_cone(max_ambient_dim=4)
547 sage: Z_of_K = Z_transformation_gens(K)
548 sage: L = ToricLattice(K.lattice_dim()**2)
549 sage: z_cone = Cone([ z.list() for z in Z_of_K ],
552 sage: z_cone.lineality() == K.lyapunov_rank()
555 The lineality spaces of pi-star and Z-star are equal:
557 sage: set_random_seed()
558 sage: K = random_cone(max_ambient_dim=4)
559 sage: pi_of_K = positive_operator_gens(K)
560 sage: Z_of_K = Z_transformation_gens(K)
561 sage: L = ToricLattice(K.lattice_dim()**2)
562 sage: pi_cone = Cone([p.list() for p in pi_of_K],
565 sage: pi_star = pi_cone.dual()
566 sage: z_cone = Cone([ z.list() for z in Z_of_K],
569 sage: z_star = z_cone.dual()
570 sage: pi_star.linear_subspace() == z_star.linear_subspace()
573 # Matrices are not vectors in Sage, so we have to convert them
574 # to vectors explicitly before we can find a basis. We need these
575 # two values to construct the appropriate "long vector" space.
576 F
= K
.lattice().base_field()
579 # These tensor products contain generators for the dual cone of
580 # the cross-positive transformations.
581 tensor_products
= [ s
.tensor_product(x
)
582 for (x
,s
) in K
.discrete_complementarity_set() ]
584 # Turn our matrices into long vectors...
585 W
= VectorSpace(F
, n
**2)
586 vectors
= [ W(m
.list()) for m
in tensor_products
]
589 if K
.is_solid() or K
.is_strictly_convex():
590 # The lineality space of either ``K`` or ``K.dual()`` is
591 # trivial and it's easy to show that our generating set is
592 # minimal. I would love a proof that this works when ``K`` is
593 # neither pointed nor solid.
595 # Note that in that case we can get *duplicates*, since the
596 # tensor product of (x,s) is the same as that of (-x,-s).
599 # Create the dual cone of the cross-positive operators,
600 # expressed as long vectors.
601 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()), check
=check
)
603 # Now compute the desired cone from its dual...
604 Sigma_cone
= Sigma_dual
.dual()
606 # And finally convert its rays back to matrix representations.
607 # But first, make them negative, so we get Z-transformations and
608 # not cross-positive ones.
609 M
= MatrixSpace(F
, n
)
610 return [ -M(v
.list()) for v
in Sigma_cone
]
614 gens
= Z_transformation_gens(K
)
615 L
= ToricLattice(K
.lattice_dim()**2)
616 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)
619 gens
= positive_operator_gens(K
)
620 L
= ToricLattice(K
.lattice_dim()**2)
621 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)