]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
3f5a4fed4e1c49853f00eafcf6084744223ca296
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()
40 The isomorphism should be an inner product space isomorphism, and
41 thus it should preserve dual cones (and commute with the "dual"
42 operation). But beware the automatic renaming of the dual lattice.
43 OH AND YOU HAVE TO SORT THE CONES::
45 sage: K = random_cone(max_dim=10, strictly_convex=False, solid=True)
47 sage: rename_lattice(L, 'L')
48 sage: (phi, phi_inv) = span_iso(K)
49 sage: sorted(phi_inv( phi(K).dual() )) == sorted(K.dual())
52 We may need to isomorph twice to make sure we stop moving down to
53 smaller spaces. (Once you've done this on a cone and its dual, the
54 result should be proper.) OH AND YOU HAVE TO SORT THE CONES::
56 sage: K = random_cone(max_dim=10, strictly_convex=False, solid=False)
58 sage: rename_lattice(L, 'L')
59 sage: (phi, phi_inv) = span_iso(K)
61 sage: (phi_dual, phi_dual_inv) = span_iso(K_S.dual())
62 sage: J_T = phi_dual(K_S.dual()).dual()
63 sage: phi_inv(phi_dual_inv(J_T)) == K
67 phi_domain
= K
.sublattice().vector_space()
68 phi_codo
= VectorSpace(phi_domain
.base_field(), phi_domain
.dimension())
70 # S goes from the new space to the cone space.
71 S
= linear_transformation(phi_codo
, phi_domain
, phi_domain
.basis())
73 # phi goes from the cone space to the new space.
76 Takes a cone ``J`` and sends it into the new space.
78 newrays
= map(S
.inverse(), J_orig
.rays())
83 return Cone(newrays
, lattice
=L
)
85 def phi_inverse(J_sub
):
87 The inverse to phi which goes from the new space to the cone space.
89 newrays
= map(S
, J_sub
.rays())
90 return Cone(newrays
, lattice
=K
.lattice())
93 return (phi
, phi_inverse
)
97 def discrete_complementarity_set(K
):
99 Compute the discrete complementarity set of this cone.
101 The complementarity set of this cone is the set of all orthogonal
102 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
103 dual. The discrete complementarity set restricts `x` and `s` to be
104 generators of their respective cones.
108 A list of pairs `(x,s)` such that,
110 * `x` is in this cone.
111 * `x` is a generator of this cone.
112 * `s` is in this cone's dual.
113 * `s` is a generator of this cone's dual.
114 * `x` and `s` are orthogonal.
118 The discrete complementarity set of the nonnegative orthant consists
119 of pairs of standard basis vectors::
121 sage: K = Cone([(1,0),(0,1)])
122 sage: discrete_complementarity_set(K)
123 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
125 If the cone consists of a single ray, the second components of the
126 discrete complementarity set should generate the orthogonal
127 complement of that ray::
129 sage: K = Cone([(1,0)])
130 sage: discrete_complementarity_set(K)
131 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
132 sage: K = Cone([(1,0,0)])
133 sage: discrete_complementarity_set(K)
134 [((1, 0, 0), (0, 1, 0)),
135 ((1, 0, 0), (0, -1, 0)),
136 ((1, 0, 0), (0, 0, 1)),
137 ((1, 0, 0), (0, 0, -1))]
139 When the cone is the entire space, its dual is the trivial cone, so
140 the discrete complementarity set is empty::
142 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
143 sage: discrete_complementarity_set(K)
148 The complementarity set of the dual can be obtained by switching the
149 components of the complementarity set of the original cone::
151 sage: K1 = random_cone(max_dim=10, max_rays=10)
153 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
154 sage: actual = discrete_complementarity_set(K1)
155 sage: actual == expected
159 V
= K
.lattice().vector_space()
161 # Convert the rays to vectors so that we can compute inner
163 xs
= [V(x
) for x
in K
.rays()]
164 ss
= [V(s
) for s
in K
.dual().rays()]
166 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
171 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
176 A list of matrices forming a basis for the space of all
177 Lyapunov-like transformations on the given cone.
181 The trivial cone has no Lyapunov-like transformations::
183 sage: L = ToricLattice(0)
184 sage: K = Cone([], lattice=L)
188 The Lyapunov-like transformations on the nonnegative orthant are
189 simply diagonal matrices::
191 sage: K = Cone([(1,)])
195 sage: K = Cone([(1,0),(0,1)])
202 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
205 [1 0 0] [0 0 0] [0 0 0]
206 [0 0 0] [0 1 0] [0 0 0]
207 [0 0 0], [0 0 0], [0 0 1]
210 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
211 `L^{3}_{\infty}` cones [Rudolf et al.]_::
213 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
221 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
231 The inner product `\left< L\left(x\right), s \right>` is zero for
232 every pair `\left( x,s \right)` in the discrete complementarity set
235 sage: K = random_cone(max_dim=8, max_rays=10)
236 sage: C_of_K = discrete_complementarity_set(K)
237 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
238 sage: sum(map(abs, l))
242 V
= K
.lattice().vector_space()
244 C_of_K
= discrete_complementarity_set(K
)
246 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
248 # Sage doesn't think matrices are vectors, so we have to convert
249 # our matrices to vectors explicitly before we can figure out how
250 # many are linearly-indepenedent.
252 # The space W has the same base ring as V, but dimension
253 # dim(V)^2. So it has the same dimension as the space of linear
254 # transformations on V. In other words, it's just the right size
255 # to create an isomorphism between it and our matrices.
256 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
258 # Turn our matrices into long vectors...
259 vectors
= [ W(m
.list()) for m
in tensor_products
]
261 # Vector space representation of Lyapunov-like matrices
262 # (i.e. vec(L) where L is Luapunov-like).
263 LL_vector
= W
.span(vectors
).complement()
265 # Now construct an ambient MatrixSpace in which to stick our
267 M
= MatrixSpace(V
.base_ring(), V
.dimension())
269 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
275 def lyapunov_rank(K
):
277 Compute the Lyapunov (or bilinearity) rank of this cone.
279 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
281 1. The dimension of the Lie algebra of the automorphism group of the
284 2. The dimension of the linear space of all Lyapunov-like
285 transformations on the cone.
289 A closed, convex polyhedral cone.
293 An integer representing the Lyapunov rank of the cone. If the
294 dimension of the ambient vector space is `n`, then the Lyapunov rank
295 will be between `1` and `n` inclusive; however a rank of `n-1` is
296 not possible (see the first reference).
300 In the references, the cones are always assumed to be proper. We
301 do not impose this restriction.
309 The codimension formula from the second reference is used. We find
310 all pairs `(x,s)` in the complementarity set of `K` such that `x`
311 and `s` are rays of our cone. It is known that these vectors are
312 sufficient to apply the codimension formula. Once we have all such
313 pairs, we "brute force" the codimension formula by finding all
314 linearly-independent `xs^{T}`.
318 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
319 cone and Lyapunov-like transformations, Mathematical Programming, 147
322 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
323 Improper Cone. Work in-progress.
325 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
326 optimality constraints for the cone of positive polynomials,
327 Mathematical Programming, Series B, 129 (2011) 5-31.
331 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
334 sage: positives = Cone([(1,)])
335 sage: lyapunov_rank(positives)
337 sage: quadrant = Cone([(1,0), (0,1)])
338 sage: lyapunov_rank(quadrant)
340 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
341 sage: lyapunov_rank(octant)
344 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
347 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
348 sage: lyapunov_rank(L31)
351 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
353 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
354 sage: lyapunov_rank(L3infty)
357 The Lyapunov rank should be additive on a product of cones
360 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
361 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
362 sage: K = L31.cartesian_product(octant)
363 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
366 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
367 The cone ``K`` in the following example is isomorphic to the nonnegative
368 octant in `\mathbb{R}^{3}`::
370 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
371 sage: lyapunov_rank(K)
374 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
375 itself [Rudolf et al.]_::
377 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
378 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
383 The Lyapunov rank should be additive on a product of cones
386 sage: K1 = random_cone(max_dim=10, max_rays=10)
387 sage: K2 = random_cone(max_dim=10, max_rays=10)
388 sage: K = K1.cartesian_product(K2)
389 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
392 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
393 itself [Rudolf et al.]_::
395 sage: K = random_cone(max_dim=10, max_rays=10)
396 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
399 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
400 be any number between `1` and `n` inclusive, excluding `n-1`
401 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
402 trivial cone in a trivial space as well. However, in zero dimensions,
403 the Lyapunov rank of the trivial cone will be zero::
405 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
406 sage: b = lyapunov_rank(K)
407 sage: n = K.lattice_dim()
408 sage: (n == 0 or 1 <= b) and b <= n
413 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
414 Lyapunov rank `n-1` in `n` dimensions::
416 sage: K = random_cone(max_dim=10)
417 sage: b = lyapunov_rank(K)
418 sage: n = K.lattice_dim()
422 The calculation of the Lyapunov rank of an improper cone can be
423 reduced to that of a proper cone [Orlitzky/Gowda]_::
425 sage: K = random_cone(max_dim=15, solid=False, strictly_convex=False)
426 sage: actual = lyapunov_rank(K)
427 sage: (phi1, phi1_inv) = span_iso(K)
429 sage: (phi2, phi2_inv) = span_iso(K_S.dual())
430 sage: J_T = phi2(K_S.dual()).dual()
431 sage: phi1_inv(phi2_inv(J_T)) == K
433 sage: l = K.linear_subspace().dimension()
434 sage: codim = K.lattice_dim() - K.dim()
435 sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
436 sage: actual == expected