]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
b4d2be0d73e10ed9f58d1189c84ed4bec3d0ad94
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()):
94 def _restrict_to_space(K
, W
):
96 Restrict this cone a subspace of its ambient space.
100 - ``W`` -- The subspace into which this cone will be restricted.
104 A new cone in a sublattice corresponding to ``W``.
108 When this cone is solid, restricting it into its own span should do
111 sage: K = Cone([(1,)])
112 sage: _restrict_to_space(K, K.span()) == K
115 A single ray restricted into its own span gives the same output
116 regardless of the ambient space::
118 sage: K2 = Cone([(1,0)])
119 sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
123 sage: K3 = Cone([(1,0,0)])
124 sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
133 The projected cone should always be solid::
135 sage: set_random_seed()
136 sage: K = random_cone(max_ambient_dim = 8)
137 sage: _restrict_to_space(K, K.span()).is_solid()
140 And the resulting cone should live in a space having the same
141 dimension as the space we restricted it to::
143 sage: set_random_seed()
144 sage: K = random_cone(max_ambient_dim = 8)
145 sage: K_P = _restrict_to_space(K, K.dual().span())
146 sage: K_P.lattice_dim() == K.dual().dim()
149 This function should not affect the dimension of a cone::
151 sage: set_random_seed()
152 sage: K = random_cone(max_ambient_dim = 8)
153 sage: K.dim() == _restrict_to_space(K,K.span()).dim()
156 Nor should it affect the lineality of a cone::
158 sage: set_random_seed()
159 sage: K = random_cone(max_ambient_dim = 8)
160 sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
163 No matter which space we restrict to, the lineality should not
166 sage: set_random_seed()
167 sage: K = random_cone(max_ambient_dim = 8)
168 sage: S = K.span(); P = K.dual().span()
169 sage: K.lineality() >= _restrict_to_space(K,S).lineality()
171 sage: K.lineality() >= _restrict_to_space(K,P).lineality()
174 If we do this according to our paper, then the result is proper::
176 sage: set_random_seed()
177 sage: K = random_cone(max_ambient_dim = 8)
178 sage: K_S = _restrict_to_space(K, K.span())
179 sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
180 sage: K_SP.is_proper()
182 sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
183 sage: K_SP.is_proper()
186 Test the proposition in our paper concerning the duals and
187 restrictions. Generate a random cone, then create a subcone of
188 it. The operation of dual-taking should then commute with
191 sage: set_random_seed()
192 sage: J = random_cone(max_ambient_dim = 8)
193 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
194 sage: K_W_star = _restrict_to_space(K, J.span()).dual()
195 sage: K_star_W = _restrict_to_space(K.dual(), J.span())
196 sage: _basically_the_same(K_W_star, K_star_W)
200 # First we want to intersect ``K`` with ``W``. The easiest way to
201 # do this is via cone intersection, so we turn the subspace ``W``
203 W_cone
= Cone(W
.basis() + [-b
for b
in W
.basis()], lattice
=K
.lattice())
204 K
= K
.intersection(W_cone
)
206 # We've already intersected K with the span of K2, so every
207 # generator of K should belong to W now.
208 K_W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
210 L
= ToricLattice(W
.dimension())
211 return Cone(K_W_rays
, lattice
=L
)
215 def discrete_complementarity_set(K
):
217 Compute a discrete complementarity set of this cone.
219 A discrete complementarity set of `K` is the set of all orthogonal
220 pairs `(x,s)` such that `x \in G_{1}` and `s \in G_{2}` for some
221 generating sets `G_{1}` of `K` and `G_{2}` of its dual. Polyhedral
222 convex cones are input in terms of their generators, so "the" (this
223 particular) discrete complementarity set corresponds to ``G1
224 == K.rays()`` and ``G2 == K.dual().rays()``.
228 A list of pairs `(x,s)` such that,
230 * Both `x` and `s` are vectors (not rays).
231 * `x` is one of ``K.rays()``.
232 * `s` is one of ``K.dual().rays()``.
233 * `x` and `s` are orthogonal.
237 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
238 Improper Cone. Work in-progress.
242 The discrete complementarity set of the nonnegative orthant consists
243 of pairs of standard basis vectors::
245 sage: K = Cone([(1,0),(0,1)])
246 sage: discrete_complementarity_set(K)
247 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
249 If the cone consists of a single ray, the second components of the
250 discrete complementarity set should generate the orthogonal
251 complement of that ray::
253 sage: K = Cone([(1,0)])
254 sage: discrete_complementarity_set(K)
255 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
256 sage: K = Cone([(1,0,0)])
257 sage: discrete_complementarity_set(K)
258 [((1, 0, 0), (0, 1, 0)),
259 ((1, 0, 0), (0, -1, 0)),
260 ((1, 0, 0), (0, 0, 1)),
261 ((1, 0, 0), (0, 0, -1))]
263 When the cone is the entire space, its dual is the trivial cone, so
264 the discrete complementarity set is empty::
266 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
267 sage: discrete_complementarity_set(K)
270 Likewise when this cone is trivial (its dual is the entire space)::
272 sage: L = ToricLattice(0)
273 sage: K = Cone([], ToricLattice(0))
274 sage: discrete_complementarity_set(K)
279 The complementarity set of the dual can be obtained by switching the
280 components of the complementarity set of the original cone::
282 sage: set_random_seed()
283 sage: K1 = random_cone(max_ambient_dim=6)
285 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
286 sage: actual = discrete_complementarity_set(K1)
287 sage: sorted(actual) == sorted(expected)
290 The pairs in the discrete complementarity set are in fact
293 sage: set_random_seed()
294 sage: K = random_cone(max_ambient_dim=6)
295 sage: dcs = discrete_complementarity_set(K)
296 sage: sum([x.inner_product(s).abs() for (x,s) in dcs])
300 V
= K
.lattice().vector_space()
302 # Convert rays to vectors so that we can compute inner products.
303 xs
= [V(x
) for x
in K
.rays()]
305 # We also convert the generators of the dual cone so that we
306 # return pairs of vectors and not (vector, ray) pairs.
307 ss
= [V(s
) for s
in K
.dual().rays()]
309 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
314 Compute a basis of Lyapunov-like transformations on this cone.
318 A list of matrices forming a basis for the space of all
319 Lyapunov-like transformations on the given cone.
323 The trivial cone has no Lyapunov-like transformations::
325 sage: L = ToricLattice(0)
326 sage: K = Cone([], lattice=L)
330 The Lyapunov-like transformations on the nonnegative orthant are
331 simply diagonal matrices::
333 sage: K = Cone([(1,)])
337 sage: K = Cone([(1,0),(0,1)])
344 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
347 [1 0 0] [0 0 0] [0 0 0]
348 [0 0 0] [0 1 0] [0 0 0]
349 [0 0 0], [0 0 0], [0 0 1]
352 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
353 `L^{3}_{\infty}` cones [Rudolf et al.]_::
355 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
363 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
371 If our cone is the entire space, then every transformation on it is
374 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
375 sage: M = MatrixSpace(QQ,2)
376 sage: M.basis() == LL(K)
381 The inner product `\left< L\left(x\right), s \right>` is zero for
382 every pair `\left( x,s \right)` in the discrete complementarity set
385 sage: set_random_seed()
386 sage: K = random_cone(max_ambient_dim=8)
387 sage: C_of_K = discrete_complementarity_set(K)
388 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
389 sage: sum(map(abs, l))
392 The Lyapunov-like transformations on a cone and its dual are related
393 by transposition, but we're not guaranteed to compute transposed
394 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
397 sage: set_random_seed()
398 sage: K = random_cone(max_ambient_dim=8)
399 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
400 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
401 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
402 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
403 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
407 V
= K
.lattice().vector_space()
409 C_of_K
= discrete_complementarity_set(K
)
411 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
413 # Sage doesn't think matrices are vectors, so we have to convert
414 # our matrices to vectors explicitly before we can figure out how
415 # many are linearly-indepenedent.
417 # The space W has the same base ring as V, but dimension
418 # dim(V)^2. So it has the same dimension as the space of linear
419 # transformations on V. In other words, it's just the right size
420 # to create an isomorphism between it and our matrices.
421 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
423 # Turn our matrices into long vectors...
424 vectors
= [ W(m
.list()) for m
in tensor_products
]
426 # Vector space representation of Lyapunov-like matrices
427 # (i.e. vec(L) where L is Luapunov-like).
428 LL_vector
= W
.span(vectors
).complement()
430 # Now construct an ambient MatrixSpace in which to stick our
432 M
= MatrixSpace(V
.base_ring(), V
.dimension())
434 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
440 def lyapunov_rank(K
):
442 Compute the Lyapunov rank (or bilinearity rank) of this cone.
444 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
446 1. The dimension of the Lie algebra of the automorphism group of the
449 2. The dimension of the linear space of all Lyapunov-like
450 transformations on the cone.
454 A closed, convex polyhedral cone.
458 An integer representing the Lyapunov rank of the cone. If the
459 dimension of the ambient vector space is `n`, then the Lyapunov rank
460 will be between `1` and `n` inclusive; however a rank of `n-1` is
461 not possible (see [Orlitzky/Gowda]_).
465 The codimension formula from the second reference is used. We find
466 all pairs `(x,s)` in the complementarity set of `K` such that `x`
467 and `s` are rays of our cone. It is known that these vectors are
468 sufficient to apply the codimension formula. Once we have all such
469 pairs, we "brute force" the codimension formula by finding all
470 linearly-independent `xs^{T}`.
474 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
475 cone and Lyapunov-like transformations, Mathematical Programming, 147
478 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
479 Improper Cone. Work in-progress.
481 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
482 optimality constraints for the cone of positive polynomials,
483 Mathematical Programming, Series B, 129 (2011) 5-31.
487 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
490 sage: positives = Cone([(1,)])
491 sage: lyapunov_rank(positives)
493 sage: quadrant = Cone([(1,0), (0,1)])
494 sage: lyapunov_rank(quadrant)
496 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
497 sage: lyapunov_rank(octant)
500 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
503 sage: R5 = VectorSpace(QQ, 5)
504 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
506 sage: lyapunov_rank(K)
509 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
512 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
513 sage: lyapunov_rank(L31)
516 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
518 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
519 sage: lyapunov_rank(L3infty)
522 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
523 + 1` [Orlitzky/Gowda]_::
525 sage: K = Cone([(1,0,0,0,0)])
526 sage: lyapunov_rank(K)
528 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
531 A subspace (of dimension `m`) in `n` dimensions should have a
532 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
534 sage: e1 = (1,0,0,0,0)
535 sage: neg_e1 = (-1,0,0,0,0)
536 sage: e2 = (0,1,0,0,0)
537 sage: neg_e2 = (0,-1,0,0,0)
538 sage: z = (0,0,0,0,0)
539 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
540 sage: lyapunov_rank(K)
542 sage: K.lattice_dim()**2 - K.dim()*K.codim()
545 The Lyapunov rank should be additive on a product of proper cones
548 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
549 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
550 sage: K = L31.cartesian_product(octant)
551 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
554 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
555 The cone ``K`` in the following example is isomorphic to the nonnegative
556 octant in `\mathbb{R}^{3}`::
558 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
559 sage: lyapunov_rank(K)
562 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
563 itself [Rudolf et al.]_::
565 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
566 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
571 The Lyapunov rank should be additive on a product of proper cones
574 sage: set_random_seed()
575 sage: K1 = random_cone(max_ambient_dim=8,
576 ....: strictly_convex=True,
578 sage: K2 = random_cone(max_ambient_dim=8,
579 ....: strictly_convex=True,
581 sage: K = K1.cartesian_product(K2)
582 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
585 The Lyapunov rank is invariant under a linear isomorphism
588 sage: K1 = random_cone(max_ambient_dim = 8)
589 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
590 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
591 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
594 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
595 itself [Rudolf et al.]_::
597 sage: set_random_seed()
598 sage: K = random_cone(max_ambient_dim=8)
599 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
602 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
603 be any number between `1` and `n` inclusive, excluding `n-1`
604 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
605 trivial cone in a trivial space as well. However, in zero dimensions,
606 the Lyapunov rank of the trivial cone will be zero::
608 sage: set_random_seed()
609 sage: K = random_cone(max_ambient_dim=8,
610 ....: strictly_convex=True,
612 sage: b = lyapunov_rank(K)
613 sage: n = K.lattice_dim()
614 sage: (n == 0 or 1 <= b) and b <= n
619 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
620 Lyapunov rank `n-1` in `n` dimensions::
622 sage: set_random_seed()
623 sage: K = random_cone(max_ambient_dim=8)
624 sage: b = lyapunov_rank(K)
625 sage: n = K.lattice_dim()
629 The calculation of the Lyapunov rank of an improper cone can be
630 reduced to that of a proper cone [Orlitzky/Gowda]_::
632 sage: set_random_seed()
633 sage: K = random_cone(max_ambient_dim=8)
634 sage: actual = lyapunov_rank(K)
635 sage: K_S = _restrict_to_space(K, K.span())
636 sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
637 sage: l = K.lineality()
639 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
640 sage: actual == expected
643 The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
645 sage: set_random_seed()
646 sage: K = random_cone(max_ambient_dim=8)
647 sage: lyapunov_rank(K) == len(LL(K))
650 We can make an imperfect cone perfect by adding a slack variable
651 (a Theorem in [Orlitzky/Gowda]_)::
653 sage: set_random_seed()
654 sage: K = random_cone(max_ambient_dim=8,
655 ....: strictly_convex=True,
657 sage: L = ToricLattice(K.lattice_dim() + 1)
658 sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
659 sage: lyapunov_rank(K) >= K.lattice_dim()
670 # K is not solid, restrict to its span.
671 K
= _restrict_to_space(K
, K
.span())
673 # Non-solid reduction lemma.
677 # K is not pointed, restrict to the span of its dual. Uses a
678 # proposition from our paper, i.e. this is equivalent to K =
679 # _rho(K.dual()).dual().
680 K
= _restrict_to_space(K
, K
.dual().span())
682 # Non-pointed reduction lemma.
690 def is_lyapunov_like(L
,K
):
692 Determine whether or not ``L`` is Lyapunov-like on ``K``.
694 We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
695 L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
696 `\left\langle x,s \right\rangle` in the complementarity set of
697 ``K``. It is known [Orlitzky]_ that this property need only be
698 checked for generators of ``K`` and its dual.
702 - ``L`` -- A linear transformation or matrix.
704 - ``K`` -- A polyhedral closed convex cone.
708 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
709 and ``False`` otherwise.
713 If this function returns ``True``, then ``L`` is Lyapunov-like
714 on ``K``. However, if ``False`` is returned, that could mean one
715 of two things. The first is that ``L`` is definitely not
716 Lyapunov-like on ``K``. The second is more of an "I don't know"
717 answer, returned (for example) if we cannot prove that an inner
722 .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
723 improper cone (preprint).
727 The identity is always Lyapunov-like in a nontrivial space::
729 sage: set_random_seed()
730 sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
731 sage: L = identity_matrix(K.lattice_dim())
732 sage: is_lyapunov_like(L,K)
735 As is the "zero" transformation::
737 sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
738 sage: R = K.lattice().vector_space().base_ring()
739 sage: L = zero_matrix(R, K.lattice_dim())
740 sage: is_lyapunov_like(L,K)
743 Everything in ``LL(K)`` should be Lyapunov-like on ``K``::
745 sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
746 sage: all([is_lyapunov_like(L,K) for L in LL(K)])
750 return all([(L
*x
).inner_product(s
) == 0
751 for (x
,s
) in discrete_complementarity_set(K
)])
754 def random_element(K
):
756 Return a random element of ``K`` from its ambient vector space.
760 The cone ``K`` is specified in terms of its generators, so that
761 ``K`` is equal to the convex conic combination of those generators.
762 To choose a random element of ``K``, we assign random nonnegative
763 coefficients to each generator of ``K`` and construct a new vector
764 from the scaled rays.
766 A vector, rather than a ray, is returned so that the element may
767 have non-integer coordinates. Thus the element may have an
768 arbitrarily small norm.
772 A random element of the trivial cone is zero::
774 sage: set_random_seed()
775 sage: K = Cone([], ToricLattice(0))
776 sage: random_element(K)
778 sage: K = Cone([(0,)])
779 sage: random_element(K)
781 sage: K = Cone([(0,0)])
782 sage: random_element(K)
784 sage: K = Cone([(0,0,0)])
785 sage: random_element(K)
790 Any cone should contain an element of itself::
792 sage: set_random_seed()
793 sage: K = random_cone(max_rays = 8)
794 sage: K.contains(random_element(K))
798 V
= K
.lattice().vector_space()
800 coefficients
= [ F
.random_element().abs() for i
in range(K
.nrays()) ]
801 vector_gens
= map(V
, K
.rays())
802 scaled_gens
= [ coefficients
[i
]*vector_gens
[i
]
803 for i
in range(len(vector_gens
)) ]
805 # Make sure we return a vector. Without the coercion, we might
806 # return ``0`` when ``K`` has no rays.
807 v
= V(sum(scaled_gens
))
811 def positive_operators(K
):
813 Compute generators of the cone of positive operators on this cone.
817 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
818 Each matrix ``P`` in the list should have the property that ``P*x``
819 is an element of ``K`` whenever ``x`` is an element of
820 ``K``. Moreover, any nonnegative linear combination of these
821 matrices shares the same property.
825 The trivial cone in a trivial space has no positive operators::
827 sage: K = Cone([], ToricLattice(0))
828 sage: positive_operators(K)
831 Positive operators on the nonnegative orthant are nonnegative matrices::
833 sage: K = Cone([(1,)])
834 sage: positive_operators(K)
837 sage: K = Cone([(1,0),(0,1)])
838 sage: positive_operators(K)
840 [1 0] [0 1] [0 0] [0 0]
841 [0 0], [0 0], [1 0], [0 1]
844 Every operator is positive on the ambient vector space::
846 sage: K = Cone([(1,),(-1,)])
847 sage: K.is_full_space()
849 sage: positive_operators(K)
852 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
853 sage: K.is_full_space()
855 sage: positive_operators(K)
857 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
858 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
863 A positive operator on a cone should send its generators into the cone::
865 sage: K = random_cone(max_ambient_dim = 6)
866 sage: pi_of_k = positive_operators(K)
867 sage: all([K.contains(p*x) for p in pi_of_k for x in K.rays()])
871 V
= K
.lattice().vector_space()
873 # Sage doesn't think matrices are vectors, so we have to convert
874 # our matrices to vectors explicitly before we can figure out how
875 # many are linearly-indepenedent.
877 # The space W has the same base ring as V, but dimension
878 # dim(V)^2. So it has the same dimension as the space of linear
879 # transformations on V. In other words, it's just the right size
880 # to create an isomorphism between it and our matrices.
881 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
883 G1
= [ V(x
) for x
in K
.rays() ]
884 G2
= [ V(s
) for s
in K
.dual().rays() ]
886 tensor_products
= [ s
.tensor_product(x
) for x
in G1
for s
in G2
]
888 # Turn our matrices into long vectors...
889 vectors
= [ W(m
.list()) for m
in tensor_products
]
891 # Create the *dual* cone of the positive operators, expressed as
893 L
= ToricLattice(W
.dimension())
894 pi_dual
= Cone(vectors
, lattice
=L
)
896 # Now compute the desired cone from its dual...
897 pi_cone
= pi_dual
.dual()
899 # And finally convert its rays back to matrix representations.
900 M
= MatrixSpace(V
.base_ring(), V
.dimension())
902 return [ M(v
.list()) for v
in pi_cone
.rays() ]