]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
f4b2244c3486ad6d9ec417bdcdddcc9ffd0f2a62
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 K1
.lineality() != K2
.lineality():
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: K.lineality() == rho(K).lineality()
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: K.lineality() >= rho(K).lineality()
121 sage: K.lineality() >= rho(K, K.dual()).lineality()
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
)
238 def discrete_complementarity_set(K
):
240 Compute the discrete complementarity set of this cone.
242 The complementarity set of this cone is the set of all orthogonal
243 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
244 dual. The discrete complementarity set restricts `x` and `s` to be
245 generators of their respective cones.
249 A list of pairs `(x,s)` such that,
251 * `x` is a generator of this cone.
252 * `s` is a generator of this cone's dual.
253 * `x` and `s` are orthogonal.
257 The discrete complementarity set of the nonnegative orthant consists
258 of pairs of standard basis vectors::
260 sage: K = Cone([(1,0),(0,1)])
261 sage: discrete_complementarity_set(K)
262 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
264 If the cone consists of a single ray, the second components of the
265 discrete complementarity set should generate the orthogonal
266 complement of that ray::
268 sage: K = Cone([(1,0)])
269 sage: discrete_complementarity_set(K)
270 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
271 sage: K = Cone([(1,0,0)])
272 sage: discrete_complementarity_set(K)
273 [((1, 0, 0), (0, 1, 0)),
274 ((1, 0, 0), (0, -1, 0)),
275 ((1, 0, 0), (0, 0, 1)),
276 ((1, 0, 0), (0, 0, -1))]
278 When the cone is the entire space, its dual is the trivial cone, so
279 the discrete complementarity set is empty::
281 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
282 sage: discrete_complementarity_set(K)
287 The complementarity set of the dual can be obtained by switching the
288 components of the complementarity set of the original cone::
290 sage: set_random_seed()
291 sage: K1 = random_cone(max_dim=6)
293 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
294 sage: actual = discrete_complementarity_set(K1)
295 sage: sorted(actual) == sorted(expected)
299 V
= K
.lattice().vector_space()
301 # Convert the rays to vectors so that we can compute inner
303 xs
= [V(x
) for x
in K
.rays()]
304 ss
= [V(s
) for s
in K
.dual().rays()]
306 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
311 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
316 A list of matrices forming a basis for the space of all
317 Lyapunov-like transformations on the given cone.
321 The trivial cone has no Lyapunov-like transformations::
323 sage: L = ToricLattice(0)
324 sage: K = Cone([], lattice=L)
328 The Lyapunov-like transformations on the nonnegative orthant are
329 simply diagonal matrices::
331 sage: K = Cone([(1,)])
335 sage: K = Cone([(1,0),(0,1)])
342 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
345 [1 0 0] [0 0 0] [0 0 0]
346 [0 0 0] [0 1 0] [0 0 0]
347 [0 0 0], [0 0 0], [0 0 1]
350 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
351 `L^{3}_{\infty}` cones [Rudolf et al.]_::
353 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
361 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
369 If our cone is the entire space, then every transformation on it is
372 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
373 sage: M = MatrixSpace(QQ,2)
374 sage: M.basis() == LL(K)
379 The inner product `\left< L\left(x\right), s \right>` is zero for
380 every pair `\left( x,s \right)` in the discrete complementarity set
383 sage: set_random_seed()
384 sage: K = random_cone(max_dim=8)
385 sage: C_of_K = discrete_complementarity_set(K)
386 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
387 sage: sum(map(abs, l))
390 The Lyapunov-like transformations on a cone and its dual are related
391 by transposition, but we're not guaranteed to compute transposed
392 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
395 sage: set_random_seed()
396 sage: K = random_cone(max_dim=8)
397 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
398 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
399 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
400 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
401 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
405 V
= K
.lattice().vector_space()
407 C_of_K
= discrete_complementarity_set(K
)
409 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
411 # Sage doesn't think matrices are vectors, so we have to convert
412 # our matrices to vectors explicitly before we can figure out how
413 # many are linearly-indepenedent.
415 # The space W has the same base ring as V, but dimension
416 # dim(V)^2. So it has the same dimension as the space of linear
417 # transformations on V. In other words, it's just the right size
418 # to create an isomorphism between it and our matrices.
419 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
421 # Turn our matrices into long vectors...
422 vectors
= [ W(m
.list()) for m
in tensor_products
]
424 # Vector space representation of Lyapunov-like matrices
425 # (i.e. vec(L) where L is Luapunov-like).
426 LL_vector
= W
.span(vectors
).complement()
428 # Now construct an ambient MatrixSpace in which to stick our
430 M
= MatrixSpace(V
.base_ring(), V
.dimension())
432 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
438 def lyapunov_rank(K
):
440 Compute the Lyapunov (or bilinearity) rank of this cone.
442 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
444 1. The dimension of the Lie algebra of the automorphism group of the
447 2. The dimension of the linear space of all Lyapunov-like
448 transformations on the cone.
452 A closed, convex polyhedral cone.
456 An integer representing the Lyapunov rank of the cone. If the
457 dimension of the ambient vector space is `n`, then the Lyapunov rank
458 will be between `1` and `n` inclusive; however a rank of `n-1` is
459 not possible (see the first reference).
463 In the references, the cones are always assumed to be proper. We
464 do not impose this restriction.
472 The codimension formula from the second reference is used. We find
473 all pairs `(x,s)` in the complementarity set of `K` such that `x`
474 and `s` are rays of our cone. It is known that these vectors are
475 sufficient to apply the codimension formula. Once we have all such
476 pairs, we "brute force" the codimension formula by finding all
477 linearly-independent `xs^{T}`.
481 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
482 cone and Lyapunov-like transformations, Mathematical Programming, 147
485 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
486 Improper Cone. Work in-progress.
488 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
489 optimality constraints for the cone of positive polynomials,
490 Mathematical Programming, Series B, 129 (2011) 5-31.
494 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
497 sage: positives = Cone([(1,)])
498 sage: lyapunov_rank(positives)
500 sage: quadrant = Cone([(1,0), (0,1)])
501 sage: lyapunov_rank(quadrant)
503 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
504 sage: lyapunov_rank(octant)
507 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
510 sage: R5 = VectorSpace(QQ, 5)
511 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
513 sage: lyapunov_rank(K)
516 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
519 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
520 sage: lyapunov_rank(L31)
523 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
525 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
526 sage: lyapunov_rank(L3infty)
529 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
530 + 1` [Orlitzky/Gowda]_::
532 sage: K = Cone([(1,0,0,0,0)])
533 sage: lyapunov_rank(K)
535 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
538 A subspace (of dimension `m`) in `n` dimensions should have a
539 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
541 sage: e1 = (1,0,0,0,0)
542 sage: neg_e1 = (-1,0,0,0,0)
543 sage: e2 = (0,1,0,0,0)
544 sage: neg_e2 = (0,-1,0,0,0)
545 sage: z = (0,0,0,0,0)
546 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
547 sage: lyapunov_rank(K)
549 sage: K.lattice_dim()**2 - K.dim()*K.codim()
552 The Lyapunov rank should be additive on a product of proper cones
555 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
556 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
557 sage: K = L31.cartesian_product(octant)
558 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
561 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
562 The cone ``K`` in the following example is isomorphic to the nonnegative
563 octant in `\mathbb{R}^{3}`::
565 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
566 sage: lyapunov_rank(K)
569 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
570 itself [Rudolf et al.]_::
572 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
573 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
578 The Lyapunov rank should be additive on a product of proper cones
581 sage: set_random_seed()
582 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
583 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
584 sage: K = K1.cartesian_product(K2)
585 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
588 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
589 itself [Rudolf et al.]_::
591 sage: set_random_seed()
592 sage: K = random_cone(max_dim=8)
593 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
596 Make sure we exercise the non-strictly-convex/non-solid case::
598 sage: set_random_seed()
599 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
600 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
603 Let's check the other permutations as well, just to be sure::
605 sage: set_random_seed()
606 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
607 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
612 sage: set_random_seed()
613 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
614 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
619 sage: set_random_seed()
620 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
621 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
624 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
625 be any number between `1` and `n` inclusive, excluding `n-1`
626 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
627 trivial cone in a trivial space as well. However, in zero dimensions,
628 the Lyapunov rank of the trivial cone will be zero::
630 sage: set_random_seed()
631 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
632 sage: b = lyapunov_rank(K)
633 sage: n = K.lattice_dim()
634 sage: (n == 0 or 1 <= b) and b <= n
639 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
640 Lyapunov rank `n-1` in `n` dimensions::
642 sage: set_random_seed()
643 sage: K = random_cone(max_dim=8)
644 sage: b = lyapunov_rank(K)
645 sage: n = K.lattice_dim()
649 The calculation of the Lyapunov rank of an improper cone can be
650 reduced to that of a proper cone [Orlitzky/Gowda]_::
652 sage: set_random_seed()
653 sage: K = random_cone(max_dim=8)
654 sage: actual = lyapunov_rank(K)
656 sage: K_SP = rho(K_S.dual()).dual()
657 sage: l = K.lineality()
659 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
660 sage: actual == expected
663 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
665 sage: set_random_seed()
666 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
667 sage: lyapunov_rank(K) == len(LL(K))
670 In fact the same can be said of any cone. These additional tests
671 just increase our confidence that the reduction scheme works::
673 sage: set_random_seed()
674 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
675 sage: lyapunov_rank(K) == len(LL(K))
680 sage: set_random_seed()
681 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
682 sage: lyapunov_rank(K) == len(LL(K))
687 sage: set_random_seed()
688 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
689 sage: lyapunov_rank(K) == len(LL(K))
700 # K is not solid, restrict to its span.
704 beta
+= m
*(n
- m
) + (n
- m
)**2
707 # K is not pointed, restrict to the span of its dual. Uses a
708 # proposition from our paper, i.e. this is equivalent to K =
709 # rho(K.dual()).dual().