]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
55d9a06d55313c94f71be29783171c96233ccc31
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 Project ``K`` into its own span.
17 sage: K = Cone([(1,)])
18 sage: project_span(K) == K
21 sage: K2 = Cone([(1,0)])
22 sage: project_span(K2).rays()
25 sage: K3 = Cone([(1,0,0)])
26 sage: project_span(K3).rays()
29 sage: project_span(K2) == project_span(K3)
34 The projected cone should always be solid::
36 sage: K = random_cone(max_dim = 10)
37 sage: K_S = project_span(K)
41 If we do this according to our paper, then the result is proper::
43 sage: K = random_cone(max_dim = 10)
44 sage: K_S = project_span(K)
45 sage: P = project_span(K_S.dual()).dual()
52 Q
= L
.quotient(K
.sublattice_complement())
53 vecs
= [ vector(F
, reversed(list(Q(r
)))) for r
in K
.rays() ]
57 newL
= ToricLattice(0)
59 return Cone(vecs
, lattice
=newL
)
65 Compute the lineality of this cone.
67 The lineality of a cone is the dimension of the largest linear
68 subspace contained in that cone.
72 A nonnegative integer; the dimension of the largest subspace
73 contained within this cone.
77 .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
78 University Press, Princeton, 1970.
82 The lineality of the nonnegative orthant is zero, since it clearly
85 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
89 However, if we add another ray so that the entire `x`-axis belongs
90 to the cone, then the resulting cone will have lineality one::
92 sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
96 If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
97 to the dimension of the ambient space (i.e. two)::
99 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
103 Per the definition, the lineality of the trivial cone in a trivial
106 sage: K = Cone([], lattice=ToricLattice(0))
112 The lineality of a cone should be an integer between zero and the
113 dimension of the ambient space, inclusive::
115 sage: K = random_cone(max_dim = 10)
116 sage: l = lineality(K)
119 sage: (0 <= l) and (l <= K.lattice_dim())
122 A strictly convex cone should have lineality zero::
124 sage: K = random_cone(max_dim = 10, strictly_convex = True)
129 return K
.linear_subspace().dimension()
134 Compute the codimension of this cone.
136 The codimension of a cone is the dimension of the space of all
137 elements perpendicular to every element of the cone. In other words,
138 the codimension is the difference between the dimension of the
139 ambient space and the dimension of the cone itself.
143 A nonnegative integer representing the dimension of the space of all
144 elements perpendicular to this cone.
148 :meth:`dim`, :meth:`lattice_dim`
152 The codimension of the nonnegative orthant is zero, since the span of
153 its generators equals the entire ambient space::
155 sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
159 However, if we remove a ray so that the entire cone is contained
160 within the `x-y`-plane, then the resulting cone will have
161 codimension one, because the `z`-axis is perpendicular to every
162 element of the cone::
164 sage: K = Cone([(1,0,0), (0,1,0)])
168 If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
170 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
174 And if the cone is trivial in any space, then its codimension is
175 equal to the dimension of the ambient space::
177 sage: K = Cone([], lattice=ToricLattice(0))
181 sage: K = Cone([(0,)])
185 sage: K = Cone([(0,0)])
191 The codimension of a cone should be an integer between zero and
192 the dimension of the ambient space, inclusive::
194 sage: K = random_cone(max_dim = 10)
198 sage: (0 <= c) and (c <= K.lattice_dim())
201 A solid cone should have codimension zero::
203 sage: K = random_cone(max_dim = 10, solid = True)
207 The codimension of a cone is equal to the lineality of its dual::
209 sage: K = random_cone(max_dim = 10, solid = True)
210 sage: codim(K) == lineality(K.dual())
214 return (K
.lattice_dim() - K
.dim())
217 def discrete_complementarity_set(K
):
219 Compute the discrete complementarity set of this cone.
221 The complementarity set of this cone is the set of all orthogonal
222 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
223 dual. The discrete complementarity set restricts `x` and `s` to be
224 generators of their respective cones.
228 A list of pairs `(x,s)` such that,
230 * `x` is in this cone.
231 * `x` is a generator of this cone.
232 * `s` is in this cone's dual.
233 * `s` is a generator of this cone's dual.
234 * `x` and `s` are orthogonal.
238 The discrete complementarity set of the nonnegative orthant consists
239 of pairs of standard basis vectors::
241 sage: K = Cone([(1,0),(0,1)])
242 sage: discrete_complementarity_set(K)
243 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
245 If the cone consists of a single ray, the second components of the
246 discrete complementarity set should generate the orthogonal
247 complement of that ray::
249 sage: K = Cone([(1,0)])
250 sage: discrete_complementarity_set(K)
251 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
252 sage: K = Cone([(1,0,0)])
253 sage: discrete_complementarity_set(K)
254 [((1, 0, 0), (0, 1, 0)),
255 ((1, 0, 0), (0, -1, 0)),
256 ((1, 0, 0), (0, 0, 1)),
257 ((1, 0, 0), (0, 0, -1))]
259 When the cone is the entire space, its dual is the trivial cone, so
260 the discrete complementarity set is empty::
262 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
263 sage: discrete_complementarity_set(K)
268 The complementarity set of the dual can be obtained by switching the
269 components of the complementarity set of the original cone::
271 sage: K1 = random_cone(max_dim=10, max_rays=10)
273 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
274 sage: actual = discrete_complementarity_set(K1)
275 sage: actual == expected
279 V
= K
.lattice().vector_space()
281 # Convert the rays to vectors so that we can compute inner
283 xs
= [V(x
) for x
in K
.rays()]
284 ss
= [V(s
) for s
in K
.dual().rays()]
286 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
291 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
296 A list of matrices forming a basis for the space of all
297 Lyapunov-like transformations on the given cone.
301 The trivial cone has no Lyapunov-like transformations::
303 sage: L = ToricLattice(0)
304 sage: K = Cone([], lattice=L)
308 The Lyapunov-like transformations on the nonnegative orthant are
309 simply diagonal matrices::
311 sage: K = Cone([(1,)])
315 sage: K = Cone([(1,0),(0,1)])
322 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
325 [1 0 0] [0 0 0] [0 0 0]
326 [0 0 0] [0 1 0] [0 0 0]
327 [0 0 0], [0 0 0], [0 0 1]
330 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
331 `L^{3}_{\infty}` cones [Rudolf et al.]_::
333 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
341 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
351 The inner product `\left< L\left(x\right), s \right>` is zero for
352 every pair `\left( x,s \right)` in the discrete complementarity set
355 sage: K = random_cone(max_dim=8, max_rays=10)
356 sage: C_of_K = discrete_complementarity_set(K)
357 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
358 sage: sum(map(abs, l))
362 V
= K
.lattice().vector_space()
364 C_of_K
= discrete_complementarity_set(K
)
366 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
368 # Sage doesn't think matrices are vectors, so we have to convert
369 # our matrices to vectors explicitly before we can figure out how
370 # many are linearly-indepenedent.
372 # The space W has the same base ring as V, but dimension
373 # dim(V)^2. So it has the same dimension as the space of linear
374 # transformations on V. In other words, it's just the right size
375 # to create an isomorphism between it and our matrices.
376 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
378 # Turn our matrices into long vectors...
379 vectors
= [ W(m
.list()) for m
in tensor_products
]
381 # Vector space representation of Lyapunov-like matrices
382 # (i.e. vec(L) where L is Luapunov-like).
383 LL_vector
= W
.span(vectors
).complement()
385 # Now construct an ambient MatrixSpace in which to stick our
387 M
= MatrixSpace(V
.base_ring(), V
.dimension())
389 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
395 def lyapunov_rank(K
):
397 Compute the Lyapunov (or bilinearity) rank of this cone.
399 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
401 1. The dimension of the Lie algebra of the automorphism group of the
404 2. The dimension of the linear space of all Lyapunov-like
405 transformations on the cone.
409 A closed, convex polyhedral cone.
413 An integer representing the Lyapunov rank of the cone. If the
414 dimension of the ambient vector space is `n`, then the Lyapunov rank
415 will be between `1` and `n` inclusive; however a rank of `n-1` is
416 not possible (see the first reference).
420 In the references, the cones are always assumed to be proper. We
421 do not impose this restriction.
429 The codimension formula from the second reference is used. We find
430 all pairs `(x,s)` in the complementarity set of `K` such that `x`
431 and `s` are rays of our cone. It is known that these vectors are
432 sufficient to apply the codimension formula. Once we have all such
433 pairs, we "brute force" the codimension formula by finding all
434 linearly-independent `xs^{T}`.
438 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
439 cone and Lyapunov-like transformations, Mathematical Programming, 147
442 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
443 Improper Cone. Work in-progress.
445 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
446 optimality constraints for the cone of positive polynomials,
447 Mathematical Programming, Series B, 129 (2011) 5-31.
451 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
454 sage: positives = Cone([(1,)])
455 sage: lyapunov_rank(positives)
457 sage: quadrant = Cone([(1,0), (0,1)])
458 sage: lyapunov_rank(quadrant)
460 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
461 sage: lyapunov_rank(octant)
464 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
467 sage: R5 = VectorSpace(QQ, 5)
468 sage: gens = R5.basis() + [ -r for r in R5.basis() ]
470 sage: lyapunov_rank(K)
473 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
476 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
477 sage: lyapunov_rank(L31)
480 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
482 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
483 sage: lyapunov_rank(L3infty)
486 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
487 + 1` [Orlitzky/Gowda]_::
489 sage: K = Cone([(1,0,0,0,0)])
490 sage: lyapunov_rank(K)
492 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
495 A subspace (of dimension `m`) in `n` dimensions should have a
496 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
498 sage: e1 = (1,0,0,0,0)
499 sage: neg_e1 = (-1,0,0,0,0)
500 sage: e2 = (0,1,0,0,0)
501 sage: neg_e2 = (0,-1,0,0,0)
502 sage: zero = (0,0,0,0,0)
503 sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
504 sage: lyapunov_rank(K)
506 sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
509 The Lyapunov rank should be additive on a product of proper cones
512 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
513 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
514 sage: K = L31.cartesian_product(octant)
515 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
518 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
519 The cone ``K`` in the following example is isomorphic to the nonnegative
520 octant in `\mathbb{R}^{3}`::
522 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
523 sage: lyapunov_rank(K)
526 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
527 itself [Rudolf et al.]_::
529 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
530 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
535 The Lyapunov rank should be additive on a product of proper cones
538 sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
539 sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
540 sage: K = K1.cartesian_product(K2)
541 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
544 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
545 itself [Rudolf et al.]_::
547 sage: K = random_cone(max_dim=10, max_rays=10)
548 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
551 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
552 be any number between `1` and `n` inclusive, excluding `n-1`
553 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
554 trivial cone in a trivial space as well. However, in zero dimensions,
555 the Lyapunov rank of the trivial cone will be zero::
557 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
558 sage: b = lyapunov_rank(K)
559 sage: n = K.lattice_dim()
560 sage: (n == 0 or 1 <= b) and b <= n
565 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
566 Lyapunov rank `n-1` in `n` dimensions::
568 sage: K = random_cone(max_dim=10)
569 sage: b = lyapunov_rank(K)
570 sage: n = K.lattice_dim()
574 The calculation of the Lyapunov rank of an improper cone can be
575 reduced to that of a proper cone [Orlitzky/Gowda]_::
577 sage: K = random_cone(max_dim=10)
578 sage: actual = lyapunov_rank(K)
579 sage: K_S = project_span(K)
580 sage: P = project_span(K_S.dual()).dual()
581 sage: l = lineality(K)
582 sage: codim = K.lattice_dim() - K.dim()
583 sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
584 sage: actual == expected
587 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
589 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
590 sage: lyapunov_rank(K) == len(LL(K))
601 # K is not solid, project onto its span.
605 beta
+= m
*(n
- m
) + (n
- m
)**2
608 # K is not pointed, project its dual onto its span.
609 K
= project_span(K
.dual()).dual()