]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
e40579fa634dc51afc8eaf83ff5b3e68025180b5
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
= discrete_complementarity_set(K1
)
83 C_of_K2
= discrete_complementarity_set(K2
)
84 if len(C_of_K1
) != len(C_of_K2
):
87 if len(K1
.facets()) != len(K2
.facets()):
96 Restrict ``K`` into its own span, or the span of another cone.
100 - ``K2`` -- another cone whose lattice has the same rank as this
105 A new cone in a sublattice.
109 sage: K = Cone([(1,)])
113 sage: K2 = Cone([(1,0)])
114 sage: _rho(K2).rays()
117 sage: K3 = Cone([(1,0,0)])
118 sage: _rho(K3).rays()
121 sage: _rho(K2) == _rho(K3)
126 The projected cone should always be solid::
128 sage: set_random_seed()
129 sage: K = random_cone(max_ambient_dim = 8)
134 And the resulting cone should live in a space having the same
135 dimension as the space we restricted it to::
137 sage: set_random_seed()
138 sage: K = random_cone(max_ambient_dim = 8)
139 sage: K_S = _rho(K, K.dual() )
140 sage: K_S.lattice_dim() == K.dual().dim()
143 This function should not affect the dimension of a cone::
145 sage: set_random_seed()
146 sage: K = random_cone(max_ambient_dim = 8)
147 sage: K.dim() == _rho(K).dim()
150 Nor should it affect the lineality of a cone::
152 sage: set_random_seed()
153 sage: K = random_cone(max_ambient_dim = 8)
154 sage: K.lineality() == _rho(K).lineality()
157 No matter which space we restrict to, the lineality should not
160 sage: set_random_seed()
161 sage: K = random_cone(max_ambient_dim = 8)
162 sage: K.lineality() >= _rho(K).lineality()
164 sage: K.lineality() >= _rho(K, K.dual()).lineality()
167 If we do this according to our paper, then the result is proper::
169 sage: set_random_seed()
170 sage: K = random_cone(max_ambient_dim = 8,
171 ....: strictly_convex=False,
174 sage: K_SP = _rho(K_S.dual()).dual()
175 sage: K_SP.is_proper()
177 sage: K_SP = _rho(K_S, K_S.dual())
178 sage: K_SP.is_proper()
183 sage: set_random_seed()
184 sage: K = random_cone(max_ambient_dim = 8,
185 ....: strictly_convex=True,
188 sage: K_SP = _rho(K_S.dual()).dual()
189 sage: K_SP.is_proper()
191 sage: K_SP = _rho(K_S, K_S.dual())
192 sage: K_SP.is_proper()
197 sage: set_random_seed()
198 sage: K = random_cone(max_ambient_dim = 8,
199 ....: strictly_convex=False,
202 sage: K_SP = _rho(K_S.dual()).dual()
203 sage: K_SP.is_proper()
205 sage: K_SP = _rho(K_S, K_S.dual())
206 sage: K_SP.is_proper()
211 sage: set_random_seed()
212 sage: K = random_cone(max_ambient_dim = 8,
213 ....: strictly_convex=True,
216 sage: K_SP = _rho(K_S.dual()).dual()
217 sage: K_SP.is_proper()
219 sage: K_SP = _rho(K_S, K_S.dual())
220 sage: K_SP.is_proper()
223 Test Proposition 7 in our paper concerning the duals and
224 restrictions. Generate a random cone, then create a subcone of
225 it. The operation of dual-taking should then commute with rho::
227 sage: set_random_seed()
228 sage: J = random_cone(max_ambient_dim = 8,
230 ....: strictly_convex=False)
231 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
232 sage: K_W_star = _rho(K, J).dual()
233 sage: K_star_W = _rho(K.dual(), J)
234 sage: _basically_the_same(K_W_star, K_star_W)
239 sage: set_random_seed()
240 sage: J = random_cone(max_ambient_dim = 8,
242 ....: strictly_convex=False)
243 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
244 sage: K_W_star = _rho(K, J).dual()
245 sage: K_star_W = _rho(K.dual(), J)
246 sage: _basically_the_same(K_W_star, K_star_W)
251 sage: set_random_seed()
252 sage: J = random_cone(max_ambient_dim = 8,
254 ....: strictly_convex=True)
255 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
256 sage: K_W_star = _rho(K, J).dual()
257 sage: K_star_W = _rho(K.dual(), J)
258 sage: _basically_the_same(K_W_star, K_star_W)
263 sage: set_random_seed()
264 sage: J = random_cone(max_ambient_dim = 8,
266 ....: strictly_convex=True)
267 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
268 sage: K_W_star = _rho(K, J).dual()
269 sage: K_star_W = _rho(K.dual(), J)
270 sage: _basically_the_same(K_W_star, K_star_W)
277 # First we project K onto the span of K2. This will explode if the
278 # rank of ``K2.lattice()`` doesn't match ours.
279 span_K2
= Cone(K2
.rays() + (-K2
).rays(), lattice
=K
.lattice())
280 K
= K
.intersection(span_K2
)
282 # Cheat a little to get the subspace span(K2). The paper uses the
283 # rays of K2 as a basis, but everything is invariant under linear
284 # isomorphism (i.e. a change of basis), and this is a little
286 W
= span_K2
.linear_subspace()
288 # We've already intersected K with the span of K2, so every
289 # generator of K should belong to W now.
290 W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
292 L
= ToricLattice(K2
.dim())
293 return Cone(W_rays
, lattice
=L
)
297 def discrete_complementarity_set(K
):
299 Compute the discrete complementarity set of this cone.
301 The complementarity set of a cone is the set of all orthogonal pairs
302 `(x,s)` such that `x` is in the cone, and `s` is in its dual. The
303 discrete complementarity set is a subset of the complementarity set
304 where `x` and `s` are required to be generators of their respective
307 For polyhedral cones, the discrete complementarity set is always
312 A list of pairs `(x,s)` such that,
314 * Both `x` and `s` are vectors (not rays).
315 * `x` is a generator of this cone.
316 * `s` is a generator of this cone's dual.
317 * `x` and `s` are orthogonal.
321 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
322 Improper Cone. Work in-progress.
326 The discrete complementarity set of the nonnegative orthant consists
327 of pairs of standard basis vectors::
329 sage: K = Cone([(1,0),(0,1)])
330 sage: discrete_complementarity_set(K)
331 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
333 If the cone consists of a single ray, the second components of the
334 discrete complementarity set should generate the orthogonal
335 complement of that ray::
337 sage: K = Cone([(1,0)])
338 sage: discrete_complementarity_set(K)
339 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
340 sage: K = Cone([(1,0,0)])
341 sage: discrete_complementarity_set(K)
342 [((1, 0, 0), (0, 1, 0)),
343 ((1, 0, 0), (0, -1, 0)),
344 ((1, 0, 0), (0, 0, 1)),
345 ((1, 0, 0), (0, 0, -1))]
347 When the cone is the entire space, its dual is the trivial cone, so
348 the discrete complementarity set is empty::
350 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
351 sage: discrete_complementarity_set(K)
354 Likewise when this cone is trivial (its dual is the entire space)::
356 sage: L = ToricLattice(0)
357 sage: K = Cone([], ToricLattice(0))
358 sage: discrete_complementarity_set(K)
363 The complementarity set of the dual can be obtained by switching the
364 components of the complementarity set of the original cone::
366 sage: set_random_seed()
367 sage: K1 = random_cone(max_ambient_dim=6)
369 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
370 sage: actual = discrete_complementarity_set(K1)
371 sage: sorted(actual) == sorted(expected)
374 The pairs in the discrete complementarity set are in fact
377 sage: set_random_seed()
378 sage: K = random_cone(max_ambient_dim=6)
379 sage: dcs = discrete_complementarity_set(K)
380 sage: sum([x.inner_product(s).abs() for (x,s) in dcs])
384 V
= K
.lattice().vector_space()
386 # Convert rays to vectors so that we can compute inner products.
387 xs
= [V(x
) for x
in K
.rays()]
389 # We also convert the generators of the dual cone so that we
390 # return pairs of vectors and not (vector, ray) pairs.
391 ss
= [V(s
) for s
in K
.dual().rays()]
393 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
398 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
403 A list of matrices forming a basis for the space of all
404 Lyapunov-like transformations on the given cone.
408 The trivial cone has no Lyapunov-like transformations::
410 sage: L = ToricLattice(0)
411 sage: K = Cone([], lattice=L)
415 The Lyapunov-like transformations on the nonnegative orthant are
416 simply diagonal matrices::
418 sage: K = Cone([(1,)])
422 sage: K = Cone([(1,0),(0,1)])
429 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
432 [1 0 0] [0 0 0] [0 0 0]
433 [0 0 0] [0 1 0] [0 0 0]
434 [0 0 0], [0 0 0], [0 0 1]
437 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
438 `L^{3}_{\infty}` cones [Rudolf et al.]_::
440 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
448 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
456 If our cone is the entire space, then every transformation on it is
459 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
460 sage: M = MatrixSpace(QQ,2)
461 sage: M.basis() == LL(K)
466 The inner product `\left< L\left(x\right), s \right>` is zero for
467 every pair `\left( x,s \right)` in the discrete complementarity set
470 sage: set_random_seed()
471 sage: K = random_cone(max_ambient_dim=8)
472 sage: C_of_K = discrete_complementarity_set(K)
473 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
474 sage: sum(map(abs, l))
477 The Lyapunov-like transformations on a cone and its dual are related
478 by transposition, but we're not guaranteed to compute transposed
479 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
482 sage: set_random_seed()
483 sage: K = random_cone(max_ambient_dim=8)
484 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
485 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
486 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
487 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
488 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
492 V
= K
.lattice().vector_space()
494 C_of_K
= discrete_complementarity_set(K
)
496 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
498 # Sage doesn't think matrices are vectors, so we have to convert
499 # our matrices to vectors explicitly before we can figure out how
500 # many are linearly-indepenedent.
502 # The space W has the same base ring as V, but dimension
503 # dim(V)^2. So it has the same dimension as the space of linear
504 # transformations on V. In other words, it's just the right size
505 # to create an isomorphism between it and our matrices.
506 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
508 # Turn our matrices into long vectors...
509 vectors
= [ W(m
.list()) for m
in tensor_products
]
511 # Vector space representation of Lyapunov-like matrices
512 # (i.e. vec(L) where L is Luapunov-like).
513 LL_vector
= W
.span(vectors
).complement()
515 # Now construct an ambient MatrixSpace in which to stick our
517 M
= MatrixSpace(V
.base_ring(), V
.dimension())
519 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
525 def lyapunov_rank(K
):
527 Compute the Lyapunov (or bilinearity) rank of this cone.
529 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
531 1. The dimension of the Lie algebra of the automorphism group of the
534 2. The dimension of the linear space of all Lyapunov-like
535 transformations on the cone.
539 A closed, convex polyhedral cone.
543 An integer representing the Lyapunov rank of the cone. If the
544 dimension of the ambient vector space is `n`, then the Lyapunov rank
545 will be between `1` and `n` inclusive; however a rank of `n-1` is
546 not possible (see the first reference).
550 In the references, the cones are always assumed to be proper. We
551 do not impose this restriction.
559 The codimension formula from the second reference is used. We find
560 all pairs `(x,s)` in the complementarity set of `K` such that `x`
561 and `s` are rays of our cone. It is known that these vectors are
562 sufficient to apply the codimension formula. Once we have all such
563 pairs, we "brute force" the codimension formula by finding all
564 linearly-independent `xs^{T}`.
568 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
569 cone and Lyapunov-like transformations, Mathematical Programming, 147
572 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
573 Improper Cone. Work in-progress.
575 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
576 optimality constraints for the cone of positive polynomials,
577 Mathematical Programming, Series B, 129 (2011) 5-31.
581 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
584 sage: positives = Cone([(1,)])
585 sage: lyapunov_rank(positives)
587 sage: quadrant = Cone([(1,0), (0,1)])
588 sage: lyapunov_rank(quadrant)
590 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
591 sage: lyapunov_rank(octant)
594 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
597 sage: R5 = VectorSpace(QQ, 5)
598 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
600 sage: lyapunov_rank(K)
603 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
606 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
607 sage: lyapunov_rank(L31)
610 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
612 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
613 sage: lyapunov_rank(L3infty)
616 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
617 + 1` [Orlitzky/Gowda]_::
619 sage: K = Cone([(1,0,0,0,0)])
620 sage: lyapunov_rank(K)
622 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
625 A subspace (of dimension `m`) in `n` dimensions should have a
626 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
628 sage: e1 = (1,0,0,0,0)
629 sage: neg_e1 = (-1,0,0,0,0)
630 sage: e2 = (0,1,0,0,0)
631 sage: neg_e2 = (0,-1,0,0,0)
632 sage: z = (0,0,0,0,0)
633 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
634 sage: lyapunov_rank(K)
636 sage: K.lattice_dim()**2 - K.dim()*K.codim()
639 The Lyapunov rank should be additive on a product of proper cones
642 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
643 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
644 sage: K = L31.cartesian_product(octant)
645 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
648 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
649 The cone ``K`` in the following example is isomorphic to the nonnegative
650 octant in `\mathbb{R}^{3}`::
652 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
653 sage: lyapunov_rank(K)
656 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
657 itself [Rudolf et al.]_::
659 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
660 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
665 The Lyapunov rank should be additive on a product of proper cones
668 sage: set_random_seed()
669 sage: K1 = random_cone(max_ambient_dim=8,
670 ....: strictly_convex=True,
672 sage: K2 = random_cone(max_ambient_dim=8,
673 ....: strictly_convex=True,
675 sage: K = K1.cartesian_product(K2)
676 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
679 The Lyapunov rank is invariant under a linear isomorphism
682 sage: K1 = random_cone(max_ambient_dim = 8)
683 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
684 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
685 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
688 Just to be sure, test a few more::
690 sage: K1 = random_cone(max_ambient_dim=8,
691 ....: strictly_convex=True,
693 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
694 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
695 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
700 sage: K1 = random_cone(max_ambient_dim=8,
701 ....: strictly_convex=True,
703 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
704 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
705 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
710 sage: K1 = random_cone(max_ambient_dim=8,
711 ....: strictly_convex=False,
713 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
714 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
715 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
720 sage: K1 = random_cone(max_ambient_dim=8,
721 ....: strictly_convex=False,
723 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
724 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
725 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
728 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
729 itself [Rudolf et al.]_::
731 sage: set_random_seed()
732 sage: K = random_cone(max_ambient_dim=8)
733 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
736 Make sure we exercise the non-strictly-convex/non-solid case::
738 sage: set_random_seed()
739 sage: K = random_cone(max_ambient_dim=8,
740 ....: strictly_convex=False,
742 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
745 Let's check the other permutations as well, just to be sure::
747 sage: set_random_seed()
748 sage: K = random_cone(max_ambient_dim=8,
749 ....: strictly_convex=False,
751 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
756 sage: set_random_seed()
757 sage: K = random_cone(max_ambient_dim=8,
758 ....: strictly_convex=True,
760 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
765 sage: set_random_seed()
766 sage: K = random_cone(max_ambient_dim=8,
767 ....: strictly_convex=True,
769 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
772 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
773 be any number between `1` and `n` inclusive, excluding `n-1`
774 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
775 trivial cone in a trivial space as well. However, in zero dimensions,
776 the Lyapunov rank of the trivial cone will be zero::
778 sage: set_random_seed()
779 sage: K = random_cone(max_ambient_dim=8,
780 ....: strictly_convex=True,
782 sage: b = lyapunov_rank(K)
783 sage: n = K.lattice_dim()
784 sage: (n == 0 or 1 <= b) and b <= n
789 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
790 Lyapunov rank `n-1` in `n` dimensions::
792 sage: set_random_seed()
793 sage: K = random_cone(max_ambient_dim=8)
794 sage: b = lyapunov_rank(K)
795 sage: n = K.lattice_dim()
799 The calculation of the Lyapunov rank of an improper cone can be
800 reduced to that of a proper cone [Orlitzky/Gowda]_::
802 sage: set_random_seed()
803 sage: K = random_cone(max_ambient_dim=8)
804 sage: actual = lyapunov_rank(K)
806 sage: K_SP = _rho(K_S.dual()).dual()
807 sage: l = K.lineality()
809 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
810 sage: actual == expected
813 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
815 sage: set_random_seed()
816 sage: K = random_cone(max_ambient_dim=8,
817 ....: strictly_convex=True,
819 sage: lyapunov_rank(K) == len(LL(K))
822 In fact the same can be said of any cone. These additional tests
823 just increase our confidence that the reduction scheme works::
825 sage: set_random_seed()
826 sage: K = random_cone(max_ambient_dim=8,
827 ....: strictly_convex=True,
829 sage: lyapunov_rank(K) == len(LL(K))
834 sage: set_random_seed()
835 sage: K = random_cone(max_ambient_dim=8,
836 ....: strictly_convex=False,
838 sage: lyapunov_rank(K) == len(LL(K))
843 sage: set_random_seed()
844 sage: K = random_cone(max_ambient_dim=8,
845 ....: strictly_convex=False,
847 sage: lyapunov_rank(K) == len(LL(K))
850 Test Theorem 3 in [Orlitzky/Gowda]_::
852 sage: set_random_seed()
853 sage: K = random_cone(max_ambient_dim=8,
854 ....: strictly_convex=True,
856 sage: L = ToricLattice(K.lattice_dim() + 1)
857 sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
858 sage: lyapunov_rank(K) >= K.lattice_dim()
869 # K is not solid, restrict to its span.
873 beta
+= m
*(n
- m
) + (n
- m
)**2
876 # K is not pointed, restrict to the span of its dual. Uses a
877 # proposition from our paper, i.e. this is equivalent to K =
878 # _rho(K.dual()).dual().
879 K
= _rho(K
, K
.dual())