]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
f2e8b2e9ee104e6cbd216e7473fd65dd76947d98
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('../../'))
13 Construct the space `W \times W^{\perp}` isomorphic to the ambient space
14 of ``K`` where `W` is equal to the span of ``K``.
16 V
= K
.lattice().vector_space()
18 # Create the space W \times W^{\perp} isomorphic to V.
19 # First we get an orthogonal (but not normal) basis...
20 M
= matrix(V
.base_field(), K
.rays())
21 W_basis
,_
= M
.gram_schmidt()
23 W
= V
.subspace_with_basis(W_basis
)
24 W_perp
= W
.complement()
26 return W
.cartesian_product(W_perp
)
31 Construct the IPS isomorphism and its inverse from our paper.
33 Given a cone ``K``, the returned isomorphism will split its ambient
34 vector space `V` into a cartesian product `W \times W^{\perp}` where
35 `W` equals the span of ``K``.
37 V
= K
.lattice().vector_space()
39 (W
, W_perp
) = V_iso
.cartesian_factors()
41 # A space equivalent to V, but using our basis.
42 V_user
= V
.subspace_with_basis( W
.basis() + W_perp
.basis() )
45 # Write v in terms of our custom basis, where the first dim(W)
46 # coordinates are for the W-part of the basis.
47 cs
= V_user
.coordinates(v
)
49 w1
= sum([ V_user
.basis()[idx
]*cs
[idx
]
50 for idx
in range(0, W
.dimension()) ])
51 w2
= sum([ V_user
.basis()[idx
]*cs
[idx
]
52 for idx
in range(W
.dimension(), V
.dimension()) ])
54 return V_iso( (w1
, w2
) )
58 # Crash if the arguments are in the wrong spaces.
61 #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
62 #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
63 # for idx in range(0,n-m) ])
65 return sum( pair
.cartesian_factors() )
72 def unrestrict_span(K
, K2
=None):
76 _
,phi_inv
= ips_iso(K2
)
78 (W
, W_perp
) = V_iso
.cartesian_factors()
82 w
= sum([ r
[idx
]*W
.basis()[idx
] for idx
in range(0,len(r
)) ])
83 pair
= V_iso( (w
, W_perp
.zero()) )
84 rays
.append( phi_inv(pair
) )
86 L
= ToricLattice(W
.dimension() + W_perp
.dimension())
88 return Cone(rays
, lattice
=L
)
92 def intersect_span(K1
, K2
):
94 Return a new cone obtained by intersecting ``K1`` with the span of ``K2``.
98 if L
.rank() != K2
.lattice().rank():
99 raise ValueError('K1 and K2 must belong to lattices of the same rank.')
101 SL_gens
= list(K2
.rays())
102 span_K2_gens
= SL_gens
+ [ -g
for g
in SL_gens
]
104 # The lattices have the same rank (see above) so this should work.
105 span_K2
= Cone(span_K2_gens
, L
)
106 return K1
.intersection(span_K2
)
110 def restrict_span(K
, K2
=None):
112 Restrict ``K`` into its own span, or the span of another cone.
116 - ``K2`` -- another cone whose lattice has the same rank as this cone.
120 A new cone in a sublattice.
124 sage: K = Cone([(1,)])
125 sage: restrict_span(K) == K
128 sage: K2 = Cone([(1,0)])
129 sage: restrict_span(K2).rays()
132 sage: K3 = Cone([(1,0,0)])
133 sage: restrict_span(K3).rays()
136 sage: restrict_span(K2) == restrict_span(K3)
141 The projected cone should always be solid::
143 sage: set_random_seed()
144 sage: K = random_cone(max_dim = 10)
145 sage: K_S = restrict_span(K)
149 And the resulting cone should live in a space having the same
150 dimension as the space we restricted it to::
152 sage: set_random_seed()
153 sage: K = random_cone(max_dim = 10)
154 sage: K_S = restrict_span( intersect_span(K, K.dual()), K.dual() )
155 sage: K_S.lattice_dim() == K.dual().dim()
158 This function has ``unrestrict_span()`` as its inverse::
160 sage: set_random_seed()
161 sage: K = random_cone(max_dim = 10, solid=True)
162 sage: J = restrict_span(K)
163 sage: K == unrestrict_span(J,K)
166 This function should not affect the dimension of a cone::
168 sage: set_random_seed()
169 sage: K = random_cone(max_dim = 10)
170 sage: K.dim() == restrict_span(K).dim()
173 Nor should it affect the lineality of a cone::
175 sage: set_random_seed()
176 sage: K = random_cone(max_dim = 10)
177 sage: lineality(K) == lineality(restrict_span(K))
180 No matter which space we restrict to, the lineality should not
183 sage: set_random_seed()
184 sage: K = random_cone(max_dim = 10)
185 sage: J = intersect_span(K, K.dual())
186 sage: lineality(K) >= lineality(restrict_span(J, K.dual()))
189 If we do this according to our paper, then the result is proper::
191 sage: set_random_seed()
192 sage: K = random_cone(max_dim = 10)
193 sage: K_S = restrict_span(K)
194 sage: P = restrict_span(K_S.dual()).dual()
198 If ``K`` is strictly convex, then both ``K_W`` and
199 ``K_star_W.dual()`` should equal ``K`` (after we unrestrict)::
201 sage: set_random_seed()
202 sage: K = random_cone(max_dim = 10, strictly_convex=True)
203 sage: K_W = restrict_span(intersect_span(K,K.dual()), K.dual())
204 sage: K_star_W_star = restrict_span(K.dual()).dual()
205 sage: j1 = unrestrict_span(K_W, K.dual())
206 sage: j2 = unrestrict_span(K_star_W_star, K.dual())
211 sage: K; [ list(r) for r in K.rays() ]
213 Test the proposition in our paper concerning the duals, where the
214 subspace `W` is the span of `K^{*}`::
216 sage: set_random_seed()
217 sage: K = random_cone(max_dim = 10, solid=False, strictly_convex=False)
218 sage: K_W = restrict_span(intersect_span(K,K.dual()), K.dual())
219 sage: K_star_W_star = restrict_span(K.dual(), K.dual()).dual()
220 sage: K_W.nrays() == K_star_W_star.nrays()
222 sage: K_W.dim() == K_star_W_star.dim()
224 sage: lineality(K_W) == lineality(K_star_W_star)
226 sage: K_W.is_solid() == K_star_W_star.is_solid()
228 sage: K_W.is_strictly_convex() == K_star_W_star.is_strictly_convex()
236 (W
, W_perp
) = iso_space(K2
).cartesian_factors()
238 ray_pairs
= [ phi(r
) for r
in K
.rays() ]
240 if any([ w2
!= W_perp
.zero() for (_
, w2
) in ray_pairs
]):
241 msg
= 'Cone has nonzero components in W-perp!'
242 raise ValueError(msg
)
244 # Represent the cone in terms of a basis for W, i.e. with smaller
246 ws
= [ W
.coordinate_vector(w1
) for (w1
, _
) in ray_pairs
]
248 L
= ToricLattice(W
.dimension())
250 return Cone(ws
, lattice
=L
)
256 Compute the lineality of this cone.
258 The lineality of a cone is the dimension of the largest linear
259 subspace contained in that cone.
263 A nonnegative integer; the dimension of the largest subspace
264 contained within this cone.
268 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
269 University Press, Princeton, 1970.
273 The lineality of the nonnegative orthant is zero, since it clearly
276 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
280 However, if we add another ray so that the entire `x`-axis belongs
281 to the cone, then the resulting cone will have lineality one::
283 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
287 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
288 to the dimension of the ambient space (i.e. two)::
290 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
294 Per the definition, the lineality of the trivial cone in a trivial
297 sage: K = Cone([], lattice=ToricLattice(0))
303 The lineality of a cone should be an integer between zero and the
304 dimension of the ambient space, inclusive::
306 sage: set_random_seed()
307 sage: K = random_cone(max_dim = 10)
308 sage: l = lineality(K)
311 sage: (0 <= l) and (l <= K.lattice_dim())
314 A strictly convex cone should have lineality zero::
316 sage: set_random_seed()
317 sage: K = random_cone(max_dim = 10, strictly_convex = True)
322 return K
.linear_subspace().dimension()
327 Compute the codimension of this cone.
329 The codimension of a cone is the dimension of the space of all
330 elements perpendicular to every element of the cone. In other words,
331 the codimension is the difference between the dimension of the
332 ambient space and the dimension of the cone itself.
336 A nonnegative integer representing the dimension of the space of all
337 elements perpendicular to this cone.
341 :meth:`dim`, :meth:`lattice_dim`
345 The codimension of the nonnegative orthant is zero, since the span of
346 its generators equals the entire ambient space::
348 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
352 However, if we remove a ray so that the entire cone is contained
353 within the `x-y`-plane, then the resulting cone will have
354 codimension one, because the `z`-axis is perpendicular to every
355 element of the cone::
357 sage: K = Cone([(1,0,0), (0,1,0)])
361 If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
363 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
367 And if the cone is trivial in any space, then its codimension is
368 equal to the dimension of the ambient space::
370 sage: K = Cone([], lattice=ToricLattice(0))
374 sage: K = Cone([(0,)])
378 sage: K = Cone([(0,0)])
384 The codimension of a cone should be an integer between zero and
385 the dimension of the ambient space, inclusive::
387 sage: set_random_seed()
388 sage: K = random_cone(max_dim = 10)
392 sage: (0 <= c) and (c <= K.lattice_dim())
395 A solid cone should have codimension zero::
397 sage: set_random_seed()
398 sage: K = random_cone(max_dim = 10, solid = True)
402 The codimension of a cone is equal to the lineality of its dual::
404 sage: set_random_seed()
405 sage: K = random_cone(max_dim = 10, solid = True)
406 sage: codim(K) == lineality(K.dual())
410 return (K
.lattice_dim() - K
.dim())
413 def discrete_complementarity_set(K
):
415 Compute the discrete complementarity set of this cone.
417 The complementarity set of this cone is the set of all orthogonal
418 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
419 dual. The discrete complementarity set restricts `x` and `s` to be
420 generators of their respective cones.
424 A list of pairs `(x,s)` such that,
426 * `x` is in this cone.
427 * `x` is a generator of this cone.
428 * `s` is in this cone's dual.
429 * `s` is a generator of this cone's dual.
430 * `x` and `s` are orthogonal.
434 The discrete complementarity set of the nonnegative orthant consists
435 of pairs of standard basis vectors::
437 sage: K = Cone([(1,0),(0,1)])
438 sage: discrete_complementarity_set(K)
439 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
441 If the cone consists of a single ray, the second components of the
442 discrete complementarity set should generate the orthogonal
443 complement of that ray::
445 sage: K = Cone([(1,0)])
446 sage: discrete_complementarity_set(K)
447 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
448 sage: K = Cone([(1,0,0)])
449 sage: discrete_complementarity_set(K)
450 [((1, 0, 0), (0, 1, 0)),
451 ((1, 0, 0), (0, -1, 0)),
452 ((1, 0, 0), (0, 0, 1)),
453 ((1, 0, 0), (0, 0, -1))]
455 When the cone is the entire space, its dual is the trivial cone, so
456 the discrete complementarity set is empty::
458 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
459 sage: discrete_complementarity_set(K)
464 The complementarity set of the dual can be obtained by switching the
465 components of the complementarity set of the original cone::
467 sage: set_random_seed()
468 sage: K1 = random_cone(max_dim=6)
470 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
471 sage: actual = discrete_complementarity_set(K1)
472 sage: sorted(actual) == sorted(expected)
476 V
= K
.lattice().vector_space()
478 # Convert the rays to vectors so that we can compute inner
480 xs
= [V(x
) for x
in K
.rays()]
481 ss
= [V(s
) for s
in K
.dual().rays()]
483 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
488 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
493 A list of matrices forming a basis for the space of all
494 Lyapunov-like transformations on the given cone.
498 The trivial cone has no Lyapunov-like transformations::
500 sage: L = ToricLattice(0)
501 sage: K = Cone([], lattice=L)
505 The Lyapunov-like transformations on the nonnegative orthant are
506 simply diagonal matrices::
508 sage: K = Cone([(1,)])
512 sage: K = Cone([(1,0),(0,1)])
519 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
522 [1 0 0] [0 0 0] [0 0 0]
523 [0 0 0] [0 1 0] [0 0 0]
524 [0 0 0], [0 0 0], [0 0 1]
527 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
528 `L^{3}_{\infty}` cones [Rudolf et al.]_::
530 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
538 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
548 The inner product `\left< L\left(x\right), s \right>` is zero for
549 every pair `\left( x,s \right)` in the discrete complementarity set
552 sage: set_random_seed()
553 sage: K = random_cone(max_dim=8, max_rays=10)
554 sage: C_of_K = discrete_complementarity_set(K)
555 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
556 sage: sum(map(abs, l))
559 The Lyapunov-like transformations on a cone and its dual are related
560 by transposition, but we're not guaranteed to compute transposed
561 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
564 sage: set_random_seed()
565 sage: K = random_cone(max_dim=8, max_rays=10)
566 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
567 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
568 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
569 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
570 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
574 V
= K
.lattice().vector_space()
576 C_of_K
= discrete_complementarity_set(K
)
578 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
580 # Sage doesn't think matrices are vectors, so we have to convert
581 # our matrices to vectors explicitly before we can figure out how
582 # many are linearly-indepenedent.
584 # The space W has the same base ring as V, but dimension
585 # dim(V)^2. So it has the same dimension as the space of linear
586 # transformations on V. In other words, it's just the right size
587 # to create an isomorphism between it and our matrices.
588 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
590 # Turn our matrices into long vectors...
591 vectors
= [ W(m
.list()) for m
in tensor_products
]
593 # Vector space representation of Lyapunov-like matrices
594 # (i.e. vec(L) where L is Luapunov-like).
595 LL_vector
= W
.span(vectors
).complement()
597 # Now construct an ambient MatrixSpace in which to stick our
599 M
= MatrixSpace(V
.base_ring(), V
.dimension())
601 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
607 def lyapunov_rank(K
):
609 Compute the Lyapunov (or bilinearity) rank of this cone.
611 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
613 1. The dimension of the Lie algebra of the automorphism group of the
616 2. The dimension of the linear space of all Lyapunov-like
617 transformations on the cone.
621 A closed, convex polyhedral cone.
625 An integer representing the Lyapunov rank of the cone. If the
626 dimension of the ambient vector space is `n`, then the Lyapunov rank
627 will be between `1` and `n` inclusive; however a rank of `n-1` is
628 not possible (see the first reference).
632 In the references, the cones are always assumed to be proper. We
633 do not impose this restriction.
641 The codimension formula from the second reference is used. We find
642 all pairs `(x,s)` in the complementarity set of `K` such that `x`
643 and `s` are rays of our cone. It is known that these vectors are
644 sufficient to apply the codimension formula. Once we have all such
645 pairs, we "brute force" the codimension formula by finding all
646 linearly-independent `xs^{T}`.
650 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
651 cone and Lyapunov-like transformations, Mathematical Programming, 147
654 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
655 Improper Cone. Work in-progress.
657 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
658 optimality constraints for the cone of positive polynomials,
659 Mathematical Programming, Series B, 129 (2011) 5-31.
663 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
666 sage: positives = Cone([(1,)])
667 sage: lyapunov_rank(positives)
669 sage: quadrant = Cone([(1,0), (0,1)])
670 sage: lyapunov_rank(quadrant)
672 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
673 sage: lyapunov_rank(octant)
676 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
679 sage: R5 = VectorSpace(QQ, 5)
680 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
682 sage: lyapunov_rank(K)
685 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
688 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
689 sage: lyapunov_rank(L31)
692 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
694 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
695 sage: lyapunov_rank(L3infty)
698 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
699 + 1` [Orlitzky/Gowda]_::
701 sage: K = Cone([(1,0,0,0,0)])
702 sage: lyapunov_rank(K)
704 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
707 A subspace (of dimension `m`) in `n` dimensions should have a
708 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
710 sage: e1 = (1,0,0,0,0)
711 sage: neg_e1 = (-1,0,0,0,0)
712 sage: e2 = (0,1,0,0,0)
713 sage: neg_e2 = (0,-1,0,0,0)
714 sage: z = (0,0,0,0,0)
715 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
716 sage: lyapunov_rank(K)
718 sage: K.lattice_dim()**2 - K.dim()*codim(K)
721 The Lyapunov rank should be additive on a product of proper cones
724 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
725 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
726 sage: K = L31.cartesian_product(octant)
727 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
730 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
731 The cone ``K`` in the following example is isomorphic to the nonnegative
732 octant in `\mathbb{R}^{3}`::
734 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
735 sage: lyapunov_rank(K)
738 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
739 itself [Rudolf et al.]_::
741 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
742 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
747 The Lyapunov rank should be additive on a product of proper cones
750 sage: set_random_seed()
751 sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
752 sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
753 sage: K = K1.cartesian_product(K2)
754 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
757 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
758 itself [Rudolf et al.]_::
760 sage: set_random_seed()
761 sage: K = random_cone(max_dim=10, max_rays=10)
762 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
765 Make sure we exercise the non-strictly-convex/non-solid case::
767 sage: set_random_seed()
768 sage: K = random_cone(max_dim=10, strictly_convex=False, solid=False)
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_dim=10, strictly_convex=True, solid=True)
780 sage: b = lyapunov_rank(K)
781 sage: n = K.lattice_dim()
782 sage: (n == 0 or 1 <= b) and b <= n
787 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
788 Lyapunov rank `n-1` in `n` dimensions::
790 sage: set_random_seed()
791 sage: K = random_cone(max_dim=10)
792 sage: b = lyapunov_rank(K)
793 sage: n = K.lattice_dim()
797 The calculation of the Lyapunov rank of an improper cone can be
798 reduced to that of a proper cone [Orlitzky/Gowda]_::
800 sage: set_random_seed()
801 sage: K = random_cone(max_dim=10)
802 sage: actual = lyapunov_rank(K)
803 sage: K_S = restrict_span(K)
804 sage: P = restrict_span(K_S.dual()).dual()
805 sage: l = lineality(K)
807 sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
808 sage: actual == expected
811 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
813 sage: set_random_seed()
814 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
815 sage: lyapunov_rank(K) == len(LL(K))
827 # K is not solid, project onto its span.
831 beta
+= m
*(n
- m
) + (n
- m
)**2
834 # K is not pointed, project its dual onto its span.
835 # Uses a proposition from our paper, i.e. this is
836 # equivalent to K = restrict_span(K.dual()).dual()
837 K
= restrict_span(intersect_span(K
,K
.dual()), K
.dual())
838 #K = restrict_span(K.dual()).dual()
840 #Ks = [ list(r) for r in sorted(K.rays()) ]
841 #Js = [ list(r) for r in sorted(J.rays()) ]
844 # print [ list(r) for r in K_orig.rays() ]