]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
2d84337fda1cb5275add8171ba8e982559dc1853
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 drop_dependent(vs
):
13 Return the largest linearly-independent subset of ``vs``.
16 # ...for lazy enough definitions of linearly-independent
20 old_V
= VectorSpace(vs
[0].parent().base_field(), 0)
23 new_V
= span(result
+ [v
])
24 if new_V
.dimension() > old_V
.dimension():
31 def basically_the_same(K1
,K2
):
33 ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
36 if K1
.lattice_dim() != K2
.lattice_dim():
39 if K1
.nrays() != K2
.nrays():
42 if K1
.dim() != K2
.dim():
45 if lineality(K1
) != lineality(K2
):
48 if K1
.is_solid() != K2
.is_solid():
51 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
54 if len(LL(K1
)) != len(LL(K2
)):
57 C_of_K1
= discrete_complementarity_set(K1
)
58 C_of_K2
= discrete_complementarity_set(K2
)
59 if len(C_of_K1
) != len(C_of_K2
):
62 if len(K1
.facets()) != len(K2
.facets()):
71 Construct the space `W \times W^{\perp}` isomorphic to the ambient space
72 of ``K`` where `W` is equal to the span of ``K``.
74 V
= K
.lattice().vector_space()
76 # Create the space W \times W^{\perp} isomorphic to V.
77 # First we get an orthogonal (but not normal) basis...
78 M
= matrix(V
.base_field(), K
.rays())
79 W_basis
= drop_dependent(K
.rays())
81 W
= V
.subspace_with_basis(W_basis
)
82 W_perp
= W
.complement()
84 return W
.cartesian_product(W_perp
)
89 Construct the IPS isomorphism and its inverse from our paper.
91 Given a cone ``K``, the returned isomorphism will split its ambient
92 vector space `V` into a cartesian product `W \times W^{\perp}` where
93 `W` equals the span of ``K``.
95 V
= K
.lattice().vector_space()
97 (W
, W_perp
) = V_iso
.cartesian_factors()
99 # A space equivalent to V, but using our basis.
100 V_user
= V
.subspace_with_basis( W
.basis() + W_perp
.basis() )
103 # Write v in terms of our custom basis, where the first dim(W)
104 # coordinates are for the W-part of the basis.
105 cs
= V_user
.coordinates(v
)
107 w1
= sum([ V_user
.basis()[idx
]*cs
[idx
]
108 for idx
in range(0, W
.dimension()) ])
109 w2
= sum([ V_user
.basis()[idx
]*cs
[idx
]
110 for idx
in range(W
.dimension(), V
.dimension()) ])
112 return V_iso( (w1
, w2
) )
116 # Crash if the arguments are in the wrong spaces.
119 #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
120 #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
121 # for idx in range(0,n-m) ])
123 return sum( pair
.cartesian_factors() )
131 Restrict ``K`` into its own span, or the span of another cone.
135 - ``K2`` -- another cone whose lattice has the same rank as this cone.
139 A new cone in a sublattice.
143 sage: K = Cone([(1,)])
144 sage: restrict_span(K) == K
147 sage: K2 = Cone([(1,0)])
148 sage: restrict_span(K2).rays()
151 sage: K3 = Cone([(1,0,0)])
152 sage: restrict_span(K3).rays()
155 sage: restrict_span(K2) == restrict_span(K3)
160 The projected cone should always be solid::
162 sage: set_random_seed()
163 sage: K = random_cone(max_dim = 8)
164 sage: K_S = restrict_span(K)
168 And the resulting cone should live in a space having the same
169 dimension as the space we restricted it to::
171 sage: set_random_seed()
172 sage: K = random_cone(max_dim = 8)
173 sage: K_S = restrict_span(K, K.dual() )
174 sage: K_S.lattice_dim() == K.dual().dim()
177 This function should not affect the dimension of a cone::
179 sage: set_random_seed()
180 sage: K = random_cone(max_dim = 8)
181 sage: K.dim() == restrict_span(K).dim()
184 Nor should it affect the lineality of a cone::
186 sage: set_random_seed()
187 sage: K = random_cone(max_dim = 8)
188 sage: lineality(K) == lineality(restrict_span(K))
191 No matter which space we restrict to, the lineality should not
194 sage: set_random_seed()
195 sage: K = random_cone(max_dim = 8)
196 sage: lineality(K) >= lineality(restrict_span(K))
198 sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
201 If we do this according to our paper, then the result is proper::
203 sage: set_random_seed()
204 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
205 sage: K_S = restrict_span(K)
206 sage: P = restrict_span(K_S.dual()).dual()
209 sage: P = restrict_span(K_S, K_S.dual())
215 sage: set_random_seed()
216 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
217 sage: K_S = restrict_span(K)
218 sage: P = restrict_span(K_S.dual()).dual()
221 sage: P = restrict_span(K_S, K_S.dual())
227 sage: set_random_seed()
228 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
229 sage: K_S = restrict_span(K)
230 sage: P = restrict_span(K_S.dual()).dual()
233 sage: P = restrict_span(K_S, K_S.dual())
239 sage: set_random_seed()
240 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
241 sage: K_S = restrict_span(K)
242 sage: P = restrict_span(K_S.dual()).dual()
245 sage: P = restrict_span(K_S, K_S.dual())
249 Test the proposition in our paper concerning the duals, where the
250 subspace `W` is the span of `K^{*}`::
252 sage: set_random_seed()
253 sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
254 sage: K_W = restrict_span(K, K.dual())
255 sage: K_star_W_star = restrict_span(K.dual()).dual()
256 sage: basically_the_same(K_W, K_star_W_star)
261 sage: set_random_seed()
262 sage: K = random_cone(max_dim = 8, solid=True, 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=False, strictly_convex=True)
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=True, 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)
291 (W
, W_perp
) = iso_space(K2
).cartesian_factors()
293 ray_pairs
= [ phi(r
) for r
in K
.rays() ]
297 #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
298 # msg = 'Cone has nonzero components in W-perp!'
299 # raise ValueError(msg)
301 # Represent the cone in terms of a basis for W, i.e. with smaller
303 ws
= [ W
.coordinate_vector(w1
) for (w1
, _
) in ray_pairs
]
305 L
= ToricLattice(W
.dimension())
307 return Cone(ws
, lattice
=L
)
313 Compute the lineality of this cone.
315 The lineality of a cone is the dimension of the largest linear
316 subspace contained in that cone.
320 A nonnegative integer; the dimension of the largest subspace
321 contained within this cone.
325 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
326 University Press, Princeton, 1970.
330 The lineality of the nonnegative orthant is zero, since it clearly
333 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
337 However, if we add another ray so that the entire `x`-axis belongs
338 to the cone, then the resulting cone will have lineality one::
340 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
344 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
345 to the dimension of the ambient space (i.e. two)::
347 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
351 Per the definition, the lineality of the trivial cone in a trivial
354 sage: K = Cone([], lattice=ToricLattice(0))
360 The lineality of a cone should be an integer between zero and the
361 dimension of the ambient space, inclusive::
363 sage: set_random_seed()
364 sage: K = random_cone(max_dim = 8)
365 sage: l = lineality(K)
368 sage: (0 <= l) and (l <= K.lattice_dim())
371 A strictly convex cone should have lineality zero::
373 sage: set_random_seed()
374 sage: K = random_cone(max_dim = 8, strictly_convex = True)
379 return K
.linear_subspace().dimension()
382 def discrete_complementarity_set(K
):
384 Compute the discrete complementarity set of this cone.
386 The complementarity set of this cone is the set of all orthogonal
387 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
388 dual. The discrete complementarity set restricts `x` and `s` to be
389 generators of their respective cones.
393 A list of pairs `(x,s)` such that,
395 * `x` is in this cone.
396 * `x` is a generator of this cone.
397 * `s` is in this cone's dual.
398 * `s` is a generator of this cone's dual.
399 * `x` and `s` are orthogonal.
403 The discrete complementarity set of the nonnegative orthant consists
404 of pairs of standard basis vectors::
406 sage: K = Cone([(1,0),(0,1)])
407 sage: discrete_complementarity_set(K)
408 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
410 If the cone consists of a single ray, the second components of the
411 discrete complementarity set should generate the orthogonal
412 complement of that ray::
414 sage: K = Cone([(1,0)])
415 sage: discrete_complementarity_set(K)
416 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
417 sage: K = Cone([(1,0,0)])
418 sage: discrete_complementarity_set(K)
419 [((1, 0, 0), (0, 1, 0)),
420 ((1, 0, 0), (0, -1, 0)),
421 ((1, 0, 0), (0, 0, 1)),
422 ((1, 0, 0), (0, 0, -1))]
424 When the cone is the entire space, its dual is the trivial cone, so
425 the discrete complementarity set is empty::
427 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
428 sage: discrete_complementarity_set(K)
433 The complementarity set of the dual can be obtained by switching the
434 components of the complementarity set of the original cone::
436 sage: set_random_seed()
437 sage: K1 = random_cone(max_dim=6)
439 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
440 sage: actual = discrete_complementarity_set(K1)
441 sage: sorted(actual) == sorted(expected)
445 V
= K
.lattice().vector_space()
447 # Convert the rays to vectors so that we can compute inner
449 xs
= [V(x
) for x
in K
.rays()]
450 ss
= [V(s
) for s
in K
.dual().rays()]
452 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
457 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
462 A list of matrices forming a basis for the space of all
463 Lyapunov-like transformations on the given cone.
467 The trivial cone has no Lyapunov-like transformations::
469 sage: L = ToricLattice(0)
470 sage: K = Cone([], lattice=L)
474 The Lyapunov-like transformations on the nonnegative orthant are
475 simply diagonal matrices::
477 sage: K = Cone([(1,)])
481 sage: K = Cone([(1,0),(0,1)])
488 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
491 [1 0 0] [0 0 0] [0 0 0]
492 [0 0 0] [0 1 0] [0 0 0]
493 [0 0 0], [0 0 0], [0 0 1]
496 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
497 `L^{3}_{\infty}` cones [Rudolf et al.]_::
499 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
507 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
515 If our cone is the entire space, then every transformation on it is
518 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
519 sage: M = MatrixSpace(QQ,2)
520 sage: M.basis() == LL(K)
525 The inner product `\left< L\left(x\right), s \right>` is zero for
526 every pair `\left( x,s \right)` in the discrete complementarity set
529 sage: set_random_seed()
530 sage: K = random_cone(max_dim=8)
531 sage: C_of_K = discrete_complementarity_set(K)
532 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
533 sage: sum(map(abs, l))
536 The Lyapunov-like transformations on a cone and its dual are related
537 by transposition, but we're not guaranteed to compute transposed
538 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
541 sage: set_random_seed()
542 sage: K = random_cone(max_dim=8)
543 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
544 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
545 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
546 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
547 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
551 V
= K
.lattice().vector_space()
553 C_of_K
= discrete_complementarity_set(K
)
555 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
557 # Sage doesn't think matrices are vectors, so we have to convert
558 # our matrices to vectors explicitly before we can figure out how
559 # many are linearly-indepenedent.
561 # The space W has the same base ring as V, but dimension
562 # dim(V)^2. So it has the same dimension as the space of linear
563 # transformations on V. In other words, it's just the right size
564 # to create an isomorphism between it and our matrices.
565 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
567 # Turn our matrices into long vectors...
568 vectors
= [ W(m
.list()) for m
in tensor_products
]
570 # Vector space representation of Lyapunov-like matrices
571 # (i.e. vec(L) where L is Luapunov-like).
572 LL_vector
= W
.span(vectors
).complement()
574 # Now construct an ambient MatrixSpace in which to stick our
576 M
= MatrixSpace(V
.base_ring(), V
.dimension())
578 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
584 def lyapunov_rank(K
):
586 Compute the Lyapunov (or bilinearity) rank of this cone.
588 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
590 1. The dimension of the Lie algebra of the automorphism group of the
593 2. The dimension of the linear space of all Lyapunov-like
594 transformations on the cone.
598 A closed, convex polyhedral cone.
602 An integer representing the Lyapunov rank of the cone. If the
603 dimension of the ambient vector space is `n`, then the Lyapunov rank
604 will be between `1` and `n` inclusive; however a rank of `n-1` is
605 not possible (see the first reference).
609 In the references, the cones are always assumed to be proper. We
610 do not impose this restriction.
618 The codimension formula from the second reference is used. We find
619 all pairs `(x,s)` in the complementarity set of `K` such that `x`
620 and `s` are rays of our cone. It is known that these vectors are
621 sufficient to apply the codimension formula. Once we have all such
622 pairs, we "brute force" the codimension formula by finding all
623 linearly-independent `xs^{T}`.
627 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
628 cone and Lyapunov-like transformations, Mathematical Programming, 147
631 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
632 Improper Cone. Work in-progress.
634 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
635 optimality constraints for the cone of positive polynomials,
636 Mathematical Programming, Series B, 129 (2011) 5-31.
640 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
643 sage: positives = Cone([(1,)])
644 sage: lyapunov_rank(positives)
646 sage: quadrant = Cone([(1,0), (0,1)])
647 sage: lyapunov_rank(quadrant)
649 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
650 sage: lyapunov_rank(octant)
653 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
656 sage: R5 = VectorSpace(QQ, 5)
657 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
659 sage: lyapunov_rank(K)
662 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
665 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
666 sage: lyapunov_rank(L31)
669 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
671 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
672 sage: lyapunov_rank(L3infty)
675 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
676 + 1` [Orlitzky/Gowda]_::
678 sage: K = Cone([(1,0,0,0,0)])
679 sage: lyapunov_rank(K)
681 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
684 A subspace (of dimension `m`) in `n` dimensions should have a
685 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
687 sage: e1 = (1,0,0,0,0)
688 sage: neg_e1 = (-1,0,0,0,0)
689 sage: e2 = (0,1,0,0,0)
690 sage: neg_e2 = (0,-1,0,0,0)
691 sage: z = (0,0,0,0,0)
692 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
693 sage: lyapunov_rank(K)
695 sage: K.lattice_dim()**2 - K.dim()*codim(K)
698 The Lyapunov rank should be additive on a product of proper cones
701 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
702 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
703 sage: K = L31.cartesian_product(octant)
704 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
707 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
708 The cone ``K`` in the following example is isomorphic to the nonnegative
709 octant in `\mathbb{R}^{3}`::
711 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
712 sage: lyapunov_rank(K)
715 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
716 itself [Rudolf et al.]_::
718 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
719 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
724 The Lyapunov rank should be additive on a product of proper cones
727 sage: set_random_seed()
728 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
729 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
730 sage: K = K1.cartesian_product(K2)
731 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
734 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
735 itself [Rudolf et al.]_::
737 sage: set_random_seed()
738 sage: K = random_cone(max_dim=8)
739 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
742 Make sure we exercise the non-strictly-convex/non-solid case::
744 sage: set_random_seed()
745 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
746 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
749 Let's check the other permutations as well, just to be sure::
751 sage: set_random_seed()
752 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
753 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
758 sage: set_random_seed()
759 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
760 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
765 sage: set_random_seed()
766 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
767 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
770 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
771 be any number between `1` and `n` inclusive, excluding `n-1`
772 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
773 trivial cone in a trivial space as well. However, in zero dimensions,
774 the Lyapunov rank of the trivial cone will be zero::
776 sage: set_random_seed()
777 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
778 sage: b = lyapunov_rank(K)
779 sage: n = K.lattice_dim()
780 sage: (n == 0 or 1 <= b) and b <= n
785 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
786 Lyapunov rank `n-1` in `n` dimensions::
788 sage: set_random_seed()
789 sage: K = random_cone(max_dim=8)
790 sage: b = lyapunov_rank(K)
791 sage: n = K.lattice_dim()
795 The calculation of the Lyapunov rank of an improper cone can be
796 reduced to that of a proper cone [Orlitzky/Gowda]_::
798 sage: set_random_seed()
799 sage: K = random_cone(max_dim=8)
800 sage: actual = lyapunov_rank(K)
801 sage: K_S = restrict_span(K)
802 sage: P = restrict_span(K_S.dual()).dual()
803 sage: l = lineality(K)
805 sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
806 sage: actual == expected
809 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
811 sage: set_random_seed()
812 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
813 sage: lyapunov_rank(K) == len(LL(K))
816 In fact the same can be said of any cone. These additional tests
817 just increase our confidence that the reduction scheme works::
819 sage: set_random_seed()
820 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
821 sage: lyapunov_rank(K) == len(LL(K))
826 sage: set_random_seed()
827 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
828 sage: lyapunov_rank(K) == len(LL(K))
833 sage: set_random_seed()
834 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
835 sage: lyapunov_rank(K) == len(LL(K))
847 # K is not solid, project onto its span.
851 beta
+= m
*(n
- m
) + (n
- m
)**2
854 # K is not pointed, project its dual onto its span.
855 # Uses a proposition from our paper, i.e. this is
856 # equivalent to K = restrict_span(K.dual()).dual()
857 #K = restrict_span(intersect_span(K,K.dual()), K.dual())
858 K
= restrict_span(K
, K
.dual())