]>
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('../../'))
11 def random_cone(min_dim
=0, max_dim
=None, min_rays
=0, max_rays
=None):
13 Generate a random rational convex polyhedral cone.
15 Lower and upper bounds may be provided for both the dimension of the
16 ambient space and the number of generating rays of the cone. If a
17 lower bound is left unspecified, it defaults to zero. Unspecified
18 upper bounds will be chosen randomly.
20 The lower bound on the number of rays is limited to twice the
21 maximum dimension of the ambient vector space. To see why, consider
22 the space $\mathbb{R}^{2}$, and suppose we have generated four rays,
23 $\left\{ \pm e_{1}, \pm e_{2} \right\}$. Clearly any other ray in
24 the space is a nonnegative linear combination of those four,
25 so it is hopeless to generate more. It is therefore an error
26 to request more in the form of ``min_rays``.
30 If you do not explicitly request more than ``2 * max_dim`` rays,
31 a larger number may still be randomly generated. In that case,
32 the returned cone will simply be equal to the entire space.
36 - ``min_dim`` (default: zero) -- A nonnegative integer representing the
37 minimum dimension of the ambient lattice.
39 - ``max_dim`` (default: random) -- A nonnegative integer representing
40 the maximum dimension of the ambient
43 - ``min_rays`` (default: zero) -- A nonnegative integer representing the
44 minimum number of generating rays of the
47 - ``max_rays`` (default: random) -- A nonnegative integer representing the
48 maximum number of generating rays of
53 A new, randomly generated cone.
55 A ``ValueError` will be thrown under the following conditions:
57 * Any of ``min_dim``, ``max_dim``, ``min_rays``, or ``max_rays``
60 * ``max_dim`` is less than ``min_dim``.
62 * ``max_rays`` is less than ``min_rays``.
64 * ``min_rays`` is greater than twice ``max_dim``.
68 If we set the lower/upper bounds to zero, then our result is
71 sage: random_cone(0,0,0,0)
72 0-d cone in 0-d lattice N
74 We can predict the dimension when ``min_dim == max_dim``::
76 sage: random_cone(min_dim=4, max_dim=4, min_rays=0, max_rays=0)
77 0-d cone in 4-d lattice N
79 Likewise for the number of rays when ``min_rays == max_rays``::
81 sage: random_cone(min_dim=10, max_dim=10, min_rays=10, max_rays=10)
82 10-d cone in 10-d lattice N
86 It's hard to test the output of a random process, but we can at
87 least make sure that we get a cone back::
89 sage: from sage.geometry.cone import is_Cone # long time
90 sage: K = random_cone() # long time
91 sage: is_Cone(K) # long time
94 The upper/lower bounds are respected::
96 sage: K = random_cone(min_dim=5, max_dim=10, min_rays=3, max_rays=4)
97 sage: 5 <= K.lattice_dim() and K.lattice_dim() <= 10
99 sage: 3 <= K.nrays() and K.nrays() <= 4
102 Ensure that an exception is raised when either lower bound is greater
103 than its respective upper bound::
105 sage: random_cone(min_dim=5, max_dim=2)
106 Traceback (most recent call last):
108 ValueError: max_dim cannot be less than min_dim.
110 sage: random_cone(min_rays=5, max_rays=2)
111 Traceback (most recent call last):
113 ValueError: max_rays cannot be less than min_rays.
115 And if we request too many rays::
117 sage: random_cone(min_rays=5, max_dim=1)
118 Traceback (most recent call last):
120 ValueError: min_rays cannot be larger than twice max_dim.
124 # Catch obvious mistakes so that we can generate clear error
128 raise ValueError('min_dim must be nonnegative.')
131 raise ValueError('min_rays must be nonnegative.')
133 if max_dim
is not None:
135 raise ValueError('max_dim must be nonnegative.')
136 if (max_dim
< min_dim
):
137 raise ValueError('max_dim cannot be less than min_dim.')
138 if min_rays
> 2*max_dim
:
139 raise ValueError('min_rays cannot be larger than twice max_dim.')
141 if max_rays
is not None:
143 raise ValueError('max_rays must be nonnegative.')
144 if (max_rays
< min_rays
):
145 raise ValueError('max_rays cannot be less than min_rays.')
148 def random_min_max(l
,u
):
150 We need to handle two cases for the upper bounds, and we need to do
151 the same thing for max_dim/max_rays. So we consolidate the logic here.
154 # The upper bound is unspecified; return a random integer
156 return l
+ ZZ
.random_element().abs()
158 # We have an upper bound, and it's greater than or equal
159 # to our lower bound. So we generate a random integer in
160 # [0,u-l], and then add it to l to get something in
161 # [l,u]. To understand the "+1", check the
162 # ZZ.random_element() docs.
163 return l
+ ZZ
.random_element(u
- l
+ 1)
165 def is_full_space(K
):
167 Is this cone equivalent to the full ambient vector space?
169 return K
.lines().dim() == K
.lattice_dim()
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
181 rays
= [L
.random_element() for i
in range(0, r
)]
183 # (The lattice parameter is required when no rays are given, so we
184 # pass it just in case ``r == 0``).
185 K
= Cone(rays
, lattice
=L
)
187 # Now if we generated two of the "same" rays, we'll have fewer
188 # generating rays than ``r``. In that case, we keep making up new
189 # rays and recreating the cone until we get the right number of
190 # independent generators. We can obviously stop if ``K`` is the
191 # entire ambient vector space.
192 while r
> K
.nrays() and not is_full_space(K
):
193 rays
.append(L
.random_element())
199 def discrete_complementarity_set(K
):
201 Compute the discrete complementarity set of this cone.
203 The complementarity set of this cone is the set of all orthogonal
204 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
205 dual. The discrete complementarity set restricts `x` and `s` to be
206 generators of their respective cones.
210 A list of pairs `(x,s)` such that,
212 * `x` is in this cone.
213 * `x` is a generator of this cone.
214 * `s` is in this cone's dual.
215 * `s` is a generator of this cone's dual.
216 * `x` and `s` are orthogonal.
220 The discrete complementarity set of the nonnegative orthant consists
221 of pairs of standard basis vectors::
223 sage: K = Cone([(1,0),(0,1)])
224 sage: discrete_complementarity_set(K)
225 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
227 If the cone consists of a single ray, the second components of the
228 discrete complementarity set should generate the orthogonal
229 complement of that ray::
231 sage: K = Cone([(1,0)])
232 sage: discrete_complementarity_set(K)
233 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
234 sage: K = Cone([(1,0,0)])
235 sage: discrete_complementarity_set(K)
236 [((1, 0, 0), (0, 1, 0)),
237 ((1, 0, 0), (0, -1, 0)),
238 ((1, 0, 0), (0, 0, 1)),
239 ((1, 0, 0), (0, 0, -1))]
241 When the cone is the entire space, its dual is the trivial cone, so
242 the discrete complementarity set is empty::
244 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
245 sage: discrete_complementarity_set(K)
250 The complementarity set of the dual can be obtained by switching the
251 components of the complementarity set of the original cone::
253 sage: K1 = random_cone(max_dim=10, max_rays=10)
255 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
256 sage: actual = discrete_complementarity_set(K1)
257 sage: actual == expected
261 V
= K
.lattice().vector_space()
263 # Convert the rays to vectors so that we can compute inner
265 xs
= [V(x
) for x
in K
.rays()]
266 ss
= [V(s
) for s
in K
.dual().rays()]
268 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
271 def lyapunov_rank(K
):
273 Compute the Lyapunov (or bilinearity) rank of this cone.
275 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
277 1. The dimension of the Lie algebra of the automorphism group of the
280 2. The dimension of the linear space of all Lyapunov-like
281 transformations on the cone.
285 A closed, convex polyhedral cone.
289 An integer representing the Lyapunov rank of the cone. If the
290 dimension of the ambient vector space is `n`, then the Lyapunov rank
291 will be between `1` and `n` inclusive; however a rank of `n-1` is
292 not possible (see the first reference).
296 In the references, the cones are always assumed to be proper. We
297 do not impose this restriction.
305 The codimension formula from the second reference is used. We find
306 all pairs `(x,s)` in the complementarity set of `K` such that `x`
307 and `s` are rays of our cone. It is known that these vectors are
308 sufficient to apply the codimension formula. Once we have all such
309 pairs, we "brute force" the codimension formula by finding all
310 linearly-independent `xs^{T}`.
314 1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
315 and Lyapunov-like transformations, Mathematical Programming, 147
318 2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
319 optimality constraints for the cone of positive polynomials,
320 Mathematical Programming, Series B, 129 (2011) 5-31.
324 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
326 sage: positives = Cone([(1,)])
327 sage: lyapunov_rank(positives)
329 sage: quadrant = Cone([(1,0), (0,1)])
330 sage: lyapunov_rank(quadrant)
332 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
333 sage: lyapunov_rank(octant)
336 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
338 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
339 sage: lyapunov_rank(L31)
342 Likewise for the `L^{3}_{\infty}` cone::
344 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
345 sage: lyapunov_rank(L3infty)
348 The Lyapunov rank should be additive on a product of cones::
350 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
351 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
352 sage: K = L31.cartesian_product(octant)
353 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
356 Two isomorphic cones should have the same Lyapunov rank. The cone
357 ``K`` in the following example is isomorphic to the nonnegative
358 octant in `\mathbb{R}^{3}`::
360 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
361 sage: lyapunov_rank(K)
364 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
367 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
368 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
373 The Lyapunov rank should be additive on a product of cones::
375 sage: K1 = random_cone(max_dim=10, max_rays=10)
376 sage: K2 = random_cone(max_dim=10, max_rays=10)
377 sage: K = K1.cartesian_product(K2)
378 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
381 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
384 sage: K = random_cone(max_dim=10, max_rays=10)
385 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
389 V
= K
.lattice().vector_space()
391 C_of_K
= discrete_complementarity_set(K
)
393 matrices
= [x
.tensor_product(s
) for (x
,s
) in C_of_K
]
395 # Sage doesn't think matrices are vectors, so we have to convert
396 # our matrices to vectors explicitly before we can figure out how
397 # many are linearly-indepenedent.
399 # The space W has the same base ring as V, but dimension
400 # dim(V)^2. So it has the same dimension as the space of linear
401 # transformations on V. In other words, it's just the right size
402 # to create an isomorphism between it and our matrices.
403 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
407 Convert a matrix to a vector isomorphically.
411 vectors
= [phi(m
) for m
in matrices
]
413 return (W
.dimension() - W
.span(vectors
).rank())