]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
6fb15ae21e242085b223eac729c857b8d04b8189
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():
33 Construct the space `W \times W^{\perp}` isomorphic to the ambient space
34 of ``K`` where `W` is equal to the span of ``K``.
36 V
= K
.lattice().vector_space()
38 # Create the space W \times W^{\perp} isomorphic to V.
39 W_basis
= drop_dependent(K
.rays())
40 W
= V
.subspace_with_basis(W_basis
)
41 W_perp
= W
.complement()
43 return W
.cartesian_product(W_perp
)
48 Construct the IPS isomorphism and its inverse from our paper.
50 Given a cone ``K``, the returned isomorphism will split its ambient
51 vector space `V` into a cartesian product `W \times W^{\perp}` where
52 `W` equals the span of ``K``.
54 V
= K
.lattice().vector_space()
56 (W
, W_perp
) = V_iso
.cartesian_factors()
58 # A space equivalent to V, but using our basis.
59 V_user
= V
.subspace_with_basis( W
.basis() + W_perp
.basis() )
62 # Write v in terms of our custom basis, where the first dim(W)
63 # coordinates are for the W-part of the basis.
64 cs
= V_user
.coordinates(v
)
66 w1
= sum([ V_user
.basis()[idx
]*cs
[idx
]
67 for idx
in range(0, W
.dimension()) ])
68 w2
= sum([ V_user
.basis()[idx
]*cs
[idx
]
69 for idx
in range(W
.dimension(), V
.dimension()) ])
71 return V_iso( (w1
, w2
) )
75 # Crash if the arguments are in the wrong spaces.
78 #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
79 #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
80 # for idx in range(0,n-m) ])
82 return sum( pair
.cartesian_factors() )
89 def unrestrict_span(K
, K2
=None):
93 _
,phi_inv
= ips_iso(K2
)
95 (W
, W_perp
) = V_iso
.cartesian_factors()
99 w
= sum([ r
[idx
]*W
.basis()[idx
] for idx
in range(0,len(r
)) ])
100 pair
= V_iso( (w
, W_perp
.zero()) )
101 rays
.append( phi_inv(pair
) )
103 L
= ToricLattice(W
.dimension() + W_perp
.dimension())
105 return Cone(rays
, lattice
=L
)
109 def intersect_span(K1
, K2
):
111 Return a new cone obtained by intersecting ``K1`` with the span of ``K2``.
115 if L
.rank() != K2
.lattice().rank():
116 raise ValueError('K1 and K2 must belong to lattices of the same rank.')
118 SL_gens
= list(K2
.rays())
119 span_K2_gens
= SL_gens
+ [ -g
for g
in SL_gens
]
121 # The lattices have the same rank (see above) so this should work.
122 span_K2
= Cone(span_K2_gens
, L
)
123 return K1
.intersection(span_K2
)
127 def restrict_span(K
, K2
=None):
129 Restrict ``K`` into its own span, or the span of another cone.
133 - ``K2`` -- another cone whose lattice has the same rank as this cone.
137 A new cone in a sublattice.
141 sage: K = Cone([(1,)])
142 sage: restrict_span(K) == K
145 sage: K2 = Cone([(1,0)])
146 sage: restrict_span(K2).rays()
149 sage: K3 = Cone([(1,0,0)])
150 sage: restrict_span(K3).rays()
153 sage: restrict_span(K2) == restrict_span(K3)
158 The projected cone should always be solid::
160 sage: set_random_seed()
161 sage: K = random_cone(max_dim = 10)
162 sage: K_S = restrict_span(K)
166 And the resulting cone should live in a space having the same
167 dimension as the space we restricted it to::
169 sage: set_random_seed()
170 sage: K = random_cone(max_dim = 10)
171 sage: K_S = restrict_span( intersect_span(K, K.dual()), K.dual() )
172 sage: K_S.lattice_dim() == K.dual().dim()
175 This function has ``unrestrict_span()`` as its inverse::
177 sage: set_random_seed()
178 sage: K = random_cone(max_dim = 10, solid=True)
179 sage: J = restrict_span(K)
180 sage: K == unrestrict_span(J,K)
183 This function should not affect the dimension of a cone::
185 sage: set_random_seed()
186 sage: K = random_cone(max_dim = 10)
187 sage: K.dim() == restrict_span(K).dim()
190 Nor should it affect the lineality of a cone::
192 sage: set_random_seed()
193 sage: K = random_cone(max_dim = 10)
194 sage: lineality(K) == lineality(restrict_span(K))
197 No matter which space we restrict to, the lineality should not
200 sage: set_random_seed()
201 sage: K = random_cone(max_dim = 10)
202 sage: J = intersect_span(K, K.dual())
203 sage: lineality(K) >= lineality(restrict_span(J, K.dual()))
206 If we do this according to our paper, then the result is proper::
208 sage: set_random_seed()
209 sage: K = random_cone(max_dim = 10)
210 sage: K_S = restrict_span(K)
211 sage: P = restrict_span(K_S.dual()).dual()
215 If ``K`` is strictly convex, then both ``K_W`` and
216 ``K_star_W.dual()`` should equal ``K`` (after we unrestrict)::
218 sage: set_random_seed()
219 sage: K = random_cone(max_dim = 10, strictly_convex=True)
220 sage: K_W = restrict_span(intersect_span(K,K.dual()), K.dual())
221 sage: K_star_W_star = restrict_span(K.dual()).dual()
222 sage: j1 = unrestrict_span(K_W, K.dual())
223 sage: j2 = unrestrict_span(K_star_W_star, K.dual())
228 sage: K; [ list(r) for r in K.rays() ]
230 Test the proposition in our paper concerning the duals, where the
231 subspace `W` is the span of `K^{*}`::
233 sage: set_random_seed()
234 sage: K = random_cone(max_dim = 10, solid=False, strictly_convex=False)
235 sage: K_W = restrict_span(intersect_span(K,K.dual()), K.dual())
236 sage: K_star_W_star = restrict_span(K.dual(), K.dual()).dual()
237 sage: K_W.nrays() == K_star_W_star.nrays()
239 sage: K_W.dim() == K_star_W_star.dim()
241 sage: lineality(K_W) == lineality(K_star_W_star)
243 sage: K_W.is_solid() == K_star_W_star.is_solid()
245 sage: K_W.is_strictly_convex() == K_star_W_star.is_strictly_convex()
253 (W
, W_perp
) = iso_space(K2
).cartesian_factors()
255 ray_pairs
= [ phi(r
) for r
in K
.rays() ]
257 if any([ w2
!= W_perp
.zero() for (_
, w2
) in ray_pairs
]):
258 msg
= 'Cone has nonzero components in W-perp!'
259 raise ValueError(msg
)
261 # Represent the cone in terms of a basis for W, i.e. with smaller
263 ws
= [ W
.coordinate_vector(w1
) for (w1
, _
) in ray_pairs
]
265 L
= ToricLattice(W
.dimension())
267 return Cone(ws
, lattice
=L
)
273 Compute the lineality of this cone.
275 The lineality of a cone is the dimension of the largest linear
276 subspace contained in that cone.
280 A nonnegative integer; the dimension of the largest subspace
281 contained within this cone.
285 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
286 University Press, Princeton, 1970.
290 The lineality of the nonnegative orthant is zero, since it clearly
293 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
297 However, if we add another ray so that the entire `x`-axis belongs
298 to the cone, then the resulting cone will have lineality one::
300 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
304 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
305 to the dimension of the ambient space (i.e. two)::
307 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
311 Per the definition, the lineality of the trivial cone in a trivial
314 sage: K = Cone([], lattice=ToricLattice(0))
320 The lineality of a cone should be an integer between zero and the
321 dimension of the ambient space, inclusive::
323 sage: set_random_seed()
324 sage: K = random_cone(max_dim = 10)
325 sage: l = lineality(K)
328 sage: (0 <= l) and (l <= K.lattice_dim())
331 A strictly convex cone should have lineality zero::
333 sage: set_random_seed()
334 sage: K = random_cone(max_dim = 10, strictly_convex = True)
339 return K
.linear_subspace().dimension()
344 Compute the codimension of this cone.
346 The codimension of a cone is the dimension of the space of all
347 elements perpendicular to every element of the cone. In other words,
348 the codimension is the difference between the dimension of the
349 ambient space and the dimension of the cone itself.
353 A nonnegative integer representing the dimension of the space of all
354 elements perpendicular to this cone.
358 :meth:`dim`, :meth:`lattice_dim`
362 The codimension of the nonnegative orthant is zero, since the span of
363 its generators equals the entire ambient space::
365 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
369 However, if we remove a ray so that the entire cone is contained
370 within the `x-y`-plane, then the resulting cone will have
371 codimension one, because the `z`-axis is perpendicular to every
372 element of the cone::
374 sage: K = Cone([(1,0,0), (0,1,0)])
378 If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
380 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
384 And if the cone is trivial in any space, then its codimension is
385 equal to the dimension of the ambient space::
387 sage: K = Cone([], lattice=ToricLattice(0))
391 sage: K = Cone([(0,)])
395 sage: K = Cone([(0,0)])
401 The codimension of a cone should be an integer between zero and
402 the dimension of the ambient space, inclusive::
404 sage: set_random_seed()
405 sage: K = random_cone(max_dim = 10)
409 sage: (0 <= c) and (c <= K.lattice_dim())
412 A solid cone should have codimension zero::
414 sage: set_random_seed()
415 sage: K = random_cone(max_dim = 10, solid = True)
419 The codimension of a cone is equal to the lineality of its dual::
421 sage: set_random_seed()
422 sage: K = random_cone(max_dim = 10, solid = True)
423 sage: codim(K) == lineality(K.dual())
427 return (K
.lattice_dim() - K
.dim())
430 def discrete_complementarity_set(K
):
432 Compute the discrete complementarity set of this cone.
434 The complementarity set of this cone is the set of all orthogonal
435 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
436 dual. The discrete complementarity set restricts `x` and `s` to be
437 generators of their respective cones.
441 A list of pairs `(x,s)` such that,
443 * `x` is in this cone.
444 * `x` is a generator of this cone.
445 * `s` is in this cone's dual.
446 * `s` is a generator of this cone's dual.
447 * `x` and `s` are orthogonal.
451 The discrete complementarity set of the nonnegative orthant consists
452 of pairs of standard basis vectors::
454 sage: K = Cone([(1,0),(0,1)])
455 sage: discrete_complementarity_set(K)
456 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
458 If the cone consists of a single ray, the second components of the
459 discrete complementarity set should generate the orthogonal
460 complement of that ray::
462 sage: K = Cone([(1,0)])
463 sage: discrete_complementarity_set(K)
464 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
465 sage: K = Cone([(1,0,0)])
466 sage: discrete_complementarity_set(K)
467 [((1, 0, 0), (0, 1, 0)),
468 ((1, 0, 0), (0, -1, 0)),
469 ((1, 0, 0), (0, 0, 1)),
470 ((1, 0, 0), (0, 0, -1))]
472 When the cone is the entire space, its dual is the trivial cone, so
473 the discrete complementarity set is empty::
475 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
476 sage: discrete_complementarity_set(K)
481 The complementarity set of the dual can be obtained by switching the
482 components of the complementarity set of the original cone::
484 sage: set_random_seed()
485 sage: K1 = random_cone(max_dim=6)
487 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
488 sage: actual = discrete_complementarity_set(K1)
489 sage: sorted(actual) == sorted(expected)
493 V
= K
.lattice().vector_space()
495 # Convert the rays to vectors so that we can compute inner
497 xs
= [V(x
) for x
in K
.rays()]
498 ss
= [V(s
) for s
in K
.dual().rays()]
500 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
505 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
510 A list of matrices forming a basis for the space of all
511 Lyapunov-like transformations on the given cone.
515 The trivial cone has no Lyapunov-like transformations::
517 sage: L = ToricLattice(0)
518 sage: K = Cone([], lattice=L)
522 The Lyapunov-like transformations on the nonnegative orthant are
523 simply diagonal matrices::
525 sage: K = Cone([(1,)])
529 sage: K = Cone([(1,0),(0,1)])
536 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
539 [1 0 0] [0 0 0] [0 0 0]
540 [0 0 0] [0 1 0] [0 0 0]
541 [0 0 0], [0 0 0], [0 0 1]
544 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
545 `L^{3}_{\infty}` cones [Rudolf et al.]_::
547 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
555 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
565 The inner product `\left< L\left(x\right), s \right>` is zero for
566 every pair `\left( x,s \right)` in the discrete complementarity set
569 sage: set_random_seed()
570 sage: K = random_cone(max_dim=8, max_rays=10)
571 sage: C_of_K = discrete_complementarity_set(K)
572 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
573 sage: sum(map(abs, l))
576 The Lyapunov-like transformations on a cone and its dual are related
577 by transposition, but we're not guaranteed to compute transposed
578 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
581 sage: set_random_seed()
582 sage: K = random_cone(max_dim=8, max_rays=10)
583 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
584 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
585 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
586 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
587 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
591 V
= K
.lattice().vector_space()
593 C_of_K
= discrete_complementarity_set(K
)
595 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
597 # Sage doesn't think matrices are vectors, so we have to convert
598 # our matrices to vectors explicitly before we can figure out how
599 # many are linearly-indepenedent.
601 # The space W has the same base ring as V, but dimension
602 # dim(V)^2. So it has the same dimension as the space of linear
603 # transformations on V. In other words, it's just the right size
604 # to create an isomorphism between it and our matrices.
605 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
607 # Turn our matrices into long vectors...
608 vectors
= [ W(m
.list()) for m
in tensor_products
]
610 # Vector space representation of Lyapunov-like matrices
611 # (i.e. vec(L) where L is Luapunov-like).
612 LL_vector
= W
.span(vectors
).complement()
614 # Now construct an ambient MatrixSpace in which to stick our
616 M
= MatrixSpace(V
.base_ring(), V
.dimension())
618 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
624 def lyapunov_rank(K
):
626 Compute the Lyapunov (or bilinearity) rank of this cone.
628 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
630 1. The dimension of the Lie algebra of the automorphism group of the
633 2. The dimension of the linear space of all Lyapunov-like
634 transformations on the cone.
638 A closed, convex polyhedral cone.
642 An integer representing the Lyapunov rank of the cone. If the
643 dimension of the ambient vector space is `n`, then the Lyapunov rank
644 will be between `1` and `n` inclusive; however a rank of `n-1` is
645 not possible (see the first reference).
649 In the references, the cones are always assumed to be proper. We
650 do not impose this restriction.
658 The codimension formula from the second reference is used. We find
659 all pairs `(x,s)` in the complementarity set of `K` such that `x`
660 and `s` are rays of our cone. It is known that these vectors are
661 sufficient to apply the codimension formula. Once we have all such
662 pairs, we "brute force" the codimension formula by finding all
663 linearly-independent `xs^{T}`.
667 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
668 cone and Lyapunov-like transformations, Mathematical Programming, 147
671 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
672 Improper Cone. Work in-progress.
674 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
675 optimality constraints for the cone of positive polynomials,
676 Mathematical Programming, Series B, 129 (2011) 5-31.
680 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
683 sage: positives = Cone([(1,)])
684 sage: lyapunov_rank(positives)
686 sage: quadrant = Cone([(1,0), (0,1)])
687 sage: lyapunov_rank(quadrant)
689 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
690 sage: lyapunov_rank(octant)
693 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
696 sage: R5 = VectorSpace(QQ, 5)
697 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
699 sage: lyapunov_rank(K)
702 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
705 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
706 sage: lyapunov_rank(L31)
709 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
711 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
712 sage: lyapunov_rank(L3infty)
715 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
716 + 1` [Orlitzky/Gowda]_::
718 sage: K = Cone([(1,0,0,0,0)])
719 sage: lyapunov_rank(K)
721 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
724 A subspace (of dimension `m`) in `n` dimensions should have a
725 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
727 sage: e1 = (1,0,0,0,0)
728 sage: neg_e1 = (-1,0,0,0,0)
729 sage: e2 = (0,1,0,0,0)
730 sage: neg_e2 = (0,-1,0,0,0)
731 sage: z = (0,0,0,0,0)
732 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
733 sage: lyapunov_rank(K)
735 sage: K.lattice_dim()**2 - K.dim()*codim(K)
738 The Lyapunov rank should be additive on a product of proper cones
741 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
742 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
743 sage: K = L31.cartesian_product(octant)
744 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
747 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
748 The cone ``K`` in the following example is isomorphic to the nonnegative
749 octant in `\mathbb{R}^{3}`::
751 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
752 sage: lyapunov_rank(K)
755 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
756 itself [Rudolf et al.]_::
758 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
759 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
764 The Lyapunov rank should be additive on a product of proper cones
767 sage: set_random_seed()
768 sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
769 sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
770 sage: K = K1.cartesian_product(K2)
771 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
774 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
775 itself [Rudolf et al.]_::
777 sage: set_random_seed()
778 sage: K = random_cone(max_dim=10, max_rays=10)
779 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
782 Make sure we exercise the non-strictly-convex/non-solid case::
784 sage: set_random_seed()
785 sage: K = random_cone(max_dim=10, strictly_convex=False, solid=False)
786 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
789 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
790 be any number between `1` and `n` inclusive, excluding `n-1`
791 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
792 trivial cone in a trivial space as well. However, in zero dimensions,
793 the Lyapunov rank of the trivial cone will be zero::
795 sage: set_random_seed()
796 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
797 sage: b = lyapunov_rank(K)
798 sage: n = K.lattice_dim()
799 sage: (n == 0 or 1 <= b) and b <= n
804 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
805 Lyapunov rank `n-1` in `n` dimensions::
807 sage: set_random_seed()
808 sage: K = random_cone(max_dim=10)
809 sage: b = lyapunov_rank(K)
810 sage: n = K.lattice_dim()
814 The calculation of the Lyapunov rank of an improper cone can be
815 reduced to that of a proper cone [Orlitzky/Gowda]_::
817 sage: set_random_seed()
818 sage: K = random_cone(max_dim=10)
819 sage: actual = lyapunov_rank(K)
820 sage: K_S = restrict_span(K)
821 sage: P = restrict_span(K_S.dual()).dual()
822 sage: l = lineality(K)
824 sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
825 sage: actual == expected
828 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
830 sage: set_random_seed()
831 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
832 sage: lyapunov_rank(K) == len(LL(K))
844 # K is not solid, project onto its span.
848 beta
+= m
*(n
- m
) + (n
- m
)**2
851 # K is not pointed, project its dual onto its span.
852 # Uses a proposition from our paper, i.e. this is
853 # equivalent to K = restrict_span(K.dual()).dual()
854 K
= restrict_span(intersect_span(K
,K
.dual()), K
.dual())
855 #K = restrict_span(K.dual()).dual()
857 #Ks = [ list(r) for r in sorted(K.rays()) ]
858 #Js = [ list(r) for r in sorted(J.rays()) ]
861 # print [ list(r) for r in K_orig.rays() ]