]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
c6d26829f32a6efedf927dc542fdcfc453c22e30
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 ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
14 otherwise. This is intended as a lazy way to check whether or not
15 ``K1`` and ``K2`` are linearly isomorphic (i.e. ``A(K1) == K2`` for
16 some invertible linear transformation ``A``).
18 if K1
.lattice_dim() != K2
.lattice_dim():
21 if K1
.nrays() != K2
.nrays():
24 if K1
.dim() != K2
.dim():
27 if lineality(K1
) != lineality(K2
):
30 if K1
.is_solid() != K2
.is_solid():
33 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
36 if len(LL(K1
)) != len(LL(K2
)):
39 C_of_K1
= discrete_complementarity_set(K1
)
40 C_of_K2
= discrete_complementarity_set(K2
)
41 if len(C_of_K1
) != len(C_of_K2
):
44 if len(K1
.facets()) != len(K2
.facets()):
53 Restrict ``K`` into its own span, or the span of another cone.
57 - ``K2`` -- another cone whose lattice has the same rank as this
62 A new cone in a sublattice.
66 sage: K = Cone([(1,)])
70 sage: K2 = Cone([(1,0)])
74 sage: K3 = Cone([(1,0,0)])
78 sage: rho(K2) == rho(K3)
83 The projected cone should always be solid::
85 sage: set_random_seed()
86 sage: K = random_cone(max_dim = 8)
91 And the resulting cone should live in a space having the same
92 dimension as the space we restricted it to::
94 sage: set_random_seed()
95 sage: K = random_cone(max_dim = 8)
96 sage: K_S = rho(K, K.dual() )
97 sage: K_S.lattice_dim() == K.dual().dim()
100 This function should not affect the dimension of a cone::
102 sage: set_random_seed()
103 sage: K = random_cone(max_dim = 8)
104 sage: K.dim() == rho(K).dim()
107 Nor should it affect the lineality of a cone::
109 sage: set_random_seed()
110 sage: K = random_cone(max_dim = 8)
111 sage: lineality(K) == lineality(rho(K))
114 No matter which space we restrict to, the lineality should not
117 sage: set_random_seed()
118 sage: K = random_cone(max_dim = 8)
119 sage: lineality(K) >= lineality(rho(K))
121 sage: lineality(K) >= lineality(rho(K, K.dual()))
124 If we do this according to our paper, then the result is proper::
126 sage: set_random_seed()
127 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
129 sage: K_SP = rho(K_S.dual()).dual()
130 sage: K_SP.is_proper()
132 sage: K_SP = rho(K_S, K_S.dual())
133 sage: K_SP.is_proper()
138 sage: set_random_seed()
139 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
141 sage: K_SP = rho(K_S.dual()).dual()
142 sage: K_SP.is_proper()
144 sage: K_SP = rho(K_S, K_S.dual())
145 sage: K_SP.is_proper()
150 sage: set_random_seed()
151 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
153 sage: K_SP = rho(K_S.dual()).dual()
154 sage: K_SP.is_proper()
156 sage: K_SP = rho(K_S, K_S.dual())
157 sage: K_SP.is_proper()
162 sage: set_random_seed()
163 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
165 sage: K_SP = rho(K_S.dual()).dual()
166 sage: K_SP.is_proper()
168 sage: K_SP = rho(K_S, K_S.dual())
169 sage: K_SP.is_proper()
172 Test the proposition in our paper concerning the duals and
173 restrictions. Generate a random cone, then create a subcone of
174 it. The operation of dual-taking should then commute with rho::
176 sage: set_random_seed()
177 sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=False)
178 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
179 sage: K_W = rho(K, J)
180 sage: K_star_W_star = rho(K.dual(), J).dual()
181 sage: basically_the_same(K_W, K_star_W_star)
186 sage: set_random_seed()
187 sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=False)
188 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
189 sage: K_W = rho(K, J)
190 sage: K_star_W_star = rho(K.dual(), J).dual()
191 sage: basically_the_same(K_W, K_star_W_star)
196 sage: set_random_seed()
197 sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=True)
198 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
199 sage: K_W = rho(K, J)
200 sage: K_star_W_star = rho(K.dual(), J).dual()
201 sage: basically_the_same(K_W, K_star_W_star)
206 sage: set_random_seed()
207 sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=True)
208 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
209 sage: K_W = rho(K, J)
210 sage: K_star_W_star = rho(K.dual(), J).dual()
211 sage: basically_the_same(K_W, K_star_W_star)
218 # First we project K onto the span of K2. This will explode if the
219 # rank of ``K2.lattice()`` doesn't match ours.
220 span_K2
= Cone(K2
.rays() + (-K2
).rays(), lattice
=K
.lattice())
221 K
= K
.intersection(span_K2
)
223 # Cheat a little to get the subspace span(K2). The paper uses the
224 # rays of K2 as a basis, but everything is invariant under linear
225 # isomorphism (i.e. a change of basis), and this is a little
227 W
= span_K2
.linear_subspace()
229 # We've already intersected K with the span of K2, so every
230 # generator of K should belong to W now.
231 W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
233 L
= ToricLattice(K2
.dim())
234 return Cone(W_rays
, lattice
=L
)
240 Compute the lineality of this cone.
242 The lineality of a cone is the dimension of the largest linear
243 subspace contained in that cone.
247 A nonnegative integer; the dimension of the largest subspace
248 contained within this cone.
252 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
253 University Press, Princeton, 1970.
257 The lineality of the nonnegative orthant is zero, since it clearly
260 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
264 However, if we add another ray so that the entire `x`-axis belongs
265 to the cone, then the resulting cone will have lineality one::
267 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
271 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
272 to the dimension of the ambient space (i.e. two)::
274 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
278 Per the definition, the lineality of the trivial cone in a trivial
281 sage: K = Cone([], lattice=ToricLattice(0))
287 The lineality of a cone should be an integer between zero and the
288 dimension of the ambient space, inclusive::
290 sage: set_random_seed()
291 sage: K = random_cone(max_dim = 8)
292 sage: l = lineality(K)
295 sage: (0 <= l) and (l <= K.lattice_dim())
298 A strictly convex cone should have lineality zero::
300 sage: set_random_seed()
301 sage: K = random_cone(max_dim = 8, strictly_convex = True)
306 return K
.linear_subspace().dimension()
311 Compute the codimension of this cone.
313 The codimension of a cone is the dimension of the space of all
314 elements perpendicular to every element of the cone. In other words,
315 the codimension is the difference between the dimension of the
316 ambient space and the dimension of the cone itself.
320 A nonnegative integer representing the dimension of the space of all
321 elements perpendicular to this cone.
325 :meth:`dim`, :meth:`lattice_dim`
329 The codimension of the nonnegative orthant is zero, since the span of
330 its generators equals the entire ambient space::
332 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
336 However, if we remove a ray so that the entire cone is contained
337 within the `x-y`-plane, then the resulting cone will have
338 codimension one, because the `z`-axis is perpendicular to every
339 element of the cone::
341 sage: K = Cone([(1,0,0), (0,1,0)])
345 If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
347 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
351 And if the cone is trivial in any space, then its codimension is
352 equal to the dimension of the ambient space::
354 sage: K = Cone([], lattice=ToricLattice(0))
355 sage: K.lattice_dim()
360 sage: K = Cone([(0,)])
361 sage: K.lattice_dim()
366 sage: K = Cone([(0,0)])
367 sage: K.lattice_dim()
374 The codimension of a cone should be an integer between zero and
375 the dimension of the ambient space, inclusive::
377 sage: set_random_seed()
378 sage: K = random_cone(max_dim = 8)
382 sage: (0 <= c) and (c <= K.lattice_dim())
385 A solid cone should have codimension zero::
387 sage: set_random_seed()
388 sage: K = random_cone(max_dim = 8, solid = True)
392 The codimension of a cone is equal to the lineality of its dual::
394 sage: set_random_seed()
395 sage: K = random_cone(max_dim = 8, solid = True)
396 sage: codim(K) == lineality(K.dual())
400 return (K
.lattice_dim() - K
.dim())
403 def discrete_complementarity_set(K
):
405 Compute the discrete complementarity set of this cone.
407 The complementarity set of this cone is the set of all orthogonal
408 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
409 dual. The discrete complementarity set restricts `x` and `s` to be
410 generators of their respective cones.
414 A list of pairs `(x,s)` such that,
416 * `x` is a generator of this cone.
417 * `s` is a generator of this cone's dual.
418 * `x` and `s` are orthogonal.
422 The discrete complementarity set of the nonnegative orthant consists
423 of pairs of standard basis vectors::
425 sage: K = Cone([(1,0),(0,1)])
426 sage: discrete_complementarity_set(K)
427 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
429 If the cone consists of a single ray, the second components of the
430 discrete complementarity set should generate the orthogonal
431 complement of that ray::
433 sage: K = Cone([(1,0)])
434 sage: discrete_complementarity_set(K)
435 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
436 sage: K = Cone([(1,0,0)])
437 sage: discrete_complementarity_set(K)
438 [((1, 0, 0), (0, 1, 0)),
439 ((1, 0, 0), (0, -1, 0)),
440 ((1, 0, 0), (0, 0, 1)),
441 ((1, 0, 0), (0, 0, -1))]
443 When the cone is the entire space, its dual is the trivial cone, so
444 the discrete complementarity set is empty::
446 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
447 sage: discrete_complementarity_set(K)
452 The complementarity set of the dual can be obtained by switching the
453 components of the complementarity set of the original cone::
455 sage: set_random_seed()
456 sage: K1 = random_cone(max_dim=6)
458 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
459 sage: actual = discrete_complementarity_set(K1)
460 sage: sorted(actual) == sorted(expected)
464 V
= K
.lattice().vector_space()
466 # Convert the rays to vectors so that we can compute inner
468 xs
= [V(x
) for x
in K
.rays()]
469 ss
= [V(s
) for s
in K
.dual().rays()]
471 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
476 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
481 A list of matrices forming a basis for the space of all
482 Lyapunov-like transformations on the given cone.
486 The trivial cone has no Lyapunov-like transformations::
488 sage: L = ToricLattice(0)
489 sage: K = Cone([], lattice=L)
493 The Lyapunov-like transformations on the nonnegative orthant are
494 simply diagonal matrices::
496 sage: K = Cone([(1,)])
500 sage: K = Cone([(1,0),(0,1)])
507 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
510 [1 0 0] [0 0 0] [0 0 0]
511 [0 0 0] [0 1 0] [0 0 0]
512 [0 0 0], [0 0 0], [0 0 1]
515 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
516 `L^{3}_{\infty}` cones [Rudolf et al.]_::
518 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
526 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
534 If our cone is the entire space, then every transformation on it is
537 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
538 sage: M = MatrixSpace(QQ,2)
539 sage: M.basis() == LL(K)
544 The inner product `\left< L\left(x\right), s \right>` is zero for
545 every pair `\left( x,s \right)` in the discrete complementarity set
548 sage: set_random_seed()
549 sage: K = random_cone(max_dim=8)
550 sage: C_of_K = discrete_complementarity_set(K)
551 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
552 sage: sum(map(abs, l))
555 The Lyapunov-like transformations on a cone and its dual are related
556 by transposition, but we're not guaranteed to compute transposed
557 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
560 sage: set_random_seed()
561 sage: K = random_cone(max_dim=8)
562 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
563 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
564 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
565 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
566 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
570 V
= K
.lattice().vector_space()
572 C_of_K
= discrete_complementarity_set(K
)
574 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
576 # Sage doesn't think matrices are vectors, so we have to convert
577 # our matrices to vectors explicitly before we can figure out how
578 # many are linearly-indepenedent.
580 # The space W has the same base ring as V, but dimension
581 # dim(V)^2. So it has the same dimension as the space of linear
582 # transformations on V. In other words, it's just the right size
583 # to create an isomorphism between it and our matrices.
584 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
586 # Turn our matrices into long vectors...
587 vectors
= [ W(m
.list()) for m
in tensor_products
]
589 # Vector space representation of Lyapunov-like matrices
590 # (i.e. vec(L) where L is Luapunov-like).
591 LL_vector
= W
.span(vectors
).complement()
593 # Now construct an ambient MatrixSpace in which to stick our
595 M
= MatrixSpace(V
.base_ring(), V
.dimension())
597 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
603 def lyapunov_rank(K
):
605 Compute the Lyapunov (or bilinearity) rank of this cone.
607 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
609 1. The dimension of the Lie algebra of the automorphism group of the
612 2. The dimension of the linear space of all Lyapunov-like
613 transformations on the cone.
617 A closed, convex polyhedral cone.
621 An integer representing the Lyapunov rank of the cone. If the
622 dimension of the ambient vector space is `n`, then the Lyapunov rank
623 will be between `1` and `n` inclusive; however a rank of `n-1` is
624 not possible (see the first reference).
628 In the references, the cones are always assumed to be proper. We
629 do not impose this restriction.
637 The codimension formula from the second reference is used. We find
638 all pairs `(x,s)` in the complementarity set of `K` such that `x`
639 and `s` are rays of our cone. It is known that these vectors are
640 sufficient to apply the codimension formula. Once we have all such
641 pairs, we "brute force" the codimension formula by finding all
642 linearly-independent `xs^{T}`.
646 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
647 cone and Lyapunov-like transformations, Mathematical Programming, 147
650 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
651 Improper Cone. Work in-progress.
653 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
654 optimality constraints for the cone of positive polynomials,
655 Mathematical Programming, Series B, 129 (2011) 5-31.
659 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
662 sage: positives = Cone([(1,)])
663 sage: lyapunov_rank(positives)
665 sage: quadrant = Cone([(1,0), (0,1)])
666 sage: lyapunov_rank(quadrant)
668 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
669 sage: lyapunov_rank(octant)
672 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
675 sage: R5 = VectorSpace(QQ, 5)
676 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
678 sage: lyapunov_rank(K)
681 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
684 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
685 sage: lyapunov_rank(L31)
688 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
690 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
691 sage: lyapunov_rank(L3infty)
694 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
695 + 1` [Orlitzky/Gowda]_::
697 sage: K = Cone([(1,0,0,0,0)])
698 sage: lyapunov_rank(K)
700 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
703 A subspace (of dimension `m`) in `n` dimensions should have a
704 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
706 sage: e1 = (1,0,0,0,0)
707 sage: neg_e1 = (-1,0,0,0,0)
708 sage: e2 = (0,1,0,0,0)
709 sage: neg_e2 = (0,-1,0,0,0)
710 sage: z = (0,0,0,0,0)
711 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
712 sage: lyapunov_rank(K)
714 sage: K.lattice_dim()**2 - K.dim()*codim(K)
717 The Lyapunov rank should be additive on a product of proper cones
720 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
721 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
722 sage: K = L31.cartesian_product(octant)
723 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
726 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
727 The cone ``K`` in the following example is isomorphic to the nonnegative
728 octant in `\mathbb{R}^{3}`::
730 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
731 sage: lyapunov_rank(K)
734 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
735 itself [Rudolf et al.]_::
737 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
738 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
743 The Lyapunov rank should be additive on a product of proper cones
746 sage: set_random_seed()
747 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
748 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
749 sage: K = K1.cartesian_product(K2)
750 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
753 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
754 itself [Rudolf et al.]_::
756 sage: set_random_seed()
757 sage: K = random_cone(max_dim=8)
758 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
761 Make sure we exercise the non-strictly-convex/non-solid case::
763 sage: set_random_seed()
764 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
765 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
768 Let's check the other permutations as well, just to be sure::
770 sage: set_random_seed()
771 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
772 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
777 sage: set_random_seed()
778 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
779 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
784 sage: set_random_seed()
785 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
786 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
789 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
790 be any number between `1` and `n` inclusive, excluding `n-1`
791 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
792 trivial cone in a trivial space as well. However, in zero dimensions,
793 the Lyapunov rank of the trivial cone will be zero::
795 sage: set_random_seed()
796 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
797 sage: b = lyapunov_rank(K)
798 sage: n = K.lattice_dim()
799 sage: (n == 0 or 1 <= b) and b <= n
804 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
805 Lyapunov rank `n-1` in `n` dimensions::
807 sage: set_random_seed()
808 sage: K = random_cone(max_dim=8)
809 sage: b = lyapunov_rank(K)
810 sage: n = K.lattice_dim()
814 The calculation of the Lyapunov rank of an improper cone can be
815 reduced to that of a proper cone [Orlitzky/Gowda]_::
817 sage: set_random_seed()
818 sage: K = random_cone(max_dim=8)
819 sage: actual = lyapunov_rank(K)
821 sage: K_SP = rho(K_S.dual()).dual()
822 sage: l = lineality(K)
824 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
825 sage: actual == expected
828 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
830 sage: set_random_seed()
831 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
832 sage: lyapunov_rank(K) == len(LL(K))
835 In fact the same can be said of any cone. These additional tests
836 just increase our confidence that the reduction scheme works::
838 sage: set_random_seed()
839 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
840 sage: lyapunov_rank(K) == len(LL(K))
845 sage: set_random_seed()
846 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
847 sage: lyapunov_rank(K) == len(LL(K))
852 sage: set_random_seed()
853 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
854 sage: lyapunov_rank(K) == len(LL(K))
865 # K is not solid, restrict to its span.
869 beta
+= m
*(n
- m
) + (n
- m
)**2
872 # K is not pointed, restrict to the span of its dual. Uses a
873 # proposition from our paper, i.e. this is equivalent to K =
874 # rho(K.dual()).dual().