]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
3b3a02814ea099a859f29ce86a40c7344c33b4cc
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.
13 There are faster ways of checking this property. For example, we
14 could compute a `lyapunov_like_basis` of the cone, and then test
15 whether or not the given matrix is contained in the span of that
16 basis. The value of this function is that it works on symbolic
21 - ``L`` -- A linear transformation or matrix.
23 - ``K`` -- A polyhedral closed convex cone.
27 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
28 and ``False`` otherwise.
32 If this function returns ``True``, then ``L`` is Lyapunov-like
33 on ``K``. However, if ``False`` is returned, that could mean one
34 of two things. The first is that ``L`` is definitely not
35 Lyapunov-like on ``K``. The second is more of an "I don't know"
36 answer, returned (for example) if we cannot prove that an inner
41 M. Orlitzky. The Lyapunov rank of an improper cone.
42 http://www.optimization-online.org/DB_HTML/2015/10/5135.html
46 The identity is always Lyapunov-like in a nontrivial space::
48 sage: set_random_seed()
49 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
50 sage: L = identity_matrix(K.lattice_dim())
51 sage: is_lyapunov_like(L,K)
54 As is the "zero" transformation::
56 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
57 sage: R = K.lattice().vector_space().base_ring()
58 sage: L = zero_matrix(R, K.lattice_dim())
59 sage: is_lyapunov_like(L,K)
62 Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
65 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
66 sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
70 return all([(L
*x
).inner_product(s
) == 0
71 for (x
,s
) in K
.discrete_complementarity_set()])
74 def positive_operator_gens(K
):
76 Compute generators of the cone of positive operators on this cone.
80 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
81 Each matrix ``P`` in the list should have the property that ``P*x``
82 is an element of ``K`` whenever ``x`` is an element of
83 ``K``. Moreover, any nonnegative linear combination of these
84 matrices shares the same property.
90 Positive and Z-operators on closed convex cones.
94 Some results of polyhedral cones and simplicial cones.
95 Linear and Multilinear Algebra, 4:4 (1977) 281--284.
99 Positive operators on the nonnegative orthant are nonnegative matrices::
101 sage: K = Cone([(1,)])
102 sage: positive_operator_gens(K)
105 sage: K = Cone([(1,0),(0,1)])
106 sage: positive_operator_gens(K)
108 [1 0] [0 1] [0 0] [0 0]
109 [0 0], [0 0], [1 0], [0 1]
112 The trivial cone in a trivial space has no positive operators::
114 sage: K = Cone([], ToricLattice(0))
115 sage: positive_operator_gens(K)
118 Every operator is positive on the trivial cone::
120 sage: K = Cone([(0,)])
121 sage: positive_operator_gens(K)
124 sage: K = Cone([(0,0)])
127 sage: positive_operator_gens(K)
129 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
130 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
133 Every operator is positive on the ambient vector space::
135 sage: K = Cone([(1,),(-1,)])
136 sage: K.is_full_space()
138 sage: positive_operator_gens(K)
141 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
142 sage: K.is_full_space()
144 sage: positive_operator_gens(K)
146 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
147 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
150 A non-obvious application is to find the positive operators on the
153 sage: K = Cone([(1,0),(0,1),(0,-1)])
154 sage: positive_operator_gens(K)
156 [1 0] [0 0] [ 0 0] [0 0] [ 0 0]
157 [0 0], [1 0], [-1 0], [0 1], [ 0 -1]
162 Each positive operator generator should send the generators of the
165 sage: set_random_seed()
166 sage: K = random_cone(max_ambient_dim=4)
167 sage: pi_of_K = positive_operator_gens(K)
168 sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
171 Each positive operator generator should send a random element of the
174 sage: set_random_seed()
175 sage: K = random_cone(max_ambient_dim=4)
176 sage: pi_of_K = positive_operator_gens(K)
177 sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
180 A random element of the positive operator cone should send the
181 generators of the cone into the cone::
183 sage: set_random_seed()
184 sage: K = random_cone(max_ambient_dim=4)
185 sage: pi_of_K = positive_operator_gens(K)
186 sage: L = ToricLattice(K.lattice_dim()**2)
187 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
190 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
191 sage: all([ K.contains(P*x) for x in K ])
194 A random element of the positive operator cone should send a random
195 element of the cone into the cone::
197 sage: set_random_seed()
198 sage: K = random_cone(max_ambient_dim=4)
199 sage: pi_of_K = positive_operator_gens(K)
200 sage: L = ToricLattice(K.lattice_dim()**2)
201 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
204 sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
205 sage: K.contains(P*K.random_element(ring=QQ))
208 The lineality space of the dual of the cone of positive operators
209 can be computed from the lineality spaces of the cone and its dual::
211 sage: set_random_seed()
212 sage: K = random_cone(max_ambient_dim=4)
213 sage: pi_of_K = positive_operator_gens(K)
214 sage: L = ToricLattice(K.lattice_dim()**2)
215 sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
218 sage: actual = pi_cone.dual().linear_subspace()
219 sage: U1 = [ vector((s.tensor_product(x)).list())
220 ....: for x in K.lines()
221 ....: for s in K.dual() ]
222 sage: U2 = [ vector((s.tensor_product(x)).list())
224 ....: for s in K.dual().lines() ]
225 sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
226 sage: actual == expected
229 The lineality of the dual of the cone of positive operators
230 is known from its lineality space::
232 sage: set_random_seed()
233 sage: K = random_cone(max_ambient_dim=4)
234 sage: n = K.lattice_dim()
236 sage: l = K.lineality()
237 sage: pi_of_K = positive_operator_gens(K)
238 sage: L = ToricLattice(n**2)
239 sage: pi_cone = Cone([p.list() for p in pi_of_K],
242 sage: actual = pi_cone.dual().lineality()
243 sage: expected = l*(m - l) + m*(n - m)
244 sage: actual == expected
247 The dimension of the cone of positive operators is given by the
248 corollary in my paper::
250 sage: set_random_seed()
251 sage: K = random_cone(max_ambient_dim=4)
252 sage: n = K.lattice_dim()
254 sage: l = K.lineality()
255 sage: pi_of_K = positive_operator_gens(K)
256 sage: L = ToricLattice(n**2)
257 sage: pi_cone = Cone([p.list() for p in pi_of_K],
260 sage: actual = pi_cone.dim()
261 sage: expected = n**2 - l*(m - l) - (n - m)*m
262 sage: actual == expected
265 The trivial cone, full space, and half-plane all give rise to the
266 expected dimensions::
268 sage: n = ZZ.random_element().abs()
269 sage: K = Cone([[0] * n], ToricLattice(n))
272 sage: L = ToricLattice(n^2)
273 sage: pi_of_K = positive_operator_gens(K)
274 sage: pi_cone = Cone([p.list() for p in pi_of_K],
277 sage: actual = pi_cone.dim()
281 sage: K.is_full_space()
283 sage: pi_of_K = positive_operator_gens(K)
284 sage: pi_cone = Cone([p.list() for p in pi_of_K],
287 sage: actual = pi_cone.dim()
290 sage: K = Cone([(1,0),(0,1),(0,-1)])
291 sage: pi_of_K = positive_operator_gens(K)
292 sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim()
296 The lineality of the cone of positive operators follows from the
297 description of its generators::
299 sage: set_random_seed()
300 sage: K = random_cone(max_ambient_dim=4)
301 sage: n = K.lattice_dim()
302 sage: pi_of_K = positive_operator_gens(K)
303 sage: L = ToricLattice(n**2)
304 sage: pi_cone = Cone([p.list() for p in pi_of_K],
307 sage: actual = pi_cone.lineality()
308 sage: expected = n**2 - K.dim()*K.dual().dim()
309 sage: actual == expected
312 The trivial cone, full space, and half-plane all give rise to the
313 expected linealities::
315 sage: n = ZZ.random_element().abs()
316 sage: K = Cone([[0] * n], ToricLattice(n))
319 sage: L = ToricLattice(n^2)
320 sage: pi_of_K = positive_operator_gens(K)
321 sage: pi_cone = Cone([p.list() for p in pi_of_K],
324 sage: actual = pi_cone.lineality()
328 sage: K.is_full_space()
330 sage: pi_of_K = positive_operator_gens(K)
331 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
332 sage: pi_cone.lineality() == n^2
334 sage: K = Cone([(1,0),(0,1),(0,-1)])
335 sage: pi_of_K = positive_operator_gens(K)
336 sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False)
337 sage: actual = pi_cone.lineality()
341 A cone is proper if and only if its cone of positive operators
344 sage: set_random_seed()
345 sage: K = random_cone(max_ambient_dim=4)
346 sage: pi_of_K = positive_operator_gens(K)
347 sage: L = ToricLattice(K.lattice_dim()**2)
348 sage: pi_cone = Cone([p.list() for p in pi_of_K],
351 sage: K.is_proper() == pi_cone.is_proper()
354 The positive operators of a permuted cone can be obtained by
357 sage: set_random_seed()
358 sage: K = random_cone(max_ambient_dim=4)
359 sage: L = ToricLattice(K.lattice_dim()**2)
360 sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
361 sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
362 sage: pi_of_pK = positive_operator_gens(pK)
363 sage: actual = Cone([t.list() for t in pi_of_pK],
366 sage: pi_of_K = positive_operator_gens(K)
367 sage: expected = Cone([(p*t*p.inverse()).list() for t in pi_of_K],
370 sage: actual.is_equivalent(expected)
373 A transformation is positive on a cone if and only if its adjoint is
374 positive on the dual of that cone::
376 sage: set_random_seed()
377 sage: K = random_cone(max_ambient_dim=4)
378 sage: F = K.lattice().vector_space().base_field()
379 sage: n = K.lattice_dim()
380 sage: L = ToricLattice(n**2)
381 sage: W = VectorSpace(F, n**2)
382 sage: pi_of_K = positive_operator_gens(K)
383 sage: pi_of_K_star = positive_operator_gens(K.dual())
384 sage: pi_cone = Cone([p.list() for p in pi_of_K],
387 sage: pi_star = Cone([p.list() for p in pi_of_K_star],
390 sage: M = MatrixSpace(F, n)
391 sage: L = M(pi_cone.random_element(ring=QQ).list())
392 sage: pi_star.contains(W(L.transpose().list()))
395 sage: L = W.random_element()
396 sage: L_star = W(M(L.list()).transpose().list())
397 sage: pi_cone.contains(L) == pi_star.contains(L_star)
400 # Matrices are not vectors in Sage, so we have to convert them
401 # to vectors explicitly before we can find a basis. We need these
402 # two values to construct the appropriate "long vector" space.
403 F
= K
.lattice().base_field()
406 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
408 # Convert those tensor products to long vectors.
409 W
= VectorSpace(F
, n
**2)
410 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
414 # All of the generators involved are extreme vectors and
415 # therefore minimal [Tam]_. If this cone is neither solid nor
416 # strictly convex, then the tensor product of ``s`` and ``x``
417 # is the same as that of ``-s`` and ``-x``. However, as a
418 # /set/, ``tensor_products`` may still be minimal.
421 # Create the dual cone of the positive operators, expressed as
423 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()), check
=check
)
425 # Now compute the desired cone from its dual...
426 pi_cone
= pi_dual
.dual()
428 # And finally convert its rays back to matrix representations.
429 M
= MatrixSpace(F
, n
)
430 return [ M(v
.list()) for v
in pi_cone
]
433 def Z_operator_gens(K
):
435 Compute generators of the cone of Z-operators on this cone.
439 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
440 Each matrix ``L`` in the list should have the property that
441 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element of
442 this cone's :meth:`discrete_complementarity_set`. Moreover, any
443 conic (nonnegative linear) combination of these matrices shares the
449 Positive and Z-operators on closed convex cones.
453 Z-operators on the nonnegative orthant are just Z-matrices.
454 That is, matrices whose off-diagonal elements are nonnegative::
456 sage: K = Cone([(1,0),(0,1)])
457 sage: Z_operator_gens(K)
459 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
460 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
462 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
463 sage: all([ z[i][j] <= 0 for z in Z_operator_gens(K)
464 ....: for i in range(z.nrows())
465 ....: for j in range(z.ncols())
469 The trivial cone in a trivial space has no Z-operators::
471 sage: K = Cone([], ToricLattice(0))
472 sage: Z_operator_gens(K)
475 Every operator is a Z-operator on the ambient vector space::
477 sage: K = Cone([(1,),(-1,)])
478 sage: K.is_full_space()
480 sage: Z_operator_gens(K)
483 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
484 sage: K.is_full_space()
486 sage: Z_operator_gens(K)
488 [-1 0] [1 0] [ 0 -1] [0 1] [ 0 0] [0 0] [ 0 0] [0 0]
489 [ 0 0], [0 0], [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
492 A non-obvious application is to find the Z-operators on the
495 sage: K = Cone([(1,0),(0,1),(0,-1)])
496 sage: Z_operator_gens(K)
498 [-1 0] [1 0] [ 0 0] [0 0] [ 0 0] [0 0]
499 [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1]
502 Z-operators on a subspace are Lyapunov-like and vice-versa::
504 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
505 sage: K.is_full_space()
507 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
508 sage: zs = span([ vector(z.list()) for z in Z_operator_gens(K) ])
514 The Z-property is possessed by every Z-operator::
516 sage: set_random_seed()
517 sage: K = random_cone(max_ambient_dim=4)
518 sage: Z_of_K = Z_operator_gens(K)
519 sage: dcs = K.discrete_complementarity_set()
520 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
521 ....: for (x,s) in dcs])
524 The lineality space of the cone of Z-operators is the space of
525 Lyapunov-like operators::
527 sage: set_random_seed()
528 sage: K = random_cone(max_ambient_dim=4)
529 sage: L = ToricLattice(K.lattice_dim()**2)
530 sage: Z_cone = Cone([ z.list() for z in Z_operator_gens(K) ],
533 sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
534 sage: lls = L.vector_space().span(ll_basis)
535 sage: Z_cone.linear_subspace() == lls
538 The lineality of the Z-operators on a cone is the Lyapunov
541 sage: set_random_seed()
542 sage: K = random_cone(max_ambient_dim=4)
543 sage: Z_of_K = Z_operator_gens(K)
544 sage: L = ToricLattice(K.lattice_dim()**2)
545 sage: Z_cone = Cone([ z.list() for z in Z_of_K ],
548 sage: Z_cone.lineality() == K.lyapunov_rank()
551 The lineality spaces of the duals of the positive and Z-operator
552 cones are equal. From this it follows that the dimensions of the
553 Z-operator cone and positive operator cone are equal::
555 sage: set_random_seed()
556 sage: K = random_cone(max_ambient_dim=4)
557 sage: pi_of_K = positive_operator_gens(K)
558 sage: Z_of_K = Z_operator_gens(K)
559 sage: L = ToricLattice(K.lattice_dim()**2)
560 sage: pi_cone = Cone([p.list() for p in pi_of_K],
563 sage: Z_cone = Cone([ z.list() for z in Z_of_K],
566 sage: pi_cone.dim() == Z_cone.dim()
568 sage: pi_star = pi_cone.dual()
569 sage: z_star = Z_cone.dual()
570 sage: pi_star.linear_subspace() == z_star.linear_subspace()
573 The trivial cone, full space, and half-plane all give rise to the
574 expected dimensions::
576 sage: n = ZZ.random_element().abs()
577 sage: K = Cone([[0] * n], ToricLattice(n))
580 sage: L = ToricLattice(n^2)
581 sage: Z_of_K = Z_operator_gens(K)
582 sage: Z_cone = Cone([z.list() for z in Z_of_K],
585 sage: actual = Z_cone.dim()
589 sage: K.is_full_space()
591 sage: Z_of_K = Z_operator_gens(K)
592 sage: Z_cone = Cone([z.list() for z in Z_of_K],
595 sage: actual = Z_cone.dim()
598 sage: K = Cone([(1,0),(0,1),(0,-1)])
599 sage: Z_of_K = Z_operator_gens(K)
600 sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
601 sage: Z_cone.dim() == 3
604 The Z-operators of a permuted cone can be obtained by conjugation::
606 sage: set_random_seed()
607 sage: K = random_cone(max_ambient_dim=4)
608 sage: L = ToricLattice(K.lattice_dim()**2)
609 sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
610 sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
611 sage: Z_of_pK = Z_operator_gens(pK)
612 sage: actual = Cone([t.list() for t in Z_of_pK],
615 sage: Z_of_K = Z_operator_gens(K)
616 sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K],
619 sage: actual.is_equivalent(expected)
622 An operator is a Z-operator on a cone if and only if its
623 adjoint is a Z-operator on the dual of that cone::
625 sage: set_random_seed()
626 sage: K = random_cone(max_ambient_dim=4)
627 sage: F = K.lattice().vector_space().base_field()
628 sage: n = K.lattice_dim()
629 sage: L = ToricLattice(n**2)
630 sage: W = VectorSpace(F, n**2)
631 sage: Z_of_K = Z_operator_gens(K)
632 sage: Z_of_K_star = Z_operator_gens(K.dual())
633 sage: Z_cone = Cone([p.list() for p in Z_of_K],
636 sage: Z_star = Cone([p.list() for p in Z_of_K_star],
639 sage: M = MatrixSpace(F, n)
640 sage: L = M(Z_cone.random_element(ring=QQ).list())
641 sage: Z_star.contains(W(L.transpose().list()))
644 sage: L = W.random_element()
645 sage: L_star = W(M(L.list()).transpose().list())
646 sage: Z_cone.contains(L) == Z_star.contains(L_star)
649 # Matrices are not vectors in Sage, so we have to convert them
650 # to vectors explicitly before we can find a basis. We need these
651 # two values to construct the appropriate "long vector" space.
652 F
= K
.lattice().base_field()
655 # These tensor products contain generators for the dual cone of
656 # the cross-positive operators.
657 tensor_products
= [ s
.tensor_product(x
)
658 for (x
,s
) in K
.discrete_complementarity_set() ]
660 # Turn our matrices into long vectors...
661 W
= VectorSpace(F
, n
**2)
662 vectors
= [ W(m
.list()) for m
in tensor_products
]
666 # All of the generators involved are extreme vectors and
667 # therefore minimal. If this cone is neither solid nor
668 # strictly convex, then the tensor product of ``s`` and ``x``
669 # is the same as that of ``-s`` and ``-x``. However, as a
670 # /set/, ``tensor_products`` may still be minimal.
673 # Create the dual cone of the cross-positive operators,
674 # expressed as long vectors.
675 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()), check
=check
)
677 # Now compute the desired cone from its dual...
678 Sigma_cone
= Sigma_dual
.dual()
680 # And finally convert its rays back to matrix representations.
681 # But first, make them negative, so we get Z-operators and
682 # not cross-positive ones.
683 M
= MatrixSpace(F
, n
)
684 return [ -M(v
.list()) for v
in Sigma_cone
]
688 gens
= Z_operator_gens(K
)
689 L
= ToricLattice(K
.lattice_dim()**2)
690 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)
693 gens
= positive_operator_gens(K
)
694 L
= ToricLattice(K
.lattice_dim()**2)
695 return Cone([ g
.list() for g
in gens
], lattice
=L
, check
=False)