]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
821b95958e81f1e9a93dcf176f774a4c82e9317f
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 Test whether or not ``K1`` and ``K2`` are "basically the same."
15 This is a hack to get around the fact that it's difficult to tell
16 when two cones are linearly isomorphic. We have a proposition that
17 equates two cones, but represented over `\mathbb{Q}`, they are
18 merely linearly isomorphic (not equal). So rather than test for
19 equality, we test a list of properties that should be preserved
20 under an invertible linear transformation.
24 ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
29 Any proper cone with three generators in `\mathbb{R}^{3}` is
30 basically the same as the nonnegative orthant::
32 sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
33 sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
34 sage: _basically_the_same(K1, K2)
37 Negating a cone gives you another cone that is basically the same::
39 sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
40 sage: _basically_the_same(K, -K)
45 Any cone is basically the same as itself::
47 sage: K = random_cone(max_ambient_dim = 8)
48 sage: _basically_the_same(K, K)
51 After applying an invertible matrix to the rows of a cone, the
52 result should be basically the same as the cone we started with::
54 sage: K1 = random_cone(max_ambient_dim = 8)
55 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
56 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
57 sage: _basically_the_same(K1, K2)
61 if K1
.lattice_dim() != K2
.lattice_dim():
64 if K1
.nrays() != K2
.nrays():
67 if K1
.dim() != K2
.dim():
70 if K1
.lineality() != K2
.lineality():
73 if K1
.is_solid() != K2
.is_solid():
76 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
79 if len(LL(K1
)) != len(LL(K2
)):
82 C_of_K1
= discrete_complementarity_set(K1
)
83 C_of_K2
= discrete_complementarity_set(K2
)
84 if len(C_of_K1
) != len(C_of_K2
):
87 if len(K1
.facets()) != len(K2
.facets()):
96 Restrict ``K`` into its own span, or the span of another cone.
100 - ``K2`` -- another cone whose lattice has the same rank as this
105 A new cone in a sublattice.
109 sage: K = Cone([(1,)])
113 sage: K2 = Cone([(1,0)])
114 sage: _rho(K2).rays()
117 sage: K3 = Cone([(1,0,0)])
118 sage: _rho(K3).rays()
121 sage: _rho(K2) == _rho(K3)
126 The projected cone should always be solid::
128 sage: set_random_seed()
129 sage: K = random_cone(max_ambient_dim = 8)
134 And the resulting cone should live in a space having the same
135 dimension as the space we restricted it to::
137 sage: set_random_seed()
138 sage: K = random_cone(max_ambient_dim = 8)
139 sage: K_S = _rho(K, K.dual() )
140 sage: K_S.lattice_dim() == K.dual().dim()
143 This function should not affect the dimension of a cone::
145 sage: set_random_seed()
146 sage: K = random_cone(max_ambient_dim = 8)
147 sage: K.dim() == _rho(K).dim()
150 Nor should it affect the lineality of a cone::
152 sage: set_random_seed()
153 sage: K = random_cone(max_ambient_dim = 8)
154 sage: K.lineality() == _rho(K).lineality()
157 No matter which space we restrict to, the lineality should not
160 sage: set_random_seed()
161 sage: K = random_cone(max_ambient_dim = 8)
162 sage: K.lineality() >= _rho(K).lineality()
164 sage: K.lineality() >= _rho(K, K.dual()).lineality()
167 If we do this according to our paper, then the result is proper::
169 sage: set_random_seed()
170 sage: K = random_cone(max_ambient_dim = 8)
172 sage: K_SP = _rho(K_S.dual()).dual()
173 sage: K_SP.is_proper()
175 sage: K_SP = _rho(K_S, K_S.dual())
176 sage: K_SP.is_proper()
179 Test the proposition in our paper concerning the duals and
180 restrictions. Generate a random cone, then create a subcone of
181 it. The operation of dual-taking should then commute with rho::
183 sage: set_random_seed()
184 sage: J = random_cone(max_ambient_dim = 8)
185 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
186 sage: K_W_star = _rho(K, J).dual()
187 sage: K_star_W = _rho(K.dual(), J)
188 sage: _basically_the_same(K_W_star, K_star_W)
195 # First we project K onto the span of K2. This will explode if the
196 # rank of ``K2.lattice()`` doesn't match ours.
197 span_K2
= Cone(K2
.rays() + (-K2
).rays(), lattice
=K
.lattice())
198 K
= K
.intersection(span_K2
)
200 # Cheat a little to get the subspace span(K2). The paper uses the
201 # rays of K2 as a basis, but everything is invariant under linear
202 # isomorphism (i.e. a change of basis), and this is a little
204 W
= span_K2
.linear_subspace()
206 # We've already intersected K with the span of K2, so every
207 # generator of K should belong to W now.
208 W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
210 L
= ToricLattice(K2
.dim())
211 return Cone(W_rays
, lattice
=L
)
215 def discrete_complementarity_set(K
):
217 Compute the discrete complementarity set of this cone.
219 The complementarity set of a cone is the set of all orthogonal pairs
220 `(x,s)` such that `x` is in the cone, and `s` is in its dual. The
221 discrete complementarity set is a subset of the complementarity set
222 where `x` and `s` are required to be generators of their respective
225 For polyhedral cones, the discrete complementarity set is always
230 A list of pairs `(x,s)` such that,
232 * Both `x` and `s` are vectors (not rays).
233 * `x` is a generator of this cone.
234 * `s` is a generator of this cone's dual.
235 * `x` and `s` are orthogonal.
239 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
240 Improper Cone. Work in-progress.
244 The discrete complementarity set of the nonnegative orthant consists
245 of pairs of standard basis vectors::
247 sage: K = Cone([(1,0),(0,1)])
248 sage: discrete_complementarity_set(K)
249 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
251 If the cone consists of a single ray, the second components of the
252 discrete complementarity set should generate the orthogonal
253 complement of that ray::
255 sage: K = Cone([(1,0)])
256 sage: discrete_complementarity_set(K)
257 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
258 sage: K = Cone([(1,0,0)])
259 sage: discrete_complementarity_set(K)
260 [((1, 0, 0), (0, 1, 0)),
261 ((1, 0, 0), (0, -1, 0)),
262 ((1, 0, 0), (0, 0, 1)),
263 ((1, 0, 0), (0, 0, -1))]
265 When the cone is the entire space, its dual is the trivial cone, so
266 the discrete complementarity set is empty::
268 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
269 sage: discrete_complementarity_set(K)
272 Likewise when this cone is trivial (its dual is the entire space)::
274 sage: L = ToricLattice(0)
275 sage: K = Cone([], ToricLattice(0))
276 sage: discrete_complementarity_set(K)
281 The complementarity set of the dual can be obtained by switching the
282 components of the complementarity set of the original cone::
284 sage: set_random_seed()
285 sage: K1 = random_cone(max_ambient_dim=6)
287 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
288 sage: actual = discrete_complementarity_set(K1)
289 sage: sorted(actual) == sorted(expected)
292 The pairs in the discrete complementarity set are in fact
295 sage: set_random_seed()
296 sage: K = random_cone(max_ambient_dim=6)
297 sage: dcs = discrete_complementarity_set(K)
298 sage: sum([x.inner_product(s).abs() for (x,s) in dcs])
302 V
= K
.lattice().vector_space()
304 # Convert rays to vectors so that we can compute inner products.
305 xs
= [V(x
) for x
in K
.rays()]
307 # We also convert the generators of the dual cone so that we
308 # return pairs of vectors and not (vector, ray) pairs.
309 ss
= [V(s
) for s
in K
.dual().rays()]
311 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
316 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
321 A list of matrices forming a basis for the space of all
322 Lyapunov-like transformations on the given cone.
326 The trivial cone has no Lyapunov-like transformations::
328 sage: L = ToricLattice(0)
329 sage: K = Cone([], lattice=L)
333 The Lyapunov-like transformations on the nonnegative orthant are
334 simply diagonal matrices::
336 sage: K = Cone([(1,)])
340 sage: K = Cone([(1,0),(0,1)])
347 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
350 [1 0 0] [0 0 0] [0 0 0]
351 [0 0 0] [0 1 0] [0 0 0]
352 [0 0 0], [0 0 0], [0 0 1]
355 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
356 `L^{3}_{\infty}` cones [Rudolf et al.]_::
358 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
366 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
374 If our cone is the entire space, then every transformation on it is
377 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
378 sage: M = MatrixSpace(QQ,2)
379 sage: M.basis() == LL(K)
384 The inner product `\left< L\left(x\right), s \right>` is zero for
385 every pair `\left( x,s \right)` in the discrete complementarity set
388 sage: set_random_seed()
389 sage: K = random_cone(max_ambient_dim=8)
390 sage: C_of_K = discrete_complementarity_set(K)
391 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
392 sage: sum(map(abs, l))
395 The Lyapunov-like transformations on a cone and its dual are related
396 by transposition, but we're not guaranteed to compute transposed
397 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
400 sage: set_random_seed()
401 sage: K = random_cone(max_ambient_dim=8)
402 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
403 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
404 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
405 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
406 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
410 V
= K
.lattice().vector_space()
412 C_of_K
= discrete_complementarity_set(K
)
414 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
416 # Sage doesn't think matrices are vectors, so we have to convert
417 # our matrices to vectors explicitly before we can figure out how
418 # many are linearly-indepenedent.
420 # The space W has the same base ring as V, but dimension
421 # dim(V)^2. So it has the same dimension as the space of linear
422 # transformations on V. In other words, it's just the right size
423 # to create an isomorphism between it and our matrices.
424 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
426 # Turn our matrices into long vectors...
427 vectors
= [ W(m
.list()) for m
in tensor_products
]
429 # Vector space representation of Lyapunov-like matrices
430 # (i.e. vec(L) where L is Luapunov-like).
431 LL_vector
= W
.span(vectors
).complement()
433 # Now construct an ambient MatrixSpace in which to stick our
435 M
= MatrixSpace(V
.base_ring(), V
.dimension())
437 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
443 def lyapunov_rank(K
):
445 Compute the Lyapunov rank (or bilinearity rank) of this cone.
447 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
449 1. The dimension of the Lie algebra of the automorphism group of the
452 2. The dimension of the linear space of all Lyapunov-like
453 transformations on the cone.
457 A closed, convex polyhedral cone.
461 An integer representing the Lyapunov rank of the cone. If the
462 dimension of the ambient vector space is `n`, then the Lyapunov rank
463 will be between `1` and `n` inclusive; however a rank of `n-1` is
464 not possible (see [Orlitzky/Gowda]_).
468 The codimension formula from the second reference is used. We find
469 all pairs `(x,s)` in the complementarity set of `K` such that `x`
470 and `s` are rays of our cone. It is known that these vectors are
471 sufficient to apply the codimension formula. Once we have all such
472 pairs, we "brute force" the codimension formula by finding all
473 linearly-independent `xs^{T}`.
477 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
478 cone and Lyapunov-like transformations, Mathematical Programming, 147
481 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
482 Improper Cone. Work in-progress.
484 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
485 optimality constraints for the cone of positive polynomials,
486 Mathematical Programming, Series B, 129 (2011) 5-31.
490 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
493 sage: positives = Cone([(1,)])
494 sage: lyapunov_rank(positives)
496 sage: quadrant = Cone([(1,0), (0,1)])
497 sage: lyapunov_rank(quadrant)
499 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
500 sage: lyapunov_rank(octant)
503 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
506 sage: R5 = VectorSpace(QQ, 5)
507 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
509 sage: lyapunov_rank(K)
512 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
515 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
516 sage: lyapunov_rank(L31)
519 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
521 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
522 sage: lyapunov_rank(L3infty)
525 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
526 + 1` [Orlitzky/Gowda]_::
528 sage: K = Cone([(1,0,0,0,0)])
529 sage: lyapunov_rank(K)
531 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
534 A subspace (of dimension `m`) in `n` dimensions should have a
535 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
537 sage: e1 = (1,0,0,0,0)
538 sage: neg_e1 = (-1,0,0,0,0)
539 sage: e2 = (0,1,0,0,0)
540 sage: neg_e2 = (0,-1,0,0,0)
541 sage: z = (0,0,0,0,0)
542 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
543 sage: lyapunov_rank(K)
545 sage: K.lattice_dim()**2 - K.dim()*K.codim()
548 The Lyapunov rank should be additive on a product of proper cones
551 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
552 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
553 sage: K = L31.cartesian_product(octant)
554 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
557 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
558 The cone ``K`` in the following example is isomorphic to the nonnegative
559 octant in `\mathbb{R}^{3}`::
561 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
562 sage: lyapunov_rank(K)
565 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
566 itself [Rudolf et al.]_::
568 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
569 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
574 The Lyapunov rank should be additive on a product of proper cones
577 sage: set_random_seed()
578 sage: K1 = random_cone(max_ambient_dim=8,
579 ....: strictly_convex=True,
581 sage: K2 = random_cone(max_ambient_dim=8,
582 ....: strictly_convex=True,
584 sage: K = K1.cartesian_product(K2)
585 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
588 The Lyapunov rank is invariant under a linear isomorphism
591 sage: K1 = random_cone(max_ambient_dim = 8)
592 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
593 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
594 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
597 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
598 itself [Rudolf et al.]_::
600 sage: set_random_seed()
601 sage: K = random_cone(max_ambient_dim=8)
602 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
605 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
606 be any number between `1` and `n` inclusive, excluding `n-1`
607 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
608 trivial cone in a trivial space as well. However, in zero dimensions,
609 the Lyapunov rank of the trivial cone will be zero::
611 sage: set_random_seed()
612 sage: K = random_cone(max_ambient_dim=8,
613 ....: strictly_convex=True,
615 sage: b = lyapunov_rank(K)
616 sage: n = K.lattice_dim()
617 sage: (n == 0 or 1 <= b) and b <= n
622 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
623 Lyapunov rank `n-1` in `n` dimensions::
625 sage: set_random_seed()
626 sage: K = random_cone(max_ambient_dim=8)
627 sage: b = lyapunov_rank(K)
628 sage: n = K.lattice_dim()
632 The calculation of the Lyapunov rank of an improper cone can be
633 reduced to that of a proper cone [Orlitzky/Gowda]_::
635 sage: set_random_seed()
636 sage: K = random_cone(max_ambient_dim=8)
637 sage: actual = lyapunov_rank(K)
639 sage: K_SP = _rho(K_S.dual()).dual()
640 sage: l = K.lineality()
642 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
643 sage: actual == expected
646 The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
648 sage: set_random_seed()
649 sage: K = random_cone(max_ambient_dim=8)
650 sage: lyapunov_rank(K) == len(LL(K))
653 We can make an imperfect cone perfect by adding a slack variable
654 (a Theorem in [Orlitzky/Gowda]_)::
656 sage: set_random_seed()
657 sage: K = random_cone(max_ambient_dim=8,
658 ....: strictly_convex=True,
660 sage: L = ToricLattice(K.lattice_dim() + 1)
661 sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
662 sage: lyapunov_rank(K) >= K.lattice_dim()
673 # K is not solid, restrict to its span.
676 # Non-solid reduction lemma.
680 # K is not pointed, restrict to the span of its dual. Uses a
681 # proposition from our paper, i.e. this is equivalent to K =
682 # _rho(K.dual()).dual().
683 K
= _rho(K
, K
.dual())
685 # Non-pointed reduction lemma.