]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
7597cd66ab5275ec914df0c4929152b2937bacf0
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('../../'))
10 def rename_lattice(L
,s
):
12 Change all names of the given lattice to ``s``.
17 L
._latex
_dual
_name
= s
21 Return an isomorphism (and its inverse) that will send ``K`` into a
22 lower-dimensional space isomorphic to its span (and back).
26 The inverse composed with the isomorphism should be the identity::
28 sage: K = random_cone(max_dim=10)
29 sage: (phi, phi_inv) = span_iso(K)
30 sage: phi_inv(phi(K)) == K
33 The image of ``K`` under the isomorphism should have full dimension::
35 sage: K = random_cone(max_dim=10)
36 sage: (phi, phi_inv) = span_iso(K)
37 sage: phi(K).dim() == phi(K).lattice_dim()
41 phi_domain
= K
.sublattice().vector_space()
42 phi_codo
= VectorSpace(phi_domain
.base_field(), phi_domain
.dimension())
44 # S goes from the new space to the cone space.
45 S
= linear_transformation(phi_codo
, phi_domain
, phi_domain
.basis())
47 # phi goes from the cone space to the new space.
50 Takes a cone ``J`` and sends it into the new space.
52 newrays
= map(S
.inverse(), J_orig
.rays())
57 return Cone(newrays
, lattice
=L
)
59 def phi_inverse(J_sub
):
61 The inverse to phi which goes from the new space to the cone space.
63 newrays
= map(S
, J_sub
.rays())
64 return Cone(newrays
, lattice
=K
.lattice())
67 return (phi
, phi_inverse
)
71 def discrete_complementarity_set(K
):
73 Compute the discrete complementarity set of this cone.
75 The complementarity set of this cone is the set of all orthogonal
76 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
77 dual. The discrete complementarity set restricts `x` and `s` to be
78 generators of their respective cones.
82 A list of pairs `(x,s)` such that,
84 * `x` is in this cone.
85 * `x` is a generator of this cone.
86 * `s` is in this cone's dual.
87 * `s` is a generator of this cone's dual.
88 * `x` and `s` are orthogonal.
92 The discrete complementarity set of the nonnegative orthant consists
93 of pairs of standard basis vectors::
95 sage: K = Cone([(1,0),(0,1)])
96 sage: discrete_complementarity_set(K)
97 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
99 If the cone consists of a single ray, the second components of the
100 discrete complementarity set should generate the orthogonal
101 complement of that ray::
103 sage: K = Cone([(1,0)])
104 sage: discrete_complementarity_set(K)
105 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
106 sage: K = Cone([(1,0,0)])
107 sage: discrete_complementarity_set(K)
108 [((1, 0, 0), (0, 1, 0)),
109 ((1, 0, 0), (0, -1, 0)),
110 ((1, 0, 0), (0, 0, 1)),
111 ((1, 0, 0), (0, 0, -1))]
113 When the cone is the entire space, its dual is the trivial cone, so
114 the discrete complementarity set is empty::
116 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
117 sage: discrete_complementarity_set(K)
122 The complementarity set of the dual can be obtained by switching the
123 components of the complementarity set of the original cone::
125 sage: K1 = random_cone(max_dim=10, max_rays=10)
127 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
128 sage: actual = discrete_complementarity_set(K1)
129 sage: actual == expected
133 V
= K
.lattice().vector_space()
135 # Convert the rays to vectors so that we can compute inner
137 xs
= [V(x
) for x
in K
.rays()]
138 ss
= [V(s
) for s
in K
.dual().rays()]
140 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
145 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
150 A list of matrices forming a basis for the space of all
151 Lyapunov-like transformations on the given cone.
155 The trivial cone has no Lyapunov-like transformations::
157 sage: L = ToricLattice(0)
158 sage: K = Cone([], lattice=L)
162 The Lyapunov-like transformations on the nonnegative orthant are
163 simply diagonal matrices::
165 sage: K = Cone([(1,)])
169 sage: K = Cone([(1,0),(0,1)])
176 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
179 [1 0 0] [0 0 0] [0 0 0]
180 [0 0 0] [0 1 0] [0 0 0]
181 [0 0 0], [0 0 0], [0 0 1]
184 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
185 `L^{3}_{\infty}` cones [Rudolf et al.]_::
187 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
195 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
205 The inner product `\left< L\left(x\right), s \right>` is zero for
206 every pair `\left( x,s \right)` in the discrete complementarity set
209 sage: K = random_cone(max_dim=8, max_rays=10)
210 sage: C_of_K = discrete_complementarity_set(K)
211 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
212 sage: sum(map(abs, l))
216 V
= K
.lattice().vector_space()
218 C_of_K
= discrete_complementarity_set(K
)
220 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
222 # Sage doesn't think matrices are vectors, so we have to convert
223 # our matrices to vectors explicitly before we can figure out how
224 # many are linearly-indepenedent.
226 # The space W has the same base ring as V, but dimension
227 # dim(V)^2. So it has the same dimension as the space of linear
228 # transformations on V. In other words, it's just the right size
229 # to create an isomorphism between it and our matrices.
230 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
232 # Turn our matrices into long vectors...
233 vectors
= [ W(m
.list()) for m
in tensor_products
]
235 # Vector space representation of Lyapunov-like matrices
236 # (i.e. vec(L) where L is Luapunov-like).
237 LL_vector
= W
.span(vectors
).complement()
239 # Now construct an ambient MatrixSpace in which to stick our
241 M
= MatrixSpace(V
.base_ring(), V
.dimension())
243 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
249 def lyapunov_rank(K
):
251 Compute the Lyapunov (or bilinearity) rank of this cone.
253 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
255 1. The dimension of the Lie algebra of the automorphism group of the
258 2. The dimension of the linear space of all Lyapunov-like
259 transformations on the cone.
263 A closed, convex polyhedral cone.
267 An integer representing the Lyapunov rank of the cone. If the
268 dimension of the ambient vector space is `n`, then the Lyapunov rank
269 will be between `1` and `n` inclusive; however a rank of `n-1` is
270 not possible (see the first reference).
274 In the references, the cones are always assumed to be proper. We
275 do not impose this restriction.
283 The codimension formula from the second reference is used. We find
284 all pairs `(x,s)` in the complementarity set of `K` such that `x`
285 and `s` are rays of our cone. It is known that these vectors are
286 sufficient to apply the codimension formula. Once we have all such
287 pairs, we "brute force" the codimension formula by finding all
288 linearly-independent `xs^{T}`.
292 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
293 cone and Lyapunov-like transformations, Mathematical Programming, 147
296 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
297 Improper Cone. Work in-progress.
299 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
300 optimality constraints for the cone of positive polynomials,
301 Mathematical Programming, Series B, 129 (2011) 5-31.
305 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
308 sage: positives = Cone([(1,)])
309 sage: lyapunov_rank(positives)
311 sage: quadrant = Cone([(1,0), (0,1)])
312 sage: lyapunov_rank(quadrant)
314 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
315 sage: lyapunov_rank(octant)
318 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
321 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
322 sage: lyapunov_rank(L31)
325 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
327 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
328 sage: lyapunov_rank(L3infty)
331 The Lyapunov rank should be additive on a product of cones
334 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
335 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
336 sage: K = L31.cartesian_product(octant)
337 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
340 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
341 The cone ``K`` in the following example is isomorphic to the nonnegative
342 octant in `\mathbb{R}^{3}`::
344 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
345 sage: lyapunov_rank(K)
348 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
349 itself [Rudolf et al.]_::
351 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
352 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
357 The Lyapunov rank should be additive on a product of cones
360 sage: K1 = random_cone(max_dim=10, max_rays=10)
361 sage: K2 = random_cone(max_dim=10, max_rays=10)
362 sage: K = K1.cartesian_product(K2)
363 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
366 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
367 itself [Rudolf et al.]_::
369 sage: K = random_cone(max_dim=10, max_rays=10)
370 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
373 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
374 be any number between `1` and `n` inclusive, excluding `n-1`
375 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
376 trivial cone in a trivial space as well. However, in zero dimensions,
377 the Lyapunov rank of the trivial cone will be zero::
379 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
380 sage: b = lyapunov_rank(K)
381 sage: n = K.lattice_dim()
382 sage: (n == 0 or 1 <= b) and b <= n
387 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
388 Lyapunov rank `n-1` in `n` dimensions::
390 sage: K = random_cone(max_dim=10)
391 sage: b = lyapunov_rank(K)
392 sage: n = K.lattice_dim()
396 The calculation of the Lyapunov rank of an improper cone can be
397 reduced to that of a proper cone [Orlitzky/Gowda]_::
399 sage: K = random_cone(max_dim=15, solid=False, strictly_convex=False)
400 sage: actual = lyapunov_rank(K)
401 sage: (phi1, _) = span_iso(K)
403 sage: (phi2, _) = span_iso(K_S.dual())
404 sage: J_T = phi2(K_S.dual()).dual()
405 sage: l = K.linear_subspace().dimension()
406 sage: codim = K.lattice_dim() - K.dim()
407 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
408 sage: actual == expected
411 Repeat the previous test with different ``random_cone()`` params::
413 sage: K = random_cone(max_dim=15, solid=False, strictly_convex=True)
414 sage: actual = lyapunov_rank(K)
415 sage: (phi1, _) = span_iso(K)
417 sage: (phi2, _) = span_iso(K_S.dual())
418 sage: J_T = phi2(K_S.dual()).dual()
419 sage: l = K.linear_subspace().dimension()
420 sage: codim = K.lattice_dim() - K.dim()
421 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
422 sage: actual == expected
425 sage: K = random_cone(max_dim=15, solid=True, strictly_convex=False)
426 sage: actual = lyapunov_rank(K)
427 sage: (phi1, _) = span_iso(K)
429 sage: (phi2, _) = span_iso(K_S.dual())
430 sage: J_T = phi2(K_S.dual()).dual()
431 sage: l = K.linear_subspace().dimension()
432 sage: codim = K.lattice_dim() - K.dim()
433 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
434 sage: actual == expected
437 sage: K = random_cone(max_dim=15, solid=True, strictly_convex=True)
438 sage: actual = lyapunov_rank(K)
439 sage: (phi1, _) = span_iso(K)
441 sage: (phi2, _) = span_iso(K_S.dual())
442 sage: J_T = phi2(K_S.dual()).dual()
443 sage: l = K.linear_subspace().dimension()
444 sage: codim = K.lattice_dim() - K.dim()
445 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
446 sage: actual == expected
449 sage: K = random_cone(max_dim=15)
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