]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
d77793550aa44f02ea7047e2049b852b40db936f
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 cone should have lineality zero::
124 sage: K = random_cone(max_dim = 10, strictly_convex = True)
129 return K
.linear_subspace().dimension()
132 def discrete_complementarity_set(K
):
134 Compute the discrete complementarity set of this cone.
136 The complementarity set of this cone is the set of all orthogonal
137 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
138 dual. The discrete complementarity set restricts `x` and `s` to be
139 generators of their respective cones.
143 A list of pairs `(x,s)` such that,
145 * `x` is in this cone.
146 * `x` is a generator of this cone.
147 * `s` is in this cone's dual.
148 * `s` is a generator of this cone's dual.
149 * `x` and `s` are orthogonal.
153 The discrete complementarity set of the nonnegative orthant consists
154 of pairs of standard basis vectors::
156 sage: K = Cone([(1,0),(0,1)])
157 sage: discrete_complementarity_set(K)
158 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
160 If the cone consists of a single ray, the second components of the
161 discrete complementarity set should generate the orthogonal
162 complement of that ray::
164 sage: K = Cone([(1,0)])
165 sage: discrete_complementarity_set(K)
166 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
167 sage: K = Cone([(1,0,0)])
168 sage: discrete_complementarity_set(K)
169 [((1, 0, 0), (0, 1, 0)),
170 ((1, 0, 0), (0, -1, 0)),
171 ((1, 0, 0), (0, 0, 1)),
172 ((1, 0, 0), (0, 0, -1))]
174 When the cone is the entire space, its dual is the trivial cone, so
175 the discrete complementarity set is empty::
177 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
178 sage: discrete_complementarity_set(K)
183 The complementarity set of the dual can be obtained by switching the
184 components of the complementarity set of the original cone::
186 sage: K1 = random_cone(max_dim=10, max_rays=10)
188 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
189 sage: actual = discrete_complementarity_set(K1)
190 sage: actual == expected
194 V
= K
.lattice().vector_space()
196 # Convert the rays to vectors so that we can compute inner
198 xs
= [V(x
) for x
in K
.rays()]
199 ss
= [V(s
) for s
in K
.dual().rays()]
201 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
206 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
211 A list of matrices forming a basis for the space of all
212 Lyapunov-like transformations on the given cone.
216 The trivial cone has no Lyapunov-like transformations::
218 sage: L = ToricLattice(0)
219 sage: K = Cone([], lattice=L)
223 The Lyapunov-like transformations on the nonnegative orthant are
224 simply diagonal matrices::
226 sage: K = Cone([(1,)])
230 sage: K = Cone([(1,0),(0,1)])
237 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
240 [1 0 0] [0 0 0] [0 0 0]
241 [0 0 0] [0 1 0] [0 0 0]
242 [0 0 0], [0 0 0], [0 0 1]
245 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
246 `L^{3}_{\infty}` cones [Rudolf et al.]_::
248 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
256 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
266 The inner product `\left< L\left(x\right), s \right>` is zero for
267 every pair `\left( x,s \right)` in the discrete complementarity set
270 sage: K = random_cone(max_dim=8, max_rays=10)
271 sage: C_of_K = discrete_complementarity_set(K)
272 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
273 sage: sum(map(abs, l))
277 V
= K
.lattice().vector_space()
279 C_of_K
= discrete_complementarity_set(K
)
281 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
283 # Sage doesn't think matrices are vectors, so we have to convert
284 # our matrices to vectors explicitly before we can figure out how
285 # many are linearly-indepenedent.
287 # The space W has the same base ring as V, but dimension
288 # dim(V)^2. So it has the same dimension as the space of linear
289 # transformations on V. In other words, it's just the right size
290 # to create an isomorphism between it and our matrices.
291 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
293 # Turn our matrices into long vectors...
294 vectors
= [ W(m
.list()) for m
in tensor_products
]
296 # Vector space representation of Lyapunov-like matrices
297 # (i.e. vec(L) where L is Luapunov-like).
298 LL_vector
= W
.span(vectors
).complement()
300 # Now construct an ambient MatrixSpace in which to stick our
302 M
= MatrixSpace(V
.base_ring(), V
.dimension())
304 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
310 def lyapunov_rank(K
):
312 Compute the Lyapunov (or bilinearity) rank of this cone.
314 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
316 1. The dimension of the Lie algebra of the automorphism group of the
319 2. The dimension of the linear space of all Lyapunov-like
320 transformations on the cone.
324 A closed, convex polyhedral cone.
328 An integer representing the Lyapunov rank of the cone. If the
329 dimension of the ambient vector space is `n`, then the Lyapunov rank
330 will be between `1` and `n` inclusive; however a rank of `n-1` is
331 not possible (see the first reference).
335 In the references, the cones are always assumed to be proper. We
336 do not impose this restriction.
344 The codimension formula from the second reference is used. We find
345 all pairs `(x,s)` in the complementarity set of `K` such that `x`
346 and `s` are rays of our cone. It is known that these vectors are
347 sufficient to apply the codimension formula. Once we have all such
348 pairs, we "brute force" the codimension formula by finding all
349 linearly-independent `xs^{T}`.
353 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
354 cone and Lyapunov-like transformations, Mathematical Programming, 147
357 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
358 Improper Cone. Work in-progress.
360 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
361 optimality constraints for the cone of positive polynomials,
362 Mathematical Programming, Series B, 129 (2011) 5-31.
366 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
369 sage: positives = Cone([(1,)])
370 sage: lyapunov_rank(positives)
372 sage: quadrant = Cone([(1,0), (0,1)])
373 sage: lyapunov_rank(quadrant)
375 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
376 sage: lyapunov_rank(octant)
379 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
382 sage: R5 = VectorSpace(QQ, 5)
383 sage: gens = R5.basis() + [ -r for r in R5.basis() ]
385 sage: lyapunov_rank(K)
388 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
391 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
392 sage: lyapunov_rank(L31)
395 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
397 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
398 sage: lyapunov_rank(L3infty)
401 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
402 + 1` [Orlitzky/Gowda]_::
404 sage: K = Cone([(1,0,0,0,0)])
405 sage: lyapunov_rank(K)
407 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
410 A subspace (of dimension `m`) in `n` dimensions should have a
411 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
413 sage: e1 = (1,0,0,0,0)
414 sage: neg_e1 = (-1,0,0,0,0)
415 sage: e2 = (0,1,0,0,0)
416 sage: neg_e2 = (0,-1,0,0,0)
417 sage: zero = (0,0,0,0,0)
418 sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
419 sage: lyapunov_rank(K)
421 sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
424 The Lyapunov rank should be additive on a product of proper cones
427 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
428 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
429 sage: K = L31.cartesian_product(octant)
430 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
433 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
434 The cone ``K`` in the following example is isomorphic to the nonnegative
435 octant in `\mathbb{R}^{3}`::
437 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
438 sage: lyapunov_rank(K)
441 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
442 itself [Rudolf et al.]_::
444 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
445 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
450 The Lyapunov rank should be additive on a product of proper cones
453 sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
454 sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
455 sage: K = K1.cartesian_product(K2)
456 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
459 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
460 itself [Rudolf et al.]_::
462 sage: K = random_cone(max_dim=10, max_rays=10)
463 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
466 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
467 be any number between `1` and `n` inclusive, excluding `n-1`
468 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
469 trivial cone in a trivial space as well. However, in zero dimensions,
470 the Lyapunov rank of the trivial cone will be zero::
472 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
473 sage: b = lyapunov_rank(K)
474 sage: n = K.lattice_dim()
475 sage: (n == 0 or 1 <= b) and b <= n
480 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
481 Lyapunov rank `n-1` in `n` dimensions::
483 sage: K = random_cone(max_dim=10)
484 sage: b = lyapunov_rank(K)
485 sage: n = K.lattice_dim()
489 The calculation of the Lyapunov rank of an improper cone can be
490 reduced to that of a proper cone [Orlitzky/Gowda]_::
492 sage: K = random_cone(max_dim=10)
493 sage: actual = lyapunov_rank(K)
494 sage: K_S = project_span(K)
495 sage: P = project_span(K_S.dual()).dual()
496 sage: l = K.linear_subspace().dimension()
497 sage: codim = K.lattice_dim() - K.dim()
498 sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
499 sage: actual == expected
502 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
504 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
505 sage: lyapunov_rank(K) == len(LL(K))
513 l
= K
.linear_subspace().dimension()
516 # K is not solid, project onto its span.
520 beta
+= m
*(n
- m
) + (n
- m
)**2
523 # K is not pointed, project its dual onto its span.
524 K
= project_span(K
.dual()).dual()