]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
421fb3c8f04abb23fd420271ed10e0a7e65815d4
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
)
63 def discrete_complementarity_set(K
):
65 Compute the discrete complementarity set of this cone.
67 The complementarity set of this cone is the set of all orthogonal
68 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
69 dual. The discrete complementarity set restricts `x` and `s` to be
70 generators of their respective cones.
74 A list of pairs `(x,s)` such that,
76 * `x` is in this cone.
77 * `x` is a generator of this cone.
78 * `s` is in this cone's dual.
79 * `s` is a generator of this cone's dual.
80 * `x` and `s` are orthogonal.
84 The discrete complementarity set of the nonnegative orthant consists
85 of pairs of standard basis vectors::
87 sage: K = Cone([(1,0),(0,1)])
88 sage: discrete_complementarity_set(K)
89 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
91 If the cone consists of a single ray, the second components of the
92 discrete complementarity set should generate the orthogonal
93 complement of that ray::
95 sage: K = Cone([(1,0)])
96 sage: discrete_complementarity_set(K)
97 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
98 sage: K = Cone([(1,0,0)])
99 sage: discrete_complementarity_set(K)
100 [((1, 0, 0), (0, 1, 0)),
101 ((1, 0, 0), (0, -1, 0)),
102 ((1, 0, 0), (0, 0, 1)),
103 ((1, 0, 0), (0, 0, -1))]
105 When the cone is the entire space, its dual is the trivial cone, so
106 the discrete complementarity set is empty::
108 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
109 sage: discrete_complementarity_set(K)
114 The complementarity set of the dual can be obtained by switching the
115 components of the complementarity set of the original cone::
117 sage: K1 = random_cone(max_dim=10, max_rays=10)
119 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
120 sage: actual = discrete_complementarity_set(K1)
121 sage: actual == expected
125 V
= K
.lattice().vector_space()
127 # Convert the rays to vectors so that we can compute inner
129 xs
= [V(x
) for x
in K
.rays()]
130 ss
= [V(s
) for s
in K
.dual().rays()]
132 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
137 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
142 A list of matrices forming a basis for the space of all
143 Lyapunov-like transformations on the given cone.
147 The trivial cone has no Lyapunov-like transformations::
149 sage: L = ToricLattice(0)
150 sage: K = Cone([], lattice=L)
154 The Lyapunov-like transformations on the nonnegative orthant are
155 simply diagonal matrices::
157 sage: K = Cone([(1,)])
161 sage: K = Cone([(1,0),(0,1)])
168 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
171 [1 0 0] [0 0 0] [0 0 0]
172 [0 0 0] [0 1 0] [0 0 0]
173 [0 0 0], [0 0 0], [0 0 1]
176 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
177 `L^{3}_{\infty}` cones [Rudolf et al.]_::
179 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
187 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
197 The inner product `\left< L\left(x\right), s \right>` is zero for
198 every pair `\left( x,s \right)` in the discrete complementarity set
201 sage: K = random_cone(max_dim=8, max_rays=10)
202 sage: C_of_K = discrete_complementarity_set(K)
203 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
204 sage: sum(map(abs, l))
208 V
= K
.lattice().vector_space()
210 C_of_K
= discrete_complementarity_set(K
)
212 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
214 # Sage doesn't think matrices are vectors, so we have to convert
215 # our matrices to vectors explicitly before we can figure out how
216 # many are linearly-indepenedent.
218 # The space W has the same base ring as V, but dimension
219 # dim(V)^2. So it has the same dimension as the space of linear
220 # transformations on V. In other words, it's just the right size
221 # to create an isomorphism between it and our matrices.
222 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
224 # Turn our matrices into long vectors...
225 vectors
= [ W(m
.list()) for m
in tensor_products
]
227 # Vector space representation of Lyapunov-like matrices
228 # (i.e. vec(L) where L is Luapunov-like).
229 LL_vector
= W
.span(vectors
).complement()
231 # Now construct an ambient MatrixSpace in which to stick our
233 M
= MatrixSpace(V
.base_ring(), V
.dimension())
235 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
241 def lyapunov_rank(K
):
243 Compute the Lyapunov (or bilinearity) rank of this cone.
245 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
247 1. The dimension of the Lie algebra of the automorphism group of the
250 2. The dimension of the linear space of all Lyapunov-like
251 transformations on the cone.
255 A closed, convex polyhedral cone.
259 An integer representing the Lyapunov rank of the cone. If the
260 dimension of the ambient vector space is `n`, then the Lyapunov rank
261 will be between `1` and `n` inclusive; however a rank of `n-1` is
262 not possible (see the first reference).
266 In the references, the cones are always assumed to be proper. We
267 do not impose this restriction.
275 The codimension formula from the second reference is used. We find
276 all pairs `(x,s)` in the complementarity set of `K` such that `x`
277 and `s` are rays of our cone. It is known that these vectors are
278 sufficient to apply the codimension formula. Once we have all such
279 pairs, we "brute force" the codimension formula by finding all
280 linearly-independent `xs^{T}`.
284 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
285 cone and Lyapunov-like transformations, Mathematical Programming, 147
288 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
289 Improper Cone. Work in-progress.
291 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
292 optimality constraints for the cone of positive polynomials,
293 Mathematical Programming, Series B, 129 (2011) 5-31.
297 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
300 sage: positives = Cone([(1,)])
301 sage: lyapunov_rank(positives)
303 sage: quadrant = Cone([(1,0), (0,1)])
304 sage: lyapunov_rank(quadrant)
306 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
307 sage: lyapunov_rank(octant)
310 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
313 sage: R5 = VectorSpace(QQ, 5)
314 sage: gens = R5.basis() + [ -r for r in R5.basis() ]
316 sage: lyapunov_rank(K)
319 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
322 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
323 sage: lyapunov_rank(L31)
326 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
328 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
329 sage: lyapunov_rank(L3infty)
332 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
333 + 1` [Orlitzky/Gowda]_::
335 sage: K = Cone([(1,0,0,0,0)])
336 sage: lyapunov_rank(K)
338 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
341 A subspace (of dimension `m`) in `n` dimensions should have a
342 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
344 sage: e1 = (1,0,0,0,0)
345 sage: neg_e1 = (-1,0,0,0,0)
346 sage: e2 = (0,1,0,0,0)
347 sage: neg_e2 = (0,-1,0,0,0)
348 sage: zero = (0,0,0,0,0)
349 sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
350 sage: lyapunov_rank(K)
352 sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
355 The Lyapunov rank should be additive on a product of proper cones
358 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
359 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
360 sage: K = L31.cartesian_product(octant)
361 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
364 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
365 The cone ``K`` in the following example is isomorphic to the nonnegative
366 octant in `\mathbb{R}^{3}`::
368 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
369 sage: lyapunov_rank(K)
372 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
373 itself [Rudolf et al.]_::
375 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
376 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
381 The Lyapunov rank should be additive on a product of proper cones
384 sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
385 sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
386 sage: K = K1.cartesian_product(K2)
387 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
390 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
391 itself [Rudolf et al.]_::
393 sage: K = random_cone(max_dim=10, max_rays=10)
394 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
397 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
398 be any number between `1` and `n` inclusive, excluding `n-1`
399 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
400 trivial cone in a trivial space as well. However, in zero dimensions,
401 the Lyapunov rank of the trivial cone will be zero::
403 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
404 sage: b = lyapunov_rank(K)
405 sage: n = K.lattice_dim()
406 sage: (n == 0 or 1 <= b) and b <= n
411 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
412 Lyapunov rank `n-1` in `n` dimensions::
414 sage: K = random_cone(max_dim=10)
415 sage: b = lyapunov_rank(K)
416 sage: n = K.lattice_dim()
420 The calculation of the Lyapunov rank of an improper cone can be
421 reduced to that of a proper cone [Orlitzky/Gowda]_::
423 sage: K = random_cone(max_dim=10)
424 sage: actual = lyapunov_rank(K)
425 sage: K_S = project_span(K)
426 sage: P = project_span(K_S.dual()).dual()
427 sage: l = K.linear_subspace().dimension()
428 sage: codim = K.lattice_dim() - K.dim()
429 sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
430 sage: actual == expected
433 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
435 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
436 sage: lyapunov_rank(K) == len(LL(K))
444 l
= K
.linear_subspace().dimension()
447 # K is not solid, project onto its span.
451 beta
+= m
*(n
- m
) + (n
- m
)**2
454 # K is not pointed, project its dual onto its span.
455 K
= project_span(K
.dual()).dual()