]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
60f9c34ec8bc271d65812859f51ca77636c8cbbc
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('../../'))
12 Project ``K`` into its own span.
16 sage: K = Cone([(1,)])
17 sage: project_span(K) == K
20 sage: K2 = Cone([(1,0)])
21 sage: project_span(K2).rays()
24 sage: K3 = Cone([(1,0,0)])
25 sage: project_span(K3).rays()
28 sage: project_span(K2) == project_span(K3)
33 The projected cone should always be solid::
35 sage: K = random_cone()
36 sage: K_S = project_span(K)
40 If we do this according to our paper, then the result is proper::
42 sage: K = random_cone()
43 sage: K_S = project_span(K)
44 sage: P = project_span(K_S.dual()).dual()
49 F
= K
.lattice().base_field()
50 Q
= K
.lattice().quotient(K
.sublattice_complement())
51 vecs
= [ vector(F
, reversed(list(Q(r
)))) for r
in K
.rays() ]
57 return Cone(vecs
, lattice
=L
)
60 def rename_lattice(L
,s
):
62 Change all names of the given lattice to ``s``.
67 L
._latex
_dual
_name
= s
71 Return an isomorphism (and its inverse) that will send ``K`` into a
72 lower-dimensional space isomorphic to its span (and back).
76 The inverse composed with the isomorphism should be the identity::
78 sage: K = random_cone(max_dim=10)
79 sage: (phi, phi_inv) = span_iso(K)
80 sage: phi_inv(phi(K)) == K
83 The image of ``K`` under the isomorphism should have full dimension::
85 sage: K = random_cone(max_dim=10)
86 sage: (phi, phi_inv) = span_iso(K)
87 sage: phi(K).dim() == phi(K).lattice_dim()
91 phi_domain
= K
.sublattice().vector_space()
92 phi_codo
= VectorSpace(phi_domain
.base_field(), phi_domain
.dimension())
94 # S goes from the new space to the cone space.
95 S
= linear_transformation(phi_codo
, phi_domain
, phi_domain
.basis())
97 # phi goes from the cone space to the new space.
100 Takes a cone ``J`` and sends it into the new space.
102 newrays
= map(S
.inverse(), J_orig
.rays())
104 if len(newrays
) == 0:
107 return Cone(newrays
, lattice
=L
)
109 def phi_inverse(J_sub
):
111 The inverse to phi which goes from the new space to the cone space.
113 newrays
= map(S
, J_sub
.rays())
114 return Cone(newrays
, lattice
=K
.lattice())
117 return (phi
, phi_inverse
)
121 def discrete_complementarity_set(K
):
123 Compute the discrete complementarity set of this cone.
125 The complementarity set of this cone is the set of all orthogonal
126 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
127 dual. The discrete complementarity set restricts `x` and `s` to be
128 generators of their respective cones.
132 A list of pairs `(x,s)` such that,
134 * `x` is in this cone.
135 * `x` is a generator of this cone.
136 * `s` is in this cone's dual.
137 * `s` is a generator of this cone's dual.
138 * `x` and `s` are orthogonal.
142 The discrete complementarity set of the nonnegative orthant consists
143 of pairs of standard basis vectors::
145 sage: K = Cone([(1,0),(0,1)])
146 sage: discrete_complementarity_set(K)
147 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
149 If the cone consists of a single ray, the second components of the
150 discrete complementarity set should generate the orthogonal
151 complement of that ray::
153 sage: K = Cone([(1,0)])
154 sage: discrete_complementarity_set(K)
155 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
156 sage: K = Cone([(1,0,0)])
157 sage: discrete_complementarity_set(K)
158 [((1, 0, 0), (0, 1, 0)),
159 ((1, 0, 0), (0, -1, 0)),
160 ((1, 0, 0), (0, 0, 1)),
161 ((1, 0, 0), (0, 0, -1))]
163 When the cone is the entire space, its dual is the trivial cone, so
164 the discrete complementarity set is empty::
166 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
167 sage: discrete_complementarity_set(K)
172 The complementarity set of the dual can be obtained by switching the
173 components of the complementarity set of the original cone::
175 sage: K1 = random_cone(max_dim=10, max_rays=10)
177 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
178 sage: actual = discrete_complementarity_set(K1)
179 sage: actual == expected
183 V
= K
.lattice().vector_space()
185 # Convert the rays to vectors so that we can compute inner
187 xs
= [V(x
) for x
in K
.rays()]
188 ss
= [V(s
) for s
in K
.dual().rays()]
190 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
195 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
200 A list of matrices forming a basis for the space of all
201 Lyapunov-like transformations on the given cone.
205 The trivial cone has no Lyapunov-like transformations::
207 sage: L = ToricLattice(0)
208 sage: K = Cone([], lattice=L)
212 The Lyapunov-like transformations on the nonnegative orthant are
213 simply diagonal matrices::
215 sage: K = Cone([(1,)])
219 sage: K = Cone([(1,0),(0,1)])
226 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
229 [1 0 0] [0 0 0] [0 0 0]
230 [0 0 0] [0 1 0] [0 0 0]
231 [0 0 0], [0 0 0], [0 0 1]
234 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
235 `L^{3}_{\infty}` cones [Rudolf et al.]_::
237 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
245 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
255 The inner product `\left< L\left(x\right), s \right>` is zero for
256 every pair `\left( x,s \right)` in the discrete complementarity set
259 sage: K = random_cone(max_dim=8, max_rays=10)
260 sage: C_of_K = discrete_complementarity_set(K)
261 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
262 sage: sum(map(abs, l))
266 V
= K
.lattice().vector_space()
268 C_of_K
= discrete_complementarity_set(K
)
270 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
272 # Sage doesn't think matrices are vectors, so we have to convert
273 # our matrices to vectors explicitly before we can figure out how
274 # many are linearly-indepenedent.
276 # The space W has the same base ring as V, but dimension
277 # dim(V)^2. So it has the same dimension as the space of linear
278 # transformations on V. In other words, it's just the right size
279 # to create an isomorphism between it and our matrices.
280 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
282 # Turn our matrices into long vectors...
283 vectors
= [ W(m
.list()) for m
in tensor_products
]
285 # Vector space representation of Lyapunov-like matrices
286 # (i.e. vec(L) where L is Luapunov-like).
287 LL_vector
= W
.span(vectors
).complement()
289 # Now construct an ambient MatrixSpace in which to stick our
291 M
= MatrixSpace(V
.base_ring(), V
.dimension())
293 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
299 def lyapunov_rank(K
):
301 Compute the Lyapunov (or bilinearity) rank of this cone.
303 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
305 1. The dimension of the Lie algebra of the automorphism group of the
308 2. The dimension of the linear space of all Lyapunov-like
309 transformations on the cone.
313 A closed, convex polyhedral cone.
317 An integer representing the Lyapunov rank of the cone. If the
318 dimension of the ambient vector space is `n`, then the Lyapunov rank
319 will be between `1` and `n` inclusive; however a rank of `n-1` is
320 not possible (see the first reference).
324 In the references, the cones are always assumed to be proper. We
325 do not impose this restriction.
333 The codimension formula from the second reference is used. We find
334 all pairs `(x,s)` in the complementarity set of `K` such that `x`
335 and `s` are rays of our cone. It is known that these vectors are
336 sufficient to apply the codimension formula. Once we have all such
337 pairs, we "brute force" the codimension formula by finding all
338 linearly-independent `xs^{T}`.
342 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
343 cone and Lyapunov-like transformations, Mathematical Programming, 147
346 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
347 Improper Cone. Work in-progress.
349 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
350 optimality constraints for the cone of positive polynomials,
351 Mathematical Programming, Series B, 129 (2011) 5-31.
355 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
358 sage: positives = Cone([(1,)])
359 sage: lyapunov_rank(positives)
361 sage: quadrant = Cone([(1,0), (0,1)])
362 sage: lyapunov_rank(quadrant)
364 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
365 sage: lyapunov_rank(octant)
368 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
371 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
372 sage: lyapunov_rank(L31)
375 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
377 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
378 sage: lyapunov_rank(L3infty)
381 The Lyapunov rank should be additive on a product of cones
384 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
385 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
386 sage: K = L31.cartesian_product(octant)
387 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
390 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
391 The cone ``K`` in the following example is isomorphic to the nonnegative
392 octant in `\mathbb{R}^{3}`::
394 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
395 sage: lyapunov_rank(K)
398 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
399 itself [Rudolf et al.]_::
401 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
402 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
407 The Lyapunov rank should be additive on a product of cones
410 sage: K1 = random_cone(max_dim=10, max_rays=10)
411 sage: K2 = random_cone(max_dim=10, max_rays=10)
412 sage: K = K1.cartesian_product(K2)
413 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
416 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
417 itself [Rudolf et al.]_::
419 sage: K = random_cone(max_dim=10, max_rays=10)
420 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
423 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
424 be any number between `1` and `n` inclusive, excluding `n-1`
425 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
426 trivial cone in a trivial space as well. However, in zero dimensions,
427 the Lyapunov rank of the trivial cone will be zero::
429 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
430 sage: b = lyapunov_rank(K)
431 sage: n = K.lattice_dim()
432 sage: (n == 0 or 1 <= b) and b <= n
437 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
438 Lyapunov rank `n-1` in `n` dimensions::
440 sage: K = random_cone(max_dim=10)
441 sage: b = lyapunov_rank(K)
442 sage: n = K.lattice_dim()
446 The calculation of the Lyapunov rank of an improper cone can be
447 reduced to that of a proper cone [Orlitzky/Gowda]_::
449 sage: K = random_cone(max_dim=15, solid=False, strictly_convex=False)
450 sage: actual = lyapunov_rank(K)
451 sage: (phi1, _) = span_iso(K)
453 sage: (phi2, _) = span_iso(K_S.dual())
454 sage: J_T = phi2(K_S.dual()).dual()
455 sage: l = K.linear_subspace().dimension()
456 sage: codim = K.lattice_dim() - K.dim()
457 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
458 sage: actual == expected
461 Repeat the previous test with different ``random_cone()`` params::
463 sage: K = random_cone(max_dim=15, solid=False, strictly_convex=True)
464 sage: actual = lyapunov_rank(K)
465 sage: (phi1, _) = span_iso(K)
467 sage: (phi2, _) = span_iso(K_S.dual())
468 sage: J_T = phi2(K_S.dual()).dual()
469 sage: l = K.linear_subspace().dimension()
470 sage: codim = K.lattice_dim() - K.dim()
471 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
472 sage: actual == expected
475 sage: K = random_cone(max_dim=15, solid=True, strictly_convex=False)
476 sage: actual = lyapunov_rank(K)
477 sage: (phi1, _) = span_iso(K)
479 sage: (phi2, _) = span_iso(K_S.dual())
480 sage: J_T = phi2(K_S.dual()).dual()
481 sage: l = K.linear_subspace().dimension()
482 sage: codim = K.lattice_dim() - K.dim()
483 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
484 sage: actual == expected
487 sage: K = random_cone(max_dim=15, solid=True, strictly_convex=True)
488 sage: actual = lyapunov_rank(K)
489 sage: (phi1, _) = span_iso(K)
491 sage: (phi2, _) = span_iso(K_S.dual())
492 sage: J_T = phi2(K_S.dual()).dual()
493 sage: l = K.linear_subspace().dimension()
494 sage: codim = K.lattice_dim() - K.dim()
495 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
496 sage: actual == expected
499 sage: K = random_cone(max_dim=15)
500 sage: actual = lyapunov_rank(K)
501 sage: (phi1, _) = span_iso(K)
503 sage: (phi2, _) = span_iso(K_S.dual())
504 sage: J_T = phi2(K_S.dual()).dual()
505 sage: l = K.linear_subspace().dimension()
506 sage: codim = K.lattice_dim() - K.dim()
507 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
508 sage: actual == expected
511 And test with the project_span function::
513 sage: K = random_cone(max_dim=15)
514 sage: actual = lyapunov_rank(K)
515 sage: K_S = project_span(K)
516 sage: P = project_span(K_S.dual()).dual()
517 sage: l = K.linear_subspace().dimension()
518 sage: codim = K.lattice_dim() - K.dim()
519 sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
520 sage: actual == expected