]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
61914fa84189f9ba37a0bc3401ff3d6bd2e46461
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 has ``unrestrict_span()`` as its inverse::
179 sage: set_random_seed()
180 sage: K = random_cone(max_dim = 8, solid=True)
181 sage: J = restrict_span(K)
182 sage: K == unrestrict_span(J,K)
185 This function should not affect the dimension of a cone::
187 sage: set_random_seed()
188 sage: K = random_cone(max_dim = 8)
189 sage: K.dim() == restrict_span(K).dim()
192 Nor should it affect the lineality of a cone::
194 sage: set_random_seed()
195 sage: K = random_cone(max_dim = 8)
196 sage: lineality(K) == lineality(restrict_span(K))
199 No matter which space we restrict to, the lineality should not
202 sage: set_random_seed()
203 sage: K = random_cone(max_dim = 8)
204 sage: lineality(K) >= lineality(restrict_span(K))
206 sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
209 If we do this according to our paper, then the result is proper::
211 sage: set_random_seed()
212 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
213 sage: K_S = restrict_span(K)
214 sage: P = restrict_span(K_S.dual()).dual()
217 sage: P = restrict_span(K_S, K_S.dual())
223 sage: set_random_seed()
224 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
225 sage: K_S = restrict_span(K)
226 sage: P = restrict_span(K_S.dual()).dual()
229 sage: P = restrict_span(K_S, K_S.dual())
235 sage: set_random_seed()
236 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
237 sage: K_S = restrict_span(K)
238 sage: P = restrict_span(K_S.dual()).dual()
241 sage: P = restrict_span(K_S, K_S.dual())
247 sage: set_random_seed()
248 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
249 sage: K_S = restrict_span(K)
250 sage: P = restrict_span(K_S.dual()).dual()
253 sage: P = restrict_span(K_S, K_S.dual())
257 Test the proposition in our paper concerning the duals, where the
258 subspace `W` is the span of `K^{*}`::
260 sage: set_random_seed()
261 sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
262 sage: K_W = restrict_span(K, K.dual())
263 sage: K_star_W_star = restrict_span(K.dual()).dual()
264 sage: basically_the_same(K_W, K_star_W_star)
269 sage: set_random_seed()
270 sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False)
271 sage: K_W = restrict_span(K, K.dual())
272 sage: K_star_W_star = restrict_span(K.dual()).dual()
273 sage: basically_the_same(K_W, K_star_W_star)
278 sage: set_random_seed()
279 sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True)
280 sage: K_W = restrict_span(K, K.dual())
281 sage: K_star_W_star = restrict_span(K.dual()).dual()
282 sage: basically_the_same(K_W, K_star_W_star)
287 sage: set_random_seed()
288 sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True)
289 sage: K_W = restrict_span(K, K.dual())
290 sage: K_star_W_star = restrict_span(K.dual()).dual()
291 sage: basically_the_same(K_W, K_star_W_star)
299 (W
, W_perp
) = iso_space(K2
).cartesian_factors()
301 ray_pairs
= [ phi(r
) for r
in K
.rays() ]
305 #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
306 # msg = 'Cone has nonzero components in W-perp!'
307 # raise ValueError(msg)
309 # Represent the cone in terms of a basis for W, i.e. with smaller
311 ws
= [ W
.coordinate_vector(w1
) for (w1
, _
) in ray_pairs
]
313 L
= ToricLattice(W
.dimension())
315 return Cone(ws
, lattice
=L
)
321 Compute the lineality of this cone.
323 The lineality of a cone is the dimension of the largest linear
324 subspace contained in that cone.
328 A nonnegative integer; the dimension of the largest subspace
329 contained within this cone.
333 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
334 University Press, Princeton, 1970.
338 The lineality of the nonnegative orthant is zero, since it clearly
341 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
345 However, if we add another ray so that the entire `x`-axis belongs
346 to the cone, then the resulting cone will have lineality one::
348 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
352 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
353 to the dimension of the ambient space (i.e. two)::
355 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
359 Per the definition, the lineality of the trivial cone in a trivial
362 sage: K = Cone([], lattice=ToricLattice(0))
368 The lineality of a cone should be an integer between zero and the
369 dimension of the ambient space, inclusive::
371 sage: set_random_seed()
372 sage: K = random_cone(max_dim = 8)
373 sage: l = lineality(K)
376 sage: (0 <= l) and (l <= K.lattice_dim())
379 A strictly convex cone should have lineality zero::
381 sage: set_random_seed()
382 sage: K = random_cone(max_dim = 8, strictly_convex = True)
387 return K
.linear_subspace().dimension()
390 def discrete_complementarity_set(K
):
392 Compute the discrete complementarity set of this cone.
394 The complementarity set of this cone is the set of all orthogonal
395 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
396 dual. The discrete complementarity set restricts `x` and `s` to be
397 generators of their respective cones.
401 A list of pairs `(x,s)` such that,
403 * `x` is in this cone.
404 * `x` is a generator of this cone.
405 * `s` is in this cone's dual.
406 * `s` is a generator of this cone's dual.
407 * `x` and `s` are orthogonal.
411 The discrete complementarity set of the nonnegative orthant consists
412 of pairs of standard basis vectors::
414 sage: K = Cone([(1,0),(0,1)])
415 sage: discrete_complementarity_set(K)
416 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
418 If the cone consists of a single ray, the second components of the
419 discrete complementarity set should generate the orthogonal
420 complement of that ray::
422 sage: K = Cone([(1,0)])
423 sage: discrete_complementarity_set(K)
424 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
425 sage: K = Cone([(1,0,0)])
426 sage: discrete_complementarity_set(K)
427 [((1, 0, 0), (0, 1, 0)),
428 ((1, 0, 0), (0, -1, 0)),
429 ((1, 0, 0), (0, 0, 1)),
430 ((1, 0, 0), (0, 0, -1))]
432 When the cone is the entire space, its dual is the trivial cone, so
433 the discrete complementarity set is empty::
435 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
436 sage: discrete_complementarity_set(K)
441 The complementarity set of the dual can be obtained by switching the
442 components of the complementarity set of the original cone::
444 sage: set_random_seed()
445 sage: K1 = random_cone(max_dim=6)
447 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
448 sage: actual = discrete_complementarity_set(K1)
449 sage: sorted(actual) == sorted(expected)
453 V
= K
.lattice().vector_space()
455 # Convert the rays to vectors so that we can compute inner
457 xs
= [V(x
) for x
in K
.rays()]
458 ss
= [V(s
) for s
in K
.dual().rays()]
460 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
465 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
470 A list of matrices forming a basis for the space of all
471 Lyapunov-like transformations on the given cone.
475 The trivial cone has no Lyapunov-like transformations::
477 sage: L = ToricLattice(0)
478 sage: K = Cone([], lattice=L)
482 The Lyapunov-like transformations on the nonnegative orthant are
483 simply diagonal matrices::
485 sage: K = Cone([(1,)])
489 sage: K = Cone([(1,0),(0,1)])
496 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
499 [1 0 0] [0 0 0] [0 0 0]
500 [0 0 0] [0 1 0] [0 0 0]
501 [0 0 0], [0 0 0], [0 0 1]
504 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
505 `L^{3}_{\infty}` cones [Rudolf et al.]_::
507 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
515 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
523 If our cone is the entire space, then every transformation on it is
526 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
527 sage: M = MatrixSpace(QQ,2)
528 sage: M.basis() == LL(K)
533 The inner product `\left< L\left(x\right), s \right>` is zero for
534 every pair `\left( x,s \right)` in the discrete complementarity set
537 sage: set_random_seed()
538 sage: K = random_cone(max_dim=8)
539 sage: C_of_K = discrete_complementarity_set(K)
540 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
541 sage: sum(map(abs, l))
544 The Lyapunov-like transformations on a cone and its dual are related
545 by transposition, but we're not guaranteed to compute transposed
546 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
549 sage: set_random_seed()
550 sage: K = random_cone(max_dim=8)
551 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
552 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
553 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
554 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
555 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
559 V
= K
.lattice().vector_space()
561 C_of_K
= discrete_complementarity_set(K
)
563 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
565 # Sage doesn't think matrices are vectors, so we have to convert
566 # our matrices to vectors explicitly before we can figure out how
567 # many are linearly-indepenedent.
569 # The space W has the same base ring as V, but dimension
570 # dim(V)^2. So it has the same dimension as the space of linear
571 # transformations on V. In other words, it's just the right size
572 # to create an isomorphism between it and our matrices.
573 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
575 # Turn our matrices into long vectors...
576 vectors
= [ W(m
.list()) for m
in tensor_products
]
578 # Vector space representation of Lyapunov-like matrices
579 # (i.e. vec(L) where L is Luapunov-like).
580 LL_vector
= W
.span(vectors
).complement()
582 # Now construct an ambient MatrixSpace in which to stick our
584 M
= MatrixSpace(V
.base_ring(), V
.dimension())
586 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
592 def lyapunov_rank(K
):
594 Compute the Lyapunov (or bilinearity) rank of this cone.
596 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
598 1. The dimension of the Lie algebra of the automorphism group of the
601 2. The dimension of the linear space of all Lyapunov-like
602 transformations on the cone.
606 A closed, convex polyhedral cone.
610 An integer representing the Lyapunov rank of the cone. If the
611 dimension of the ambient vector space is `n`, then the Lyapunov rank
612 will be between `1` and `n` inclusive; however a rank of `n-1` is
613 not possible (see the first reference).
617 In the references, the cones are always assumed to be proper. We
618 do not impose this restriction.
626 The codimension formula from the second reference is used. We find
627 all pairs `(x,s)` in the complementarity set of `K` such that `x`
628 and `s` are rays of our cone. It is known that these vectors are
629 sufficient to apply the codimension formula. Once we have all such
630 pairs, we "brute force" the codimension formula by finding all
631 linearly-independent `xs^{T}`.
635 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
636 cone and Lyapunov-like transformations, Mathematical Programming, 147
639 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
640 Improper Cone. Work in-progress.
642 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
643 optimality constraints for the cone of positive polynomials,
644 Mathematical Programming, Series B, 129 (2011) 5-31.
648 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
651 sage: positives = Cone([(1,)])
652 sage: lyapunov_rank(positives)
654 sage: quadrant = Cone([(1,0), (0,1)])
655 sage: lyapunov_rank(quadrant)
657 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
658 sage: lyapunov_rank(octant)
661 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
664 sage: R5 = VectorSpace(QQ, 5)
665 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
667 sage: lyapunov_rank(K)
670 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
673 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
674 sage: lyapunov_rank(L31)
677 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
679 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
680 sage: lyapunov_rank(L3infty)
683 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
684 + 1` [Orlitzky/Gowda]_::
686 sage: K = Cone([(1,0,0,0,0)])
687 sage: lyapunov_rank(K)
689 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
692 A subspace (of dimension `m`) in `n` dimensions should have a
693 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
695 sage: e1 = (1,0,0,0,0)
696 sage: neg_e1 = (-1,0,0,0,0)
697 sage: e2 = (0,1,0,0,0)
698 sage: neg_e2 = (0,-1,0,0,0)
699 sage: z = (0,0,0,0,0)
700 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
701 sage: lyapunov_rank(K)
703 sage: K.lattice_dim()**2 - K.dim()*codim(K)
706 The Lyapunov rank should be additive on a product of proper cones
709 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
710 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
711 sage: K = L31.cartesian_product(octant)
712 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
715 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
716 The cone ``K`` in the following example is isomorphic to the nonnegative
717 octant in `\mathbb{R}^{3}`::
719 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
720 sage: lyapunov_rank(K)
723 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
724 itself [Rudolf et al.]_::
726 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
727 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
732 The Lyapunov rank should be additive on a product of proper cones
735 sage: set_random_seed()
736 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
737 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
738 sage: K = K1.cartesian_product(K2)
739 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
742 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
743 itself [Rudolf et al.]_::
745 sage: set_random_seed()
746 sage: K = random_cone(max_dim=8)
747 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
750 Make sure we exercise the non-strictly-convex/non-solid case::
752 sage: set_random_seed()
753 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
754 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
757 Let's check the other permutations as well, just to be sure::
759 sage: set_random_seed()
760 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
761 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
766 sage: set_random_seed()
767 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
768 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
773 sage: set_random_seed()
774 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
775 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
778 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
779 be any number between `1` and `n` inclusive, excluding `n-1`
780 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
781 trivial cone in a trivial space as well. However, in zero dimensions,
782 the Lyapunov rank of the trivial cone will be zero::
784 sage: set_random_seed()
785 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
786 sage: b = lyapunov_rank(K)
787 sage: n = K.lattice_dim()
788 sage: (n == 0 or 1 <= b) and b <= n
793 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
794 Lyapunov rank `n-1` in `n` dimensions::
796 sage: set_random_seed()
797 sage: K = random_cone(max_dim=8)
798 sage: b = lyapunov_rank(K)
799 sage: n = K.lattice_dim()
803 The calculation of the Lyapunov rank of an improper cone can be
804 reduced to that of a proper cone [Orlitzky/Gowda]_::
806 sage: set_random_seed()
807 sage: K = random_cone(max_dim=8)
808 sage: actual = lyapunov_rank(K)
809 sage: K_S = restrict_span(K)
810 sage: P = restrict_span(K_S.dual()).dual()
811 sage: l = lineality(K)
813 sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
814 sage: actual == expected
817 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
819 sage: set_random_seed()
820 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
821 sage: lyapunov_rank(K) == len(LL(K))
824 In fact the same can be said of any cone. These additional tests
825 just increase our confidence that the reduction scheme works::
827 sage: set_random_seed()
828 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
829 sage: lyapunov_rank(K) == len(LL(K))
834 sage: set_random_seed()
835 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
836 sage: lyapunov_rank(K) == len(LL(K))
841 sage: set_random_seed()
842 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
843 sage: lyapunov_rank(K) == len(LL(K))
855 # K is not solid, project onto its span.
859 beta
+= m
*(n
- m
) + (n
- m
)**2
862 # K is not pointed, project its dual onto its span.
863 # Uses a proposition from our paper, i.e. this is
864 # equivalent to K = restrict_span(K.dual()).dual()
865 #K = restrict_span(intersect_span(K,K.dual()), K.dual())
866 K
= restrict_span(K
, K
.dual())