]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
66d6ccfbf1958f9df3e0452ea527a436890293ad
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 A strictly convex cone should be equal to its strictly convex component::
140 sage: set_random_seed()
141 sage: K = random_cone(max_ambient_dim=8, strictly_convex=True)
142 sage: (P,_) = motzkin_decomposition(K)
143 sage: K.is_equivalent(P)
146 The generators of the components are obtained from orthogonal
147 projections of the original generators [Stoer-Witzgall]_::
149 sage: set_random_seed()
150 sage: K = random_cone(max_ambient_dim=8)
151 sage: (P,S) = motzkin_decomposition(K)
152 sage: A = S.linear_subspace().complement().matrix()
153 sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A
154 sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice())
155 sage: P.is_equivalent(expected_P)
157 sage: A = S.linear_subspace().matrix()
158 sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A
159 sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice())
160 sage: S.is_equivalent(expected_S)
163 # The lines() method only returns one generator per line. For a true
164 # line, we also need a generator pointing in the opposite direction.
165 S_gens
= [ direction
*gen
for direction
in [1,-1] for gen
in K
.lines() ]
166 S
= Cone(S_gens
, K
.lattice(), check
=False)
168 # Since ``S`` is a subspace, the rays of its dual generate its
169 # orthogonal complement.
170 S_perp
= Cone(S
.dual(), K
.lattice(), check
=False)
171 P
= K
.intersection(S_perp
)
176 def positive_operator_gens(K
):
178 Compute generators of the cone of positive operators on this cone.
182 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
183 Each matrix ``P`` in the list should have the property that ``P*x``
184 is an element of ``K`` whenever ``x`` is an element of
185 ``K``. Moreover, any nonnegative linear combination of these
186 matrices shares the same property.
190 Positive operators on the nonnegative orthant are nonnegative matrices::
192 sage: K = Cone([(1,)])
193 sage: positive_operator_gens(K)
196 sage: K = Cone([(1,0),(0,1)])
197 sage: positive_operator_gens(K)
199 [1 0] [0 1] [0 0] [0 0]
200 [0 0], [0 0], [1 0], [0 1]
203 The trivial cone in a trivial space has no positive operators::
205 sage: K = Cone([], ToricLattice(0))
206 sage: positive_operator_gens(K)
209 Every operator is positive on the trivial cone::
211 sage: K = Cone([(0,)])
212 sage: positive_operator_gens(K)
215 sage: K = Cone([(0,0)])
218 sage: positive_operator_gens(K)
220 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
221 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
224 Every operator is positive on the ambient vector space::
226 sage: K = Cone([(1,),(-1,)])
227 sage: K.is_full_space()
229 sage: positive_operator_gens(K)
232 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
233 sage: K.is_full_space()
235 sage: positive_operator_gens(K)
237 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
238 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
241 A non-obvious application is to find the positive operators on the
244 sage: K = Cone([(1,0),(0,1),(0,-1)])
245 sage: positive_operator_gens(K)
247 [1 0] [0 0] [ 0 0] [0 0] [ 0 0]
248 [0 0], [1 0], [-1 0], [0 1], [ 0 -1]
253 Each positive operator generator should send the generators of the
256 sage: set_random_seed()
257 sage: K = random_cone(max_ambient_dim=4)
258 sage: pi_of_K = positive_operator_gens(K)
259 sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
262 Each positive operator generator should send a random element of the
265 sage: set_random_seed()
266 sage: K = random_cone(max_ambient_dim=4)
267 sage: pi_of_K = positive_operator_gens(K)
268 sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
271 A random element of the positive operator cone should send the
272 generators of the cone into the cone::
274 sage: set_random_seed()
275 sage: K = random_cone(max_ambient_dim=4)
276 sage: pi_of_K = positive_operator_gens(K)
277 sage: L = ToricLattice(K.lattice_dim()**2)
278 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
281 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
282 sage: all([ K.contains(P*x) for x in K ])
285 A random element of the positive operator cone should send a random
286 element of the cone into the cone::
288 sage: set_random_seed()
289 sage: K = random_cone(max_ambient_dim=4)
290 sage: pi_of_K = positive_operator_gens(K)
291 sage: L = ToricLattice(K.lattice_dim()**2)
292 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
295 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
296 sage: K.contains(P*K.random_element(ring=QQ))
299 The lineality space of the dual of the cone of positive operators
300 can be computed from the lineality spaces of the cone and its dual::
302 sage: set_random_seed()
303 sage: K = random_cone(max_ambient_dim=4)
304 sage: pi_of_K = positive_operator_gens(K)
305 sage: L = ToricLattice(K.lattice_dim()**2)
306 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
309 sage: actual = pi_cone.dual().linear_subspace()
310 sage: U1 = [ vector((s.tensor_product(x)).list())
311 ....: for x in K.lines()
312 ....: for s in K.dual() ]
313 sage: U2 = [ vector((s.tensor_product(x)).list())
315 ....: for s in K.dual().lines() ]
316 sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
317 sage: actual == expected
320 The lineality of the dual of the cone of positive operators
321 is known from its lineality space::
323 sage: set_random_seed()
324 sage: K = random_cone(max_ambient_dim=4)
325 sage: n = K.lattice_dim()
327 sage: l = K.lineality()
328 sage: pi_of_K = positive_operator_gens(K)
329 sage: L = ToricLattice(n**2)
330 sage: pi_cone = Cone([p.list() for p in pi_of_K],
333 sage: actual = pi_cone.dual().lineality()
334 sage: expected = l*(m - l) + m*(n - m)
335 sage: actual == expected
338 The dimension of the cone of positive operators is given by the
339 corollary in my paper::
341 sage: set_random_seed()
342 sage: K = random_cone(max_ambient_dim=4)
343 sage: n = K.lattice_dim()
345 sage: l = K.lineality()
346 sage: pi_of_K = positive_operator_gens(K)
347 sage: L = ToricLattice(n**2)
348 sage: pi_cone = Cone([p.list() for p in pi_of_K],
351 sage: actual = pi_cone.dim()
352 sage: expected = n**2 - l*(m - l) - (n - m)*m
353 sage: actual == expected
356 The trivial cone, full space, and half-plane all give rise to the
357 expected dimensions::
359 sage: n = ZZ.random_element().abs()
360 sage: K = Cone([[0] * n], ToricLattice(n))
363 sage: L = ToricLattice(n^2)
364 sage: pi_of_K = positive_operator_gens(K)
365 sage: pi_cone = Cone([p.list() for p in pi_of_K],
368 sage: actual = pi_cone.dim()
372 sage: K.is_full_space()
374 sage: pi_of_K = positive_operator_gens(K)
375 sage: pi_cone = Cone([p.list() for p in pi_of_K],
378 sage: actual = pi_cone.dim()
381 sage: K = Cone([(1,0),(0,1),(0,-1)])
382 sage: pi_of_K = positive_operator_gens(K)
383 sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim()
387 The lineality of the cone of positive operators follows from the
388 description of its generators::
390 sage: set_random_seed()
391 sage: K = random_cone(max_ambient_dim=4)
392 sage: n = K.lattice_dim()
393 sage: pi_of_K = positive_operator_gens(K)
394 sage: L = ToricLattice(n**2)
395 sage: pi_cone = Cone([p.list() for p in pi_of_K],
398 sage: actual = pi_cone.lineality()
399 sage: expected = n**2 - K.dim()*K.dual().dim()
400 sage: actual == expected
403 The trivial cone, full space, and half-plane all give rise to the
404 expected linealities::
406 sage: n = ZZ.random_element().abs()
407 sage: K = Cone([[0] * n], ToricLattice(n))
410 sage: L = ToricLattice(n^2)
411 sage: pi_of_K = positive_operator_gens(K)
412 sage: pi_cone = Cone([p.list() for p in pi_of_K],
415 sage: actual = pi_cone.lineality()
419 sage: K.is_full_space()
421 sage: pi_of_K = positive_operator_gens(K)
422 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
423 sage: pi_cone.lineality() == n^2
425 sage: K = Cone([(1,0),(0,1),(0,-1)])
426 sage: pi_of_K = positive_operator_gens(K)
427 sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False)
428 sage: actual = pi_cone.lineality()
432 A cone is proper if and only if its cone of positive operators
435 sage: set_random_seed()
436 sage: K = random_cone(max_ambient_dim=4)
437 sage: pi_of_K = positive_operator_gens(K)
438 sage: L = ToricLattice(K.lattice_dim()**2)
439 sage: pi_cone = Cone([p.list() for p in pi_of_K],
442 sage: K.is_proper() == pi_cone.is_proper()
445 # Matrices are not vectors in Sage, so we have to convert them
446 # to vectors explicitly before we can find a basis. We need these
447 # two values to construct the appropriate "long vector" space.
448 F
= K
.lattice().base_field()
451 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
453 # Convert those tensor products to long vectors.
454 W
= VectorSpace(F
, n
**2)
455 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
458 if K
.is_solid() or K
.is_strictly_convex():
459 # The lineality space of either ``K`` or ``K.dual()`` is
460 # trivial and it's easy to show that our generating set is
461 # minimal. I would love a proof that this works when ``K`` is
462 # neither pointed nor solid.
464 # Note that in that case we can get *duplicates*, since the
465 # tensor product of (x,s) is the same as that of (-x,-s).
468 # Create the dual cone of the positive operators, expressed as
470 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()), check
=check
)
472 # Now compute the desired cone from its dual...
473 pi_cone
= pi_dual
.dual()
475 # And finally convert its rays back to matrix representations.
476 M
= MatrixSpace(F
, n
)
477 return [ M(v
.list()) for v
in pi_cone
]
480 def Z_transformation_gens(K
):
482 Compute generators of the cone of Z-transformations on this cone.
486 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
487 Each matrix ``L`` in the list should have the property that
488 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
489 discrete complementarity set of ``K``. Moreover, any nonnegative
490 linear combination of these matrices shares the same property.
494 Z-transformations on the nonnegative orthant are just Z-matrices.
495 That is, matrices whose off-diagonal elements are nonnegative::
497 sage: K = Cone([(1,0),(0,1)])
498 sage: Z_transformation_gens(K)
500 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
501 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
503 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
504 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
505 ....: for i in range(z.nrows())
506 ....: for j in range(z.ncols())
510 The trivial cone in a trivial space has no Z-transformations::
512 sage: K = Cone([], ToricLattice(0))
513 sage: Z_transformation_gens(K)
516 Every operator is a Z-transformation on the ambient vector space::
518 sage: K = Cone([(1,),(-1,)])
519 sage: K.is_full_space()
521 sage: Z_transformation_gens(K)
524 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
525 sage: K.is_full_space()
527 sage: Z_transformation_gens(K)
529 [-1 0] [1 0] [ 0 -1] [0 1] [ 0 0] [0 0] [ 0 0] [0 0]
530 [ 0 0], [0 0], [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
533 A non-obvious application is to find the Z-transformations on the
536 sage: K = Cone([(1,0),(0,1),(0,-1)])
537 sage: Z_transformation_gens(K)
539 [-1 0] [1 0] [ 0 0] [0 0] [ 0 0] [0 0]
540 [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
543 Z-transformations on a subspace are Lyapunov-like and vice-versa::
545 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
546 sage: K.is_full_space()
548 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
549 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
555 The Z-property is possessed by every Z-transformation::
557 sage: set_random_seed()
558 sage: K = random_cone(max_ambient_dim=4)
559 sage: Z_of_K = Z_transformation_gens(K)
560 sage: dcs = K.discrete_complementarity_set()
561 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
562 ....: for (x,s) in dcs])
565 The lineality space of the cone of Z-transformations is the space of
566 Lyapunov-like transformations::
568 sage: set_random_seed()
569 sage: K = random_cone(max_ambient_dim=4)
570 sage: L = ToricLattice(K.lattice_dim()**2)
571 sage: Z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ],
574 sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
575 sage: lls = L.vector_space().span(ll_basis)
576 sage: Z_cone.linear_subspace() == lls
579 The lineality of the Z-transformations on a cone is the Lyapunov
582 sage: set_random_seed()
583 sage: K = random_cone(max_ambient_dim=4)
584 sage: Z_of_K = Z_transformation_gens(K)
585 sage: L = ToricLattice(K.lattice_dim()**2)
586 sage: Z_cone = Cone([ z.list() for z in Z_of_K ],
589 sage: Z_cone.lineality() == K.lyapunov_rank()
592 The lineality spaces of the duals of the positive operator and
593 Z-transformation cones are equal. From this it follows that the
594 dimensions of the Z-transformation cone and positive operator cone
597 sage: set_random_seed()
598 sage: K = random_cone(max_ambient_dim=4)
599 sage: pi_of_K = positive_operator_gens(K)
600 sage: Z_of_K = Z_transformation_gens(K)
601 sage: L = ToricLattice(K.lattice_dim()**2)
602 sage: pi_cone = Cone([p.list() for p in pi_of_K],
605 sage: Z_cone = Cone([ z.list() for z in Z_of_K],
608 sage: pi_cone.dim() == Z_cone.dim()
610 sage: pi_star = pi_cone.dual()
611 sage: z_star = Z_cone.dual()
612 sage: pi_star.linear_subspace() == z_star.linear_subspace()
615 The trivial cone, full space, and half-plane all give rise to the
616 expected dimensions::
618 sage: n = ZZ.random_element().abs()
619 sage: K = Cone([[0] * n], ToricLattice(n))
622 sage: L = ToricLattice(n^2)
623 sage: Z_of_K = Z_transformation_gens(K)
624 sage: Z_cone = Cone([z.list() for z in Z_of_K],
627 sage: actual = Z_cone.dim()
631 sage: K.is_full_space()
633 sage: Z_of_K = Z_transformation_gens(K)
634 sage: Z_cone = Cone([z.list() for z in Z_of_K],
637 sage: actual = Z_cone.dim()
640 sage: K = Cone([(1,0),(0,1),(0,-1)])
641 sage: Z_of_K = Z_transformation_gens(K)
642 sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
643 sage: Z_cone.dim() == 3
646 # Matrices are not vectors in Sage, so we have to convert them
647 # to vectors explicitly before we can find a basis. We need these
648 # two values to construct the appropriate "long vector" space.
649 F
= K
.lattice().base_field()
652 # These tensor products contain generators for the dual cone of
653 # the cross-positive transformations.
654 tensor_products
= [ s
.tensor_product(x
)
655 for (x
,s
) in K
.discrete_complementarity_set() ]
657 # Turn our matrices into long vectors...
658 W
= VectorSpace(F
, n
**2)
659 vectors
= [ W(m
.list()) for m
in tensor_products
]
662 if K
.is_solid() or K
.is_strictly_convex():
663 # The lineality space of either ``K`` or ``K.dual()`` is
664 # trivial and it's easy to show that our generating set is
665 # minimal. I would love a proof that this works when ``K`` is
666 # neither pointed nor solid.
668 # Note that in that case we can get *duplicates*, since the
669 # tensor product of (x,s) is the same as that of (-x,-s).
672 # Create the dual cone of the cross-positive operators,
673 # expressed as long vectors.
674 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()), check
=check
)
676 # Now compute the desired cone from its dual...
677 Sigma_cone
= Sigma_dual
.dual()
679 # And finally convert its rays back to matrix representations.
680 # But first, make them negative, so we get Z-transformations and
681 # not cross-positive ones.
682 M
= MatrixSpace(F
, n
)
683 return [ -M(v
.list()) for v
in Sigma_cone
]
687 gens
= Z_transformation_gens(K
)
688 L
= ToricLattice(K
.lattice_dim()**2)
689 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)
692 gens
= positive_operator_gens(K
)
693 L
= ToricLattice(K
.lattice_dim()**2)
694 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)