]>
gitweb.michael.orlitzky.com - sage.d.git/blob - cone/cone.py
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 The positive operators of a permuted cone can be obtained by
448 sage: set_random_seed()
449 sage: K = random_cone(max_ambient_dim=4)
450 sage: L = ToricLattice(K.lattice_dim()**2)
451 sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
452 sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
453 sage: pi_of_pK = positive_operator_gens(pK)
454 sage: actual = Cone([t.list() for t in pi_of_pK],
457 sage: pi_of_K = positive_operator_gens(K)
458 sage: expected = Cone([(p*t*p.inverse()).list() for t in pi_of_K],
461 sage: actual.is_equivalent(expected)
464 # Matrices are not vectors in Sage, so we have to convert them
465 # to vectors explicitly before we can find a basis. We need these
466 # two values to construct the appropriate "long vector" space.
467 F
= K
.lattice().base_field()
470 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
472 # Convert those tensor products to long vectors.
473 W
= VectorSpace(F
, n
**2)
474 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
477 if K
.is_solid() or K
.is_strictly_convex():
478 # The lineality space of either ``K`` or ``K.dual()`` is
479 # trivial and it's easy to show that our generating set is
480 # minimal. I would love a proof that this works when ``K`` is
481 # neither pointed nor solid.
483 # Note that in that case we can get *duplicates*, since the
484 # tensor product of (x,s) is the same as that of (-x,-s).
487 # Create the dual cone of the positive operators, expressed as
489 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()), check
=check
)
491 # Now compute the desired cone from its dual...
492 pi_cone
= pi_dual
.dual()
494 # And finally convert its rays back to matrix representations.
495 M
= MatrixSpace(F
, n
)
496 return [ M(v
.list()) for v
in pi_cone
]
499 def Z_transformation_gens(K
):
501 Compute generators of the cone of Z-transformations on this cone.
505 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
506 Each matrix ``L`` in the list should have the property that
507 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
508 discrete complementarity set of ``K``. Moreover, any nonnegative
509 linear combination of these matrices shares the same property.
513 Z-transformations on the nonnegative orthant are just Z-matrices.
514 That is, matrices whose off-diagonal elements are nonnegative::
516 sage: K = Cone([(1,0),(0,1)])
517 sage: Z_transformation_gens(K)
519 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
520 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
522 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
523 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
524 ....: for i in range(z.nrows())
525 ....: for j in range(z.ncols())
529 The trivial cone in a trivial space has no Z-transformations::
531 sage: K = Cone([], ToricLattice(0))
532 sage: Z_transformation_gens(K)
535 Every operator is a Z-transformation on the ambient vector space::
537 sage: K = Cone([(1,),(-1,)])
538 sage: K.is_full_space()
540 sage: Z_transformation_gens(K)
543 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
544 sage: K.is_full_space()
546 sage: Z_transformation_gens(K)
548 [-1 0] [1 0] [ 0 -1] [0 1] [ 0 0] [0 0] [ 0 0] [0 0]
549 [ 0 0], [0 0], [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
552 A non-obvious application is to find the Z-transformations on the
555 sage: K = Cone([(1,0),(0,1),(0,-1)])
556 sage: Z_transformation_gens(K)
558 [-1 0] [1 0] [ 0 0] [0 0] [ 0 0] [0 0]
559 [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
562 Z-transformations on a subspace are Lyapunov-like and vice-versa::
564 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
565 sage: K.is_full_space()
567 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
568 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
574 The Z-property is possessed by every Z-transformation::
576 sage: set_random_seed()
577 sage: K = random_cone(max_ambient_dim=4)
578 sage: Z_of_K = Z_transformation_gens(K)
579 sage: dcs = K.discrete_complementarity_set()
580 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
581 ....: for (x,s) in dcs])
584 The lineality space of the cone of Z-transformations is the space of
585 Lyapunov-like transformations::
587 sage: set_random_seed()
588 sage: K = random_cone(max_ambient_dim=4)
589 sage: L = ToricLattice(K.lattice_dim()**2)
590 sage: Z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ],
593 sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
594 sage: lls = L.vector_space().span(ll_basis)
595 sage: Z_cone.linear_subspace() == lls
598 The lineality of the Z-transformations on a cone is the Lyapunov
601 sage: set_random_seed()
602 sage: K = random_cone(max_ambient_dim=4)
603 sage: Z_of_K = Z_transformation_gens(K)
604 sage: L = ToricLattice(K.lattice_dim()**2)
605 sage: Z_cone = Cone([ z.list() for z in Z_of_K ],
608 sage: Z_cone.lineality() == K.lyapunov_rank()
611 The lineality spaces of the duals of the positive operator and
612 Z-transformation cones are equal. From this it follows that the
613 dimensions of the Z-transformation cone and positive operator cone
616 sage: set_random_seed()
617 sage: K = random_cone(max_ambient_dim=4)
618 sage: pi_of_K = positive_operator_gens(K)
619 sage: Z_of_K = Z_transformation_gens(K)
620 sage: L = ToricLattice(K.lattice_dim()**2)
621 sage: pi_cone = Cone([p.list() for p in pi_of_K],
624 sage: Z_cone = Cone([ z.list() for z in Z_of_K],
627 sage: pi_cone.dim() == Z_cone.dim()
629 sage: pi_star = pi_cone.dual()
630 sage: z_star = Z_cone.dual()
631 sage: pi_star.linear_subspace() == z_star.linear_subspace()
634 The trivial cone, full space, and half-plane all give rise to the
635 expected dimensions::
637 sage: n = ZZ.random_element().abs()
638 sage: K = Cone([[0] * n], ToricLattice(n))
641 sage: L = ToricLattice(n^2)
642 sage: Z_of_K = Z_transformation_gens(K)
643 sage: Z_cone = Cone([z.list() for z in Z_of_K],
646 sage: actual = Z_cone.dim()
650 sage: K.is_full_space()
652 sage: Z_of_K = Z_transformation_gens(K)
653 sage: Z_cone = Cone([z.list() for z in Z_of_K],
656 sage: actual = Z_cone.dim()
659 sage: K = Cone([(1,0),(0,1),(0,-1)])
660 sage: Z_of_K = Z_transformation_gens(K)
661 sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
662 sage: Z_cone.dim() == 3
665 The Z-transformations of a permuted cone can be obtained by
668 sage: set_random_seed()
669 sage: K = random_cone(max_ambient_dim=4)
670 sage: L = ToricLattice(K.lattice_dim()**2)
671 sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
672 sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
673 sage: Z_of_pK = Z_transformation_gens(pK)
674 sage: actual = Cone([t.list() for t in Z_of_pK],
677 sage: Z_of_K = Z_transformation_gens(K)
678 sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K],
681 sage: actual.is_equivalent(expected)
684 # Matrices are not vectors in Sage, so we have to convert them
685 # to vectors explicitly before we can find a basis. We need these
686 # two values to construct the appropriate "long vector" space.
687 F
= K
.lattice().base_field()
690 # These tensor products contain generators for the dual cone of
691 # the cross-positive transformations.
692 tensor_products
= [ s
.tensor_product(x
)
693 for (x
,s
) in K
.discrete_complementarity_set() ]
695 # Turn our matrices into long vectors...
696 W
= VectorSpace(F
, n
**2)
697 vectors
= [ W(m
.list()) for m
in tensor_products
]
700 if K
.is_solid() or K
.is_strictly_convex():
701 # The lineality space of either ``K`` or ``K.dual()`` is
702 # trivial and it's easy to show that our generating set is
703 # minimal. I would love a proof that this works when ``K`` is
704 # neither pointed nor solid.
706 # Note that in that case we can get *duplicates*, since the
707 # tensor product of (x,s) is the same as that of (-x,-s).
710 # Create the dual cone of the cross-positive operators,
711 # expressed as long vectors.
712 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()), check
=check
)
714 # Now compute the desired cone from its dual...
715 Sigma_cone
= Sigma_dual
.dual()
717 # And finally convert its rays back to matrix representations.
718 # But first, make them negative, so we get Z-transformations and
719 # not cross-positive ones.
720 M
= MatrixSpace(F
, n
)
721 return [ -M(v
.list()) for v
in Sigma_cone
]
725 gens
= Z_transformation_gens(K
)
726 L
= ToricLattice(K
.lattice_dim()**2)
727 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)
730 gens
= positive_operator_gens(K
)
731 L
= ToricLattice(K
.lattice_dim()**2)
732 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)