]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
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 # TODO: This test fails, maybe due to a bug in the existing cone code.
11 # If we request enough generators to span the space, then the returned
12 # cone should equal the ambient space::
14 # sage: K = random_cone(min_dim=5, max_dim=5, min_rays=10, max_rays=10)
15 # sage: K.lines().dimension() == K.lattice_dim()
18 def random_cone(min_dim
=0, max_dim
=None, min_rays
=0, max_rays
=None):
20 Generate a random rational convex polyhedral cone.
22 Lower and upper bounds may be provided for both the dimension of the
23 ambient space and the number of generating rays of the cone. If a
24 lower bound is left unspecified, it defaults to zero. Unspecified
25 upper bounds will be chosen randomly.
27 The number of generating rays is naturally limited to twice the
28 dimension of the ambient space. Take for example $\mathbb{R}^{2}$.
29 You could have the generators $\left\{ \pm e_{1}, \pm e_{2}
30 \right\}$, with cardinality $4 = 2 \cdot 2$; however any other ray
31 in the space is a nonnegative linear combination of those four.
35 If you do not explicitly request more than ``2 * max_dim`` rays,
36 a larger number may still be randomly generated. In that case,
37 the returned cone will simply be equal to the entire space.
41 - ``min_dim`` (default: zero) -- A nonnegative integer representing the
42 minimum dimension of the ambient lattice.
44 - ``max_dim`` (default: random) -- A nonnegative integer representing
45 the maximum dimension of the ambient
48 - ``min_rays`` (default: zero) -- A nonnegative integer representing the
49 minimum number of generating rays of the
52 - ``max_rays`` (default: random) -- A nonnegative integer representing the
53 maximum number of generating rays of
58 A new, randomly generated cone.
60 A ``ValueError` will be thrown under the following conditions:
62 * Any of ``min_dim``, ``max_dim``, ``min_rays``, or ``max_rays``
65 * ``max_dim`` is less than ``min_dim``.
67 * ``max_rays`` is less than ``min_rays``.
69 * ``min_rays`` is greater than twice ``max_dim``.
73 If we set the lower/upper bounds to zero, then our result is
76 sage: random_cone(0,0,0,0)
77 0-d cone in 0-d lattice N
79 We can predict the dimension when ``min_dim == max_dim``::
81 sage: random_cone(min_dim=4, max_dim=4, min_rays=0, max_rays=0)
82 0-d cone in 4-d lattice N
84 Likewise for the number of rays when ``min_rays == max_rays``::
86 sage: random_cone(min_dim=10, max_dim=10, min_rays=10, max_rays=10)
87 10-d cone in 10-d lattice N
91 It's hard to test the output of a random process, but we can at
92 least make sure that we get a cone back::
94 sage: from sage.geometry.cone import is_Cone # long time
95 sage: K = random_cone() # long time
96 sage: is_Cone(K) # long time
99 The upper/lower bounds are respected::
101 sage: K = random_cone(min_dim=5, max_dim=10, min_rays=3, max_rays=4)
102 sage: 5 <= K.lattice_dim() and K.lattice_dim() <= 10
104 sage: 3 <= K.nrays() and K.nrays() <= 4
107 Ensure that an exception is raised when either lower bound is greater
108 than its respective upper bound::
110 sage: random_cone(min_dim=5, max_dim=2)
111 Traceback (most recent call last):
113 ValueError: max_dim cannot be less than min_dim.
115 sage: random_cone(min_rays=5, max_rays=2)
116 Traceback (most recent call last):
118 ValueError: max_rays cannot be less than min_rays.
120 And if we request too many rays::
122 sage: random_cone(min_rays=5, max_dim=1)
123 Traceback (most recent call last):
125 ValueError: min_rays cannot be larger than twice max_dim.
129 # Catch obvious mistakes so that we can generate clear error
133 raise ValueError('min_dim must be nonnegative.')
136 raise ValueError('min_rays must be nonnegative.')
138 if max_dim
is not None:
140 raise ValueError('max_dim must be nonnegative.')
141 if (max_dim
< min_dim
):
142 raise ValueError('max_dim cannot be less than min_dim.')
143 if min_rays
> 2*max_dim
:
144 raise ValueError('min_rays cannot be larger than twice max_dim.')
146 if max_rays
is not None:
148 raise ValueError('max_rays must be nonnegative.')
149 if (max_rays
< min_rays
):
150 raise ValueError('max_rays cannot be less than min_rays.')
153 def random_min_max(l
,u
):
155 We need to handle two cases for the upper bounds, and we need to do
156 the same thing for max_dim/max_rays. So we consolidate the logic here.
159 # The upper bound is unspecified; return a random integer
161 return l
+ ZZ
.random_element().abs()
163 # We have an upper bound, and it's greater than or equal
164 # to our lower bound. So we generate a random integer in
165 # [0,u-l], and then add it to l to get something in
166 # [l,u]. To understand the "+1", check the
167 # ZZ.random_element() docs.
168 return l
+ ZZ
.random_element(u
- l
+ 1)
171 d
= random_min_max(min_dim
, max_dim
)
172 r
= random_min_max(min_rays
, max_rays
)
176 # The rays are trickier to generate, since we could generate v and
177 # 2*v as our "two rays." In that case, the resuting cone would
178 # have one generating ray. To avoid such a situation, we start by
179 # generating ``r`` rays where ``r`` is the number we want to end
182 # However, since we're going to *check* whether or not we actually
183 # have ``r``, we need ``r`` rays to be attainable. So we need to
184 # limit ``r`` to twice the dimension of the ambient space.
187 rays
= [L
.random_element() for i
in range(0, r
)]
189 # (The lattice parameter is required when no rays are given, so we
190 # pass it just in case ``r == 0``).
191 K
= Cone(rays
, lattice
=L
)
193 # Now if we generated two of the "same" rays, we'll have fewer
194 # generating rays than ``r``. In that case, we keep making up new
195 # rays and recreating the cone until we get the right number of
196 # independent generators.
198 rays
.append(L
.random_element())
204 def discrete_complementarity_set(K
):
206 Compute the discrete complementarity set of this cone.
208 The complementarity set of this cone is the set of all orthogonal
209 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
210 dual. The discrete complementarity set restricts `x` and `s` to be
211 generators of their respective cones.
215 A list of pairs `(x,s)` such that,
217 * `x` is in this cone.
218 * `x` is a generator of this cone.
219 * `s` is in this cone's dual.
220 * `s` is a generator of this cone's dual.
221 * `x` and `s` are orthogonal.
225 The discrete complementarity set of the nonnegative orthant consists
226 of pairs of standard basis vectors::
228 sage: K = Cone([(1,0),(0,1)])
229 sage: discrete_complementarity_set(K)
230 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
232 If the cone consists of a single ray, the second components of the
233 discrete complementarity set should generate the orthogonal
234 complement of that ray::
236 sage: K = Cone([(1,0)])
237 sage: discrete_complementarity_set(K)
238 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
239 sage: K = Cone([(1,0,0)])
240 sage: discrete_complementarity_set(K)
241 [((1, 0, 0), (0, 1, 0)),
242 ((1, 0, 0), (0, -1, 0)),
243 ((1, 0, 0), (0, 0, 1)),
244 ((1, 0, 0), (0, 0, -1))]
246 When the cone is the entire space, its dual is the trivial cone, so
247 the discrete complementarity set is empty::
249 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
250 sage: discrete_complementarity_set(K)
255 The complementarity set of the dual can be obtained by switching the
256 components of the complementarity set of the original cone::
258 sage: K1 = random_cone(max_dim=10, max_rays=10)
260 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
261 sage: actual = discrete_complementarity_set(K1)
262 sage: actual == expected
266 V
= K
.lattice().vector_space()
268 # Convert the rays to vectors so that we can compute inner
270 xs
= [V(x
) for x
in K
.rays()]
271 ss
= [V(s
) for s
in K
.dual().rays()]
273 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
276 def lyapunov_rank(K
):
278 Compute the Lyapunov (or bilinearity) rank of this cone.
280 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
282 1. The dimension of the Lie algebra of the automorphism group of the
285 2. The dimension of the linear space of all Lyapunov-like
286 transformations on the cone.
290 A closed, convex polyhedral cone.
294 An integer representing the Lyapunov rank of the cone. If the
295 dimension of the ambient vector space is `n`, then the Lyapunov rank
296 will be between `1` and `n` inclusive; however a rank of `n-1` is
297 not possible (see the first reference).
301 In the references, the cones are always assumed to be proper. We
302 do not impose this restriction.
310 The codimension formula from the second reference is used. We find
311 all pairs `(x,s)` in the complementarity set of `K` such that `x`
312 and `s` are rays of our cone. It is known that these vectors are
313 sufficient to apply the codimension formula. Once we have all such
314 pairs, we "brute force" the codimension formula by finding all
315 linearly-independent `xs^{T}`.
319 1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
320 and Lyapunov-like transformations, Mathematical Programming, 147
323 2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
324 optimality constraints for the cone of positive polynomials,
325 Mathematical Programming, Series B, 129 (2011) 5-31.
329 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
331 sage: positives = Cone([(1,)])
332 sage: lyapunov_rank(positives)
334 sage: quadrant = Cone([(1,0), (0,1)])
335 sage: lyapunov_rank(quadrant)
337 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
338 sage: lyapunov_rank(octant)
341 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
343 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
344 sage: lyapunov_rank(L31)
347 Likewise for the `L^{3}_{\infty}` cone::
349 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
350 sage: lyapunov_rank(L3infty)
353 The Lyapunov rank should be additive on a product of cones::
355 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
356 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
357 sage: K = L31.cartesian_product(octant)
358 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
361 Two isomorphic cones should have the same Lyapunov rank. The cone
362 ``K`` in the following example is isomorphic to the nonnegative
363 octant in `\mathbb{R}^{3}`::
365 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
366 sage: lyapunov_rank(K)
369 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
372 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
373 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
378 The Lyapunov rank should be additive on a product of cones::
380 sage: K1 = random_cone(max_dim=10, max_rays=10)
381 sage: K2 = random_cone(max_dim=10, max_rays=10)
382 sage: K = K1.cartesian_product(K2)
383 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
386 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
389 sage: K = random_cone(max_dim=10, max_rays=10)
390 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
394 V
= K
.lattice().vector_space()
396 C_of_K
= discrete_complementarity_set(K
)
398 matrices
= [x
.tensor_product(s
) for (x
,s
) in C_of_K
]
400 # Sage doesn't think matrices are vectors, so we have to convert
401 # our matrices to vectors explicitly before we can figure out how
402 # many are linearly-indepenedent.
404 # The space W has the same base ring as V, but dimension
405 # dim(V)^2. So it has the same dimension as the space of linear
406 # transformations on V. In other words, it's just the right size
407 # to create an isomorphism between it and our matrices.
408 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
412 Convert a matrix to a vector isomorphically.
416 vectors
= [phi(m
) for m
in matrices
]
418 return (W
.dimension() - W
.span(vectors
).rank())