]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
f8b879131908a8d1a15d5069fcc888c5db319b07
1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
2 # have to explicitly mangle our sitedir here so that "mjo.cone"
4 from os
.path
import abspath
5 from site
import addsitedir
6 addsitedir(abspath('../../'))
11 def _basically_the_same(K1
, K2
):
13 Test whether or not ``K1`` and ``K2`` are "basically the same."
15 This is a hack to get around the fact that it's difficult to tell
16 when two cones are linearly isomorphic. We have a proposition that
17 equates two cones, but represented over `\mathbb{Q}`, they are
18 merely linearly isomorphic (not equal). So rather than test for
19 equality, we test a list of properties that should be preserved
20 under an invertible linear transformation.
24 ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
29 Any proper cone with three generators in `\mathbb{R}^{3}` is
30 basically the same as the nonnegative orthant::
32 sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
33 sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
34 sage: _basically_the_same(K1, K2)
37 Negating a cone gives you another cone that is basically the same::
39 sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
40 sage: _basically_the_same(K, -K)
45 Any cone is basically the same as itself::
47 sage: K = random_cone(max_ambient_dim = 8)
48 sage: _basically_the_same(K, K)
51 After applying an invertible matrix to the rows of a cone, the
52 result should be basically the same as the cone we started with::
54 sage: K1 = random_cone(max_ambient_dim = 8)
55 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
56 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
57 sage: _basically_the_same(K1, K2)
61 if K1
.lattice_dim() != K2
.lattice_dim():
64 if K1
.nrays() != K2
.nrays():
67 if K1
.dim() != K2
.dim():
70 if K1
.lineality() != K2
.lineality():
73 if K1
.is_solid() != K2
.is_solid():
76 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
79 if len(LL(K1
)) != len(LL(K2
)):
82 C_of_K1
= K1
.discrete_complementarity_set()
83 C_of_K2
= K2
.discrete_complementarity_set()
84 if len(C_of_K1
) != len(C_of_K2
):
87 if len(K1
.facets()) != len(K2
.facets()):
94 def _restrict_to_space(K
, W
):
96 Restrict this cone a subspace of its ambient space.
100 - ``W`` -- The subspace into which this cone will be restricted.
104 A new cone in a sublattice corresponding to ``W``.
108 When this cone is solid, restricting it into its own span should do
111 sage: K = Cone([(1,)])
112 sage: _restrict_to_space(K, K.span()) == K
115 A single ray restricted into its own span gives the same output
116 regardless of the ambient space::
118 sage: K2 = Cone([(1,0)])
119 sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
123 sage: K3 = Cone([(1,0,0)])
124 sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
133 The projected cone should always be solid::
135 sage: set_random_seed()
136 sage: K = random_cone(max_ambient_dim = 8)
137 sage: _restrict_to_space(K, K.span()).is_solid()
140 And the resulting cone should live in a space having the same
141 dimension as the space we restricted it to::
143 sage: set_random_seed()
144 sage: K = random_cone(max_ambient_dim = 8)
145 sage: K_P = _restrict_to_space(K, K.dual().span())
146 sage: K_P.lattice_dim() == K.dual().dim()
149 This function should not affect the dimension of a cone::
151 sage: set_random_seed()
152 sage: K = random_cone(max_ambient_dim = 8)
153 sage: K.dim() == _restrict_to_space(K,K.span()).dim()
156 Nor should it affect the lineality of a cone::
158 sage: set_random_seed()
159 sage: K = random_cone(max_ambient_dim = 8)
160 sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
163 No matter which space we restrict to, the lineality should not
166 sage: set_random_seed()
167 sage: K = random_cone(max_ambient_dim = 8)
168 sage: S = K.span(); P = K.dual().span()
169 sage: K.lineality() >= _restrict_to_space(K,S).lineality()
171 sage: K.lineality() >= _restrict_to_space(K,P).lineality()
174 If we do this according to our paper, then the result is proper::
176 sage: set_random_seed()
177 sage: K = random_cone(max_ambient_dim = 8)
178 sage: K_S = _restrict_to_space(K, K.span())
179 sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
180 sage: K_SP.is_proper()
182 sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
183 sage: K_SP.is_proper()
186 Test the proposition in our paper concerning the duals and
187 restrictions. Generate a random cone, then create a subcone of
188 it. The operation of dual-taking should then commute with
191 sage: set_random_seed()
192 sage: J = random_cone(max_ambient_dim = 8)
193 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
194 sage: K_W_star = _restrict_to_space(K, J.span()).dual()
195 sage: K_star_W = _restrict_to_space(K.dual(), J.span())
196 sage: _basically_the_same(K_W_star, K_star_W)
200 # First we want to intersect ``K`` with ``W``. The easiest way to
201 # do this is via cone intersection, so we turn the subspace ``W``
203 W_cone
= Cone(W
.basis() + [-b
for b
in W
.basis()], lattice
=K
.lattice())
204 K
= K
.intersection(W_cone
)
206 # We've already intersected K with the span of K2, so every
207 # generator of K should belong to W now.
208 K_W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
210 L
= ToricLattice(W
.dimension())
211 return Cone(K_W_rays
, lattice
=L
)
216 Compute a basis of Lyapunov-like transformations on this cone.
220 A list of matrices forming a basis for the space of all
221 Lyapunov-like transformations on the given cone.
225 The trivial cone has no Lyapunov-like transformations::
227 sage: L = ToricLattice(0)
228 sage: K = Cone([], lattice=L)
232 The Lyapunov-like transformations on the nonnegative orthant are
233 simply diagonal matrices::
235 sage: K = Cone([(1,)])
239 sage: K = Cone([(1,0),(0,1)])
246 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
249 [1 0 0] [0 0 0] [0 0 0]
250 [0 0 0] [0 1 0] [0 0 0]
251 [0 0 0], [0 0 0], [0 0 1]
254 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
255 `L^{3}_{\infty}` cones [Rudolf et al.]_::
257 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
265 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
273 If our cone is the entire space, then every transformation on it is
276 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
277 sage: M = MatrixSpace(QQ,2)
278 sage: M.basis() == LL(K)
283 The inner product `\left< L\left(x\right), s \right>` is zero for
284 every pair `\left( x,s \right)` in the discrete complementarity set
287 sage: set_random_seed()
288 sage: K = random_cone(max_ambient_dim=8)
289 sage: C_of_K = K.discrete_complementarity_set()
290 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
291 sage: sum(map(abs, l))
294 The Lyapunov-like transformations on a cone and its dual are related
295 by transposition, but we're not guaranteed to compute transposed
296 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
299 sage: set_random_seed()
300 sage: K = random_cone(max_ambient_dim=8)
301 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
302 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
303 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
304 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
305 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
309 V
= K
.lattice().vector_space()
311 C_of_K
= K
.discrete_complementarity_set()
313 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
315 # Sage doesn't think matrices are vectors, so we have to convert
316 # our matrices to vectors explicitly before we can figure out how
317 # many are linearly-indepenedent.
319 # The space W has the same base ring as V, but dimension
320 # dim(V)^2. So it has the same dimension as the space of linear
321 # transformations on V. In other words, it's just the right size
322 # to create an isomorphism between it and our matrices.
323 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
325 # Turn our matrices into long vectors...
326 vectors
= [ W(m
.list()) for m
in tensor_products
]
328 # Vector space representation of Lyapunov-like matrices
329 # (i.e. vec(L) where L is Luapunov-like).
330 LL_vector
= W
.span(vectors
).complement()
332 # Now construct an ambient MatrixSpace in which to stick our
334 M
= MatrixSpace(V
.base_ring(), V
.dimension())
336 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
342 def lyapunov_rank(K
):
344 Compute the Lyapunov rank (or bilinearity rank) of this cone.
346 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
348 1. The dimension of the Lie algebra of the automorphism group of the
351 2. The dimension of the linear space of all Lyapunov-like
352 transformations on the cone.
356 A closed, convex polyhedral cone.
360 An integer representing the Lyapunov rank of the cone. If the
361 dimension of the ambient vector space is `n`, then the Lyapunov rank
362 will be between `1` and `n` inclusive; however a rank of `n-1` is
363 not possible (see [Orlitzky/Gowda]_).
367 The codimension formula from the second reference is used. We find
368 all pairs `(x,s)` in the complementarity set of `K` such that `x`
369 and `s` are rays of our cone. It is known that these vectors are
370 sufficient to apply the codimension formula. Once we have all such
371 pairs, we "brute force" the codimension formula by finding all
372 linearly-independent `xs^{T}`.
376 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
377 cone and Lyapunov-like transformations, Mathematical Programming, 147
380 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
381 Improper Cone. Work in-progress.
383 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
384 optimality constraints for the cone of positive polynomials,
385 Mathematical Programming, Series B, 129 (2011) 5-31.
389 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
392 sage: positives = Cone([(1,)])
393 sage: lyapunov_rank(positives)
395 sage: quadrant = Cone([(1,0), (0,1)])
396 sage: lyapunov_rank(quadrant)
398 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
399 sage: lyapunov_rank(octant)
402 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
405 sage: R5 = VectorSpace(QQ, 5)
406 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
408 sage: lyapunov_rank(K)
411 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
414 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
415 sage: lyapunov_rank(L31)
418 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
420 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
421 sage: lyapunov_rank(L3infty)
424 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
425 + 1` [Orlitzky/Gowda]_::
427 sage: K = Cone([(1,0,0,0,0)])
428 sage: lyapunov_rank(K)
430 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
433 A subspace (of dimension `m`) in `n` dimensions should have a
434 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
436 sage: e1 = (1,0,0,0,0)
437 sage: neg_e1 = (-1,0,0,0,0)
438 sage: e2 = (0,1,0,0,0)
439 sage: neg_e2 = (0,-1,0,0,0)
440 sage: z = (0,0,0,0,0)
441 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
442 sage: lyapunov_rank(K)
444 sage: K.lattice_dim()**2 - K.dim()*K.codim()
447 The Lyapunov rank should be additive on a product of proper cones
450 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
451 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
452 sage: K = L31.cartesian_product(octant)
453 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
456 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
457 The cone ``K`` in the following example is isomorphic to the nonnegative
458 octant in `\mathbb{R}^{3}`::
460 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
461 sage: lyapunov_rank(K)
464 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
465 itself [Rudolf et al.]_::
467 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
468 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
473 The Lyapunov rank should be additive on a product of proper cones
476 sage: set_random_seed()
477 sage: K1 = random_cone(max_ambient_dim=8,
478 ....: strictly_convex=True,
480 sage: K2 = random_cone(max_ambient_dim=8,
481 ....: strictly_convex=True,
483 sage: K = K1.cartesian_product(K2)
484 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
487 The Lyapunov rank is invariant under a linear isomorphism
490 sage: K1 = random_cone(max_ambient_dim = 8)
491 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
492 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
493 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
496 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
497 itself [Rudolf et al.]_::
499 sage: set_random_seed()
500 sage: K = random_cone(max_ambient_dim=8)
501 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
504 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
505 be any number between `1` and `n` inclusive, excluding `n-1`
506 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
507 trivial cone in a trivial space as well. However, in zero dimensions,
508 the Lyapunov rank of the trivial cone will be zero::
510 sage: set_random_seed()
511 sage: K = random_cone(max_ambient_dim=8,
512 ....: strictly_convex=True,
514 sage: b = lyapunov_rank(K)
515 sage: n = K.lattice_dim()
516 sage: (n == 0 or 1 <= b) and b <= n
521 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
522 Lyapunov rank `n-1` in `n` dimensions::
524 sage: set_random_seed()
525 sage: K = random_cone(max_ambient_dim=8)
526 sage: b = lyapunov_rank(K)
527 sage: n = K.lattice_dim()
531 The calculation of the Lyapunov rank of an improper cone can be
532 reduced to that of a proper cone [Orlitzky/Gowda]_::
534 sage: set_random_seed()
535 sage: K = random_cone(max_ambient_dim=8)
536 sage: actual = lyapunov_rank(K)
537 sage: K_S = _restrict_to_space(K, K.span())
538 sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
539 sage: l = K.lineality()
541 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
542 sage: actual == expected
545 The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
547 sage: set_random_seed()
548 sage: K = random_cone(max_ambient_dim=8)
549 sage: lyapunov_rank(K) == len(LL(K))
552 We can make an imperfect cone perfect by adding a slack variable
553 (a Theorem in [Orlitzky/Gowda]_)::
555 sage: set_random_seed()
556 sage: K = random_cone(max_ambient_dim=8,
557 ....: strictly_convex=True,
559 sage: L = ToricLattice(K.lattice_dim() + 1)
560 sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
561 sage: lyapunov_rank(K) >= K.lattice_dim()
572 # K is not solid, restrict to its span.
573 K
= _restrict_to_space(K
, K
.span())
575 # Non-solid reduction lemma.
579 # K is not pointed, restrict to the span of its dual. Uses a
580 # proposition from our paper, i.e. this is equivalent to K =
581 # _rho(K.dual()).dual().
582 K
= _restrict_to_space(K
, K
.dual().span())
584 # Non-pointed reduction lemma.
592 def is_lyapunov_like(L
,K
):
594 Determine whether or not ``L`` is Lyapunov-like on ``K``.
596 We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
597 L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
598 `\left\langle x,s \right\rangle` in the complementarity set of
599 ``K``. It is known [Orlitzky]_ that this property need only be
600 checked for generators of ``K`` and its dual.
604 - ``L`` -- A linear transformation or matrix.
606 - ``K`` -- A polyhedral closed convex cone.
610 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
611 and ``False`` otherwise.
615 If this function returns ``True``, then ``L`` is Lyapunov-like
616 on ``K``. However, if ``False`` is returned, that could mean one
617 of two things. The first is that ``L`` is definitely not
618 Lyapunov-like on ``K``. The second is more of an "I don't know"
619 answer, returned (for example) if we cannot prove that an inner
624 .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
625 improper cone (preprint).
629 The identity is always Lyapunov-like in a nontrivial space::
631 sage: set_random_seed()
632 sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
633 sage: L = identity_matrix(K.lattice_dim())
634 sage: is_lyapunov_like(L,K)
637 As is the "zero" transformation::
639 sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
640 sage: R = K.lattice().vector_space().base_ring()
641 sage: L = zero_matrix(R, K.lattice_dim())
642 sage: is_lyapunov_like(L,K)
645 Everything in ``LL(K)`` should be Lyapunov-like on ``K``::
647 sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
648 sage: all([is_lyapunov_like(L,K) for L in LL(K)])
652 return all([(L
*x
).inner_product(s
) == 0
653 for (x
,s
) in K
.discrete_complementarity_set()])
656 def random_element(K
):
658 Return a random element of ``K`` from its ambient vector space.
662 The cone ``K`` is specified in terms of its generators, so that
663 ``K`` is equal to the convex conic combination of those generators.
664 To choose a random element of ``K``, we assign random nonnegative
665 coefficients to each generator of ``K`` and construct a new vector
666 from the scaled rays.
668 A vector, rather than a ray, is returned so that the element may
669 have non-integer coordinates. Thus the element may have an
670 arbitrarily small norm.
674 A random element of the trivial cone is zero::
676 sage: set_random_seed()
677 sage: K = Cone([], ToricLattice(0))
678 sage: random_element(K)
680 sage: K = Cone([(0,)])
681 sage: random_element(K)
683 sage: K = Cone([(0,0)])
684 sage: random_element(K)
686 sage: K = Cone([(0,0,0)])
687 sage: random_element(K)
692 Any cone should contain an element of itself::
694 sage: set_random_seed()
695 sage: K = random_cone(max_rays = 8)
696 sage: K.contains(random_element(K))
700 V
= K
.lattice().vector_space()
702 coefficients
= [ F
.random_element().abs() for i
in range(K
.nrays()) ]
703 vector_gens
= map(V
, K
.rays())
704 scaled_gens
= [ coefficients
[i
]*vector_gens
[i
]
705 for i
in range(len(vector_gens
)) ]
707 # Make sure we return a vector. Without the coercion, we might
708 # return ``0`` when ``K`` has no rays.
709 v
= V(sum(scaled_gens
))
713 def positive_operators(K
):
715 Compute generators of the cone of positive operators on this cone.
719 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
720 Each matrix ``P`` in the list should have the property that ``P*x``
721 is an element of ``K`` whenever ``x`` is an element of
722 ``K``. Moreover, any nonnegative linear combination of these
723 matrices shares the same property.
727 The trivial cone in a trivial space has no positive operators::
729 sage: K = Cone([], ToricLattice(0))
730 sage: positive_operators(K)
733 Positive operators on the nonnegative orthant are nonnegative matrices::
735 sage: K = Cone([(1,)])
736 sage: positive_operators(K)
739 sage: K = Cone([(1,0),(0,1)])
740 sage: positive_operators(K)
742 [1 0] [0 1] [0 0] [0 0]
743 [0 0], [0 0], [1 0], [0 1]
746 Every operator is positive on the ambient vector space::
748 sage: K = Cone([(1,),(-1,)])
749 sage: K.is_full_space()
751 sage: positive_operators(K)
754 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
755 sage: K.is_full_space()
757 sage: positive_operators(K)
759 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
760 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
765 A positive operator on a cone should send its generators into the cone::
767 sage: K = random_cone(max_ambient_dim = 6)
768 sage: pi_of_k = positive_operators(K)
769 sage: all([K.contains(p*x) for p in pi_of_k for x in K.rays()])
773 V
= K
.lattice().vector_space()
775 # Sage doesn't think matrices are vectors, so we have to convert
776 # our matrices to vectors explicitly before we can figure out how
777 # many are linearly-indepenedent.
779 # The space W has the same base ring as V, but dimension
780 # dim(V)^2. So it has the same dimension as the space of linear
781 # transformations on V. In other words, it's just the right size
782 # to create an isomorphism between it and our matrices.
783 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
785 G1
= [ V(x
) for x
in K
.rays() ]
786 G2
= [ V(s
) for s
in K
.dual().rays() ]
788 tensor_products
= [ s
.tensor_product(x
) for x
in G1
for s
in G2
]
790 # Turn our matrices into long vectors...
791 vectors
= [ W(m
.list()) for m
in tensor_products
]
793 # Create the *dual* cone of the positive operators, expressed as
795 L
= ToricLattice(W
.dimension())
796 pi_dual
= Cone(vectors
, lattice
=L
)
798 # Now compute the desired cone from its dual...
799 pi_cone
= pi_dual
.dual()
801 # And finally convert its rays back to matrix representations.
802 M
= MatrixSpace(V
.base_ring(), V
.dimension())
804 return [ M(v
.list()) for v
in pi_cone
.rays() ]