]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
87cdf704580e68b55fa72cd93b8dfa6c1d08a484
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``
16 if K1
.lattice_dim() != K2
.lattice_dim():
19 if K1
.nrays() != K2
.nrays():
22 if K1
.dim() != K2
.dim():
25 if lineality(K1
) != lineality(K2
):
28 if K1
.is_solid() != K2
.is_solid():
31 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
34 if len(LL(K1
)) != len(LL(K2
)):
37 C_of_K1
= discrete_complementarity_set(K1
)
38 C_of_K2
= discrete_complementarity_set(K2
)
39 if len(C_of_K1
) != len(C_of_K2
):
42 if len(K1
.facets()) != len(K2
.facets()):
51 Construct the space `W \times W^{\perp}` isomorphic to the ambient space
52 of ``K`` where `W` is equal to the span of ``K``.
54 V
= K
.lattice().vector_space()
56 # Create the space W \times W^{\perp} isomorphic to V.
57 # First we get an orthogonal (but not normal) basis...
58 M
= matrix(V
.base_field(), K
.rays())
59 W_basis
,_
= M
.gram_schmidt()
61 W
= V
.subspace_with_basis(W_basis
)
62 W_perp
= W
.complement()
64 return W
.cartesian_product(W_perp
)
69 Construct the IPS isomorphism and its inverse from our paper.
71 Given a cone ``K``, the returned isomorphism will split its ambient
72 vector space `V` into a cartesian product `W \times W^{\perp}` where
73 `W` equals the span of ``K``.
75 V
= K
.lattice().vector_space()
77 (W
, W_perp
) = V_iso
.cartesian_factors()
79 # A space equivalent to V, but using our basis.
80 V_user
= V
.subspace_with_basis( W
.basis() + W_perp
.basis() )
83 # Write v in terms of our custom basis, where the first dim(W)
84 # coordinates are for the W-part of the basis.
85 cs
= V_user
.coordinates(v
)
87 w1
= sum([ V_user
.basis()[idx
]*cs
[idx
]
88 for idx
in range(0, W
.dimension()) ])
89 w2
= sum([ V_user
.basis()[idx
]*cs
[idx
]
90 for idx
in range(W
.dimension(), V
.dimension()) ])
92 return V_iso( (w1
, w2
) )
96 # Crash if the arguments are in the wrong spaces.
99 #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
100 #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
101 # for idx in range(0,n-m) ])
103 return sum( pair
.cartesian_factors() )
110 def unrestrict_span(K
, K2
=None):
114 _
,phi_inv
= ips_iso(K2
)
115 V_iso
= iso_space(K2
)
116 (W
, W_perp
) = V_iso
.cartesian_factors()
120 w
= sum([ r
[idx
]*W
.basis()[idx
] for idx
in range(0,len(r
)) ])
121 pair
= V_iso( (w
, W_perp
.zero()) )
122 rays
.append( phi_inv(pair
) )
124 L
= ToricLattice(W
.dimension() + W_perp
.dimension())
126 return Cone(rays
, lattice
=L
)
130 def restrict_span(K
, K2
=None):
132 Restrict ``K`` into its own span, or the span of another cone.
136 - ``K2`` -- another cone whose lattice has the same rank as this cone.
140 A new cone in a sublattice.
144 sage: K = Cone([(1,)])
145 sage: restrict_span(K) == K
148 sage: K2 = Cone([(1,0)])
149 sage: restrict_span(K2).rays()
152 sage: K3 = Cone([(1,0,0)])
153 sage: restrict_span(K3).rays()
156 sage: restrict_span(K2) == restrict_span(K3)
161 The projected cone should always be solid::
163 sage: set_random_seed()
164 sage: K = random_cone(max_dim = 8)
165 sage: K_S = restrict_span(K)
169 And the resulting cone should live in a space having the same
170 dimension as the space we restricted it to::
172 sage: set_random_seed()
173 sage: K = random_cone(max_dim = 8)
174 sage: K_S = restrict_span(K, K.dual() )
175 sage: K_S.lattice_dim() == K.dual().dim()
178 This function has ``unrestrict_span()`` as its inverse::
180 sage: set_random_seed()
181 sage: K = random_cone(max_dim = 8, solid=True)
182 sage: J = restrict_span(K)
183 sage: K == unrestrict_span(J,K)
186 This function should not affect the dimension of a cone::
188 sage: set_random_seed()
189 sage: K = random_cone(max_dim = 8)
190 sage: K.dim() == restrict_span(K).dim()
193 Nor should it affect the lineality of a cone::
195 sage: set_random_seed()
196 sage: K = random_cone(max_dim = 8)
197 sage: lineality(K) == lineality(restrict_span(K))
200 No matter which space we restrict to, the lineality should not
203 sage: set_random_seed()
204 sage: K = random_cone(max_dim = 8)
205 sage: lineality(K) >= lineality(restrict_span(K))
207 sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
210 If we do this according to our paper, then the result is proper::
212 sage: set_random_seed()
213 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
214 sage: K_S = restrict_span(K)
215 sage: P = restrict_span(K_S.dual()).dual()
218 sage: P = restrict_span(K_S, K_S.dual())
224 sage: set_random_seed()
225 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
226 sage: K_S = restrict_span(K)
227 sage: P = restrict_span(K_S.dual()).dual()
230 sage: P = restrict_span(K_S, K_S.dual())
236 sage: set_random_seed()
237 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
238 sage: K_S = restrict_span(K)
239 sage: P = restrict_span(K_S.dual()).dual()
242 sage: P = restrict_span(K_S, K_S.dual())
248 sage: set_random_seed()
249 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
250 sage: K_S = restrict_span(K)
251 sage: P = restrict_span(K_S.dual()).dual()
254 sage: P = restrict_span(K_S, K_S.dual())
258 Test the proposition in our paper concerning the duals, where the
259 subspace `W` is the span of `K^{*}`::
261 sage: set_random_seed()
262 sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
263 sage: K_W = restrict_span(K, K.dual())
264 sage: K_star_W_star = restrict_span(K.dual()).dual()
265 sage: basically_the_same(K_W, K_star_W_star)
270 sage: set_random_seed()
271 sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False)
272 sage: K_W = restrict_span(K, K.dual())
273 sage: K_star_W_star = restrict_span(K.dual()).dual()
274 sage: basically_the_same(K_W, K_star_W_star)
279 sage: set_random_seed()
280 sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True)
281 sage: K_W = restrict_span(K, K.dual())
282 sage: K_star_W_star = restrict_span(K.dual()).dual()
283 sage: basically_the_same(K_W, K_star_W_star)
288 sage: set_random_seed()
289 sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True)
290 sage: K_W = restrict_span(K, K.dual())
291 sage: K_star_W_star = restrict_span(K.dual()).dual()
292 sage: basically_the_same(K_W, K_star_W_star)
300 (W
, W_perp
) = iso_space(K2
).cartesian_factors()
302 ray_pairs
= [ phi(r
) for r
in K
.rays() ]
306 #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
307 # msg = 'Cone has nonzero components in W-perp!'
308 # raise ValueError(msg)
310 # Represent the cone in terms of a basis for W, i.e. with smaller
312 ws
= [ W
.coordinate_vector(w1
) for (w1
, _
) in ray_pairs
]
314 L
= ToricLattice(W
.dimension())
316 return Cone(ws
, lattice
=L
)
322 Compute the lineality of this cone.
324 The lineality of a cone is the dimension of the largest linear
325 subspace contained in that cone.
329 A nonnegative integer; the dimension of the largest subspace
330 contained within this cone.
334 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
335 University Press, Princeton, 1970.
339 The lineality of the nonnegative orthant is zero, since it clearly
342 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
346 However, if we add another ray so that the entire `x`-axis belongs
347 to the cone, then the resulting cone will have lineality one::
349 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
353 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
354 to the dimension of the ambient space (i.e. two)::
356 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
360 Per the definition, the lineality of the trivial cone in a trivial
363 sage: K = Cone([], lattice=ToricLattice(0))
369 The lineality of a cone should be an integer between zero and the
370 dimension of the ambient space, inclusive::
372 sage: set_random_seed()
373 sage: K = random_cone(max_dim = 8)
374 sage: l = lineality(K)
377 sage: (0 <= l) and (l <= K.lattice_dim())
380 A strictly convex cone should have lineality zero::
382 sage: set_random_seed()
383 sage: K = random_cone(max_dim = 8, strictly_convex = True)
388 return K
.linear_subspace().dimension()
393 Compute the codimension of this cone.
395 The codimension of a cone is the dimension of the space of all
396 elements perpendicular to every element of the cone. In other words,
397 the codimension is the difference between the dimension of the
398 ambient space and the dimension of the cone itself.
402 A nonnegative integer representing the dimension of the space of all
403 elements perpendicular to this cone.
407 :meth:`dim`, :meth:`lattice_dim`
411 The codimension of the nonnegative orthant is zero, since the span of
412 its generators equals the entire ambient space::
414 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
418 However, if we remove a ray so that the entire cone is contained
419 within the `x-y`-plane, then the resulting cone will have
420 codimension one, because the `z`-axis is perpendicular to every
421 element of the cone::
423 sage: K = Cone([(1,0,0), (0,1,0)])
427 If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
429 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
433 And if the cone is trivial in any space, then its codimension is
434 equal to the dimension of the ambient space::
436 sage: K = Cone([], lattice=ToricLattice(0))
437 sage: K.lattice_dim()
442 sage: K = Cone([(0,)])
443 sage: K.lattice_dim()
448 sage: K = Cone([(0,0)])
449 sage: K.lattice_dim()
456 The codimension of a cone should be an integer between zero and
457 the dimension of the ambient space, inclusive::
459 sage: set_random_seed()
460 sage: K = random_cone(max_dim = 8)
464 sage: (0 <= c) and (c <= K.lattice_dim())
467 A solid cone should have codimension zero::
469 sage: set_random_seed()
470 sage: K = random_cone(max_dim = 8, solid = True)
474 The codimension of a cone is equal to the lineality of its dual::
476 sage: set_random_seed()
477 sage: K = random_cone(max_dim = 8, solid = True)
478 sage: codim(K) == lineality(K.dual())
482 return (K
.lattice_dim() - K
.dim())
485 def discrete_complementarity_set(K
):
487 Compute the discrete complementarity set of this cone.
489 The complementarity set of this cone is the set of all orthogonal
490 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
491 dual. The discrete complementarity set restricts `x` and `s` to be
492 generators of their respective cones.
496 A list of pairs `(x,s)` such that,
498 * `x` is in this cone.
499 * `x` is a generator of this cone.
500 * `s` is in this cone's dual.
501 * `s` is a generator of this cone's dual.
502 * `x` and `s` are orthogonal.
506 The discrete complementarity set of the nonnegative orthant consists
507 of pairs of standard basis vectors::
509 sage: K = Cone([(1,0),(0,1)])
510 sage: discrete_complementarity_set(K)
511 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
513 If the cone consists of a single ray, the second components of the
514 discrete complementarity set should generate the orthogonal
515 complement of that ray::
517 sage: K = Cone([(1,0)])
518 sage: discrete_complementarity_set(K)
519 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
520 sage: K = Cone([(1,0,0)])
521 sage: discrete_complementarity_set(K)
522 [((1, 0, 0), (0, 1, 0)),
523 ((1, 0, 0), (0, -1, 0)),
524 ((1, 0, 0), (0, 0, 1)),
525 ((1, 0, 0), (0, 0, -1))]
527 When the cone is the entire space, its dual is the trivial cone, so
528 the discrete complementarity set is empty::
530 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
531 sage: discrete_complementarity_set(K)
536 The complementarity set of the dual can be obtained by switching the
537 components of the complementarity set of the original cone::
539 sage: set_random_seed()
540 sage: K1 = random_cone(max_dim=6)
542 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
543 sage: actual = discrete_complementarity_set(K1)
544 sage: sorted(actual) == sorted(expected)
548 V
= K
.lattice().vector_space()
550 # Convert the rays to vectors so that we can compute inner
552 xs
= [V(x
) for x
in K
.rays()]
553 ss
= [V(s
) for s
in K
.dual().rays()]
555 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
560 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
565 A list of matrices forming a basis for the space of all
566 Lyapunov-like transformations on the given cone.
570 The trivial cone has no Lyapunov-like transformations::
572 sage: L = ToricLattice(0)
573 sage: K = Cone([], lattice=L)
577 The Lyapunov-like transformations on the nonnegative orthant are
578 simply diagonal matrices::
580 sage: K = Cone([(1,)])
584 sage: K = Cone([(1,0),(0,1)])
591 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
594 [1 0 0] [0 0 0] [0 0 0]
595 [0 0 0] [0 1 0] [0 0 0]
596 [0 0 0], [0 0 0], [0 0 1]
599 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
600 `L^{3}_{\infty}` cones [Rudolf et al.]_::
602 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
610 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
618 If our cone is the entire space, then every transformation on it is
621 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
622 sage: M = MatrixSpace(QQ,2)
623 sage: M.basis() == LL(K)
628 The inner product `\left< L\left(x\right), s \right>` is zero for
629 every pair `\left( x,s \right)` in the discrete complementarity set
632 sage: set_random_seed()
633 sage: K = random_cone(max_dim=8)
634 sage: C_of_K = discrete_complementarity_set(K)
635 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
636 sage: sum(map(abs, l))
639 The Lyapunov-like transformations on a cone and its dual are related
640 by transposition, but we're not guaranteed to compute transposed
641 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
644 sage: set_random_seed()
645 sage: K = random_cone(max_dim=8)
646 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
647 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
648 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
649 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
650 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
654 V
= K
.lattice().vector_space()
656 C_of_K
= discrete_complementarity_set(K
)
658 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
660 # Sage doesn't think matrices are vectors, so we have to convert
661 # our matrices to vectors explicitly before we can figure out how
662 # many are linearly-indepenedent.
664 # The space W has the same base ring as V, but dimension
665 # dim(V)^2. So it has the same dimension as the space of linear
666 # transformations on V. In other words, it's just the right size
667 # to create an isomorphism between it and our matrices.
668 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
670 # Turn our matrices into long vectors...
671 vectors
= [ W(m
.list()) for m
in tensor_products
]
673 # Vector space representation of Lyapunov-like matrices
674 # (i.e. vec(L) where L is Luapunov-like).
675 LL_vector
= W
.span(vectors
).complement()
677 # Now construct an ambient MatrixSpace in which to stick our
679 M
= MatrixSpace(V
.base_ring(), V
.dimension())
681 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
687 def lyapunov_rank(K
):
689 Compute the Lyapunov (or bilinearity) rank of this cone.
691 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
693 1. The dimension of the Lie algebra of the automorphism group of the
696 2. The dimension of the linear space of all Lyapunov-like
697 transformations on the cone.
701 A closed, convex polyhedral cone.
705 An integer representing the Lyapunov rank of the cone. If the
706 dimension of the ambient vector space is `n`, then the Lyapunov rank
707 will be between `1` and `n` inclusive; however a rank of `n-1` is
708 not possible (see the first reference).
712 In the references, the cones are always assumed to be proper. We
713 do not impose this restriction.
721 The codimension formula from the second reference is used. We find
722 all pairs `(x,s)` in the complementarity set of `K` such that `x`
723 and `s` are rays of our cone. It is known that these vectors are
724 sufficient to apply the codimension formula. Once we have all such
725 pairs, we "brute force" the codimension formula by finding all
726 linearly-independent `xs^{T}`.
730 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
731 cone and Lyapunov-like transformations, Mathematical Programming, 147
734 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
735 Improper Cone. Work in-progress.
737 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
738 optimality constraints for the cone of positive polynomials,
739 Mathematical Programming, Series B, 129 (2011) 5-31.
743 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
746 sage: positives = Cone([(1,)])
747 sage: lyapunov_rank(positives)
749 sage: quadrant = Cone([(1,0), (0,1)])
750 sage: lyapunov_rank(quadrant)
752 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
753 sage: lyapunov_rank(octant)
756 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
759 sage: R5 = VectorSpace(QQ, 5)
760 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
762 sage: lyapunov_rank(K)
765 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
768 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
769 sage: lyapunov_rank(L31)
772 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
774 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
775 sage: lyapunov_rank(L3infty)
778 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
779 + 1` [Orlitzky/Gowda]_::
781 sage: K = Cone([(1,0,0,0,0)])
782 sage: lyapunov_rank(K)
784 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
787 A subspace (of dimension `m`) in `n` dimensions should have a
788 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
790 sage: e1 = (1,0,0,0,0)
791 sage: neg_e1 = (-1,0,0,0,0)
792 sage: e2 = (0,1,0,0,0)
793 sage: neg_e2 = (0,-1,0,0,0)
794 sage: z = (0,0,0,0,0)
795 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
796 sage: lyapunov_rank(K)
798 sage: K.lattice_dim()**2 - K.dim()*codim(K)
801 The Lyapunov rank should be additive on a product of proper cones
804 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
805 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
806 sage: K = L31.cartesian_product(octant)
807 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
810 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
811 The cone ``K`` in the following example is isomorphic to the nonnegative
812 octant in `\mathbb{R}^{3}`::
814 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
815 sage: lyapunov_rank(K)
818 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
819 itself [Rudolf et al.]_::
821 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
822 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
827 The Lyapunov rank should be additive on a product of proper cones
830 sage: set_random_seed()
831 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
832 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
833 sage: K = K1.cartesian_product(K2)
834 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
837 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
838 itself [Rudolf et al.]_::
840 sage: set_random_seed()
841 sage: K = random_cone(max_dim=8)
842 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
845 Make sure we exercise the non-strictly-convex/non-solid case::
847 sage: set_random_seed()
848 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
849 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
852 Let's check the other permutations as well, just to be sure::
854 sage: set_random_seed()
855 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
856 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
861 sage: set_random_seed()
862 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
863 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
868 sage: set_random_seed()
869 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
870 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
873 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
874 be any number between `1` and `n` inclusive, excluding `n-1`
875 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
876 trivial cone in a trivial space as well. However, in zero dimensions,
877 the Lyapunov rank of the trivial cone will be zero::
879 sage: set_random_seed()
880 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
881 sage: b = lyapunov_rank(K)
882 sage: n = K.lattice_dim()
883 sage: (n == 0 or 1 <= b) and b <= n
888 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
889 Lyapunov rank `n-1` in `n` dimensions::
891 sage: set_random_seed()
892 sage: K = random_cone(max_dim=8)
893 sage: b = lyapunov_rank(K)
894 sage: n = K.lattice_dim()
898 The calculation of the Lyapunov rank of an improper cone can be
899 reduced to that of a proper cone [Orlitzky/Gowda]_::
901 sage: set_random_seed()
902 sage: K = random_cone(max_dim=8)
903 sage: actual = lyapunov_rank(K)
904 sage: K_S = restrict_span(K)
905 sage: P = restrict_span(K_S.dual()).dual()
906 sage: l = lineality(K)
908 sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
909 sage: actual == expected
912 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
914 sage: set_random_seed()
915 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
916 sage: lyapunov_rank(K) == len(LL(K))
919 In fact the same can be said of any cone. These additional tests
920 just increase our confidence that the reduction scheme works::
922 sage: set_random_seed()
923 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
924 sage: lyapunov_rank(K) == len(LL(K))
929 sage: set_random_seed()
930 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
931 sage: lyapunov_rank(K) == len(LL(K))
936 sage: set_random_seed()
937 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
938 sage: lyapunov_rank(K) == len(LL(K))
950 # K is not solid, project onto its span.
954 beta
+= m
*(n
- m
) + (n
- m
)**2
957 # K is not pointed, project its dual onto its span.
958 # Uses a proposition from our paper, i.e. this is
959 # equivalent to K = restrict_span(K.dual()).dual()
960 #K = restrict_span(intersect_span(K,K.dual()), K.dual())
961 K
= restrict_span(K
, K
.dual())