]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
507b6cea5c24d0a4eb745d9245c85a1fbf4ddb90
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
=None, max_dim
=None, min_rays
=None, 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. Any
17 parameters left unspecified will be chosen randomly.
21 - ``min_dim`` (default: random) -- The minimum dimension of the ambient
24 - ``max_dim`` (default: random) -- The maximum dimension of the ambient
27 - ``min_rays`` (default: random) -- The minimum number of generating rays
30 - ``max_rays`` (default: random) -- The maximum number of generating rays
35 A new, randomly generated cone.
39 It's hard to test the output of a random process, but we can at
40 least make sure that we get a cone back::
42 sage: from sage.geometry.cone import is_Cone
43 sage: K = random_cone()
44 sage: is_Cone(K) # long time
49 def random_min_max(l
,u
):
51 We need to handle four cases to prevent us from doing
52 something stupid like having an upper bound that's lower than
53 our lower bound. And we would need to repeat all of that logic
54 for the dimension/rays, so we consolidate it here.
56 if l
is None and u
is None:
57 # They're both random, just return a random nonnegative
59 return ZZ
.random_element().abs()
61 if l
is not None and u
is not None:
62 # Both were specified. Again, just make up a number and
63 # return it. If the user wants to give us u < l then he
64 # can have an exception.
65 return ZZ
.random_element(l
,u
)
67 if l
is not None and u
is None:
68 # In this case, we're generating the upper bound randomly
69 # GIVEN A LOWER BOUND. So we add a random nonnegative
70 # integer to the given lower bound.
71 u
= l
+ ZZ
.random_element().abs()
72 return ZZ
.random_element(l
,u
)
74 # Here we must be in the only remaining case, where we are
75 # given an upper bound but no lower bound. We might as well
77 return ZZ
.random_element(0,u
)
79 d
= random_min_max(min_dim
, max_dim
)
80 r
= random_min_max(min_rays
, max_rays
)
83 rays
= [L
.random_element() for i
in range(0,r
)]
85 # We pass the lattice in case there are no rays.
86 return Cone(rays
, lattice
=L
)
89 def discrete_complementarity_set(K
):
91 Compute the discrete complementarity set of this cone.
93 The complementarity set of this cone is the set of all orthogonal
94 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
95 dual. The discrete complementarity set restricts `x` and `s` to be
96 generators of their respective cones.
100 A list of pairs `(x,s)` such that,
102 * `x` is in this cone.
103 * `x` is a generator of this cone.
104 * `s` is in this cone's dual.
105 * `s` is a generator of this cone's dual.
106 * `x` and `s` are orthogonal.
110 The discrete complementarity set of the nonnegative orthant consists
111 of pairs of standard basis vectors::
113 sage: K = Cone([(1,0),(0,1)])
114 sage: discrete_complementarity_set(K)
115 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
117 If the cone consists of a single ray, the second components of the
118 discrete complementarity set should generate the orthogonal
119 complement of that ray::
121 sage: K = Cone([(1,0)])
122 sage: discrete_complementarity_set(K)
123 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
124 sage: K = Cone([(1,0,0)])
125 sage: discrete_complementarity_set(K)
126 [((1, 0, 0), (0, 1, 0)),
127 ((1, 0, 0), (0, -1, 0)),
128 ((1, 0, 0), (0, 0, 1)),
129 ((1, 0, 0), (0, 0, -1))]
131 When the cone is the entire space, its dual is the trivial cone, so
132 the discrete complementarity set is empty::
134 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
135 sage: discrete_complementarity_set(K)
140 The complementarity set of the dual can be obtained by switching the
141 components of the complementarity set of the original cone::
143 sage: K1 = random_cone(0,10,0,10)
145 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
146 sage: actual = discrete_complementarity_set(K1)
147 sage: actual == expected
151 V
= K
.lattice().vector_space()
153 # Convert the rays to vectors so that we can compute inner
155 xs
= [V(x
) for x
in K
.rays()]
156 ss
= [V(s
) for s
in K
.dual().rays()]
158 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
161 def lyapunov_rank(K
):
163 Compute the Lyapunov (or bilinearity) rank of this cone.
165 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
167 1. The dimension of the Lie algebra of the automorphism group of the
170 2. The dimension of the linear space of all Lyapunov-like
171 transformations on the cone.
175 A closed, convex polyhedral cone.
179 An integer representing the Lyapunov rank of the cone. If the
180 dimension of the ambient vector space is `n`, then the Lyapunov rank
181 will be between `1` and `n` inclusive; however a rank of `n-1` is
182 not possible (see the first reference).
186 In the references, the cones are always assumed to be proper. We
187 do not impose this restriction.
195 The codimension formula from the second reference is used. We find
196 all pairs `(x,s)` in the complementarity set of `K` such that `x`
197 and `s` are rays of our cone. It is known that these vectors are
198 sufficient to apply the codimension formula. Once we have all such
199 pairs, we "brute force" the codimension formula by finding all
200 linearly-independent `xs^{T}`.
204 1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
205 and Lyapunov-like transformations, Mathematical Programming, 147
208 2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
209 optimality constraints for the cone of positive polynomials,
210 Mathematical Programming, Series B, 129 (2011) 5-31.
214 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
216 sage: positives = Cone([(1,)])
217 sage: lyapunov_rank(positives)
219 sage: quadrant = Cone([(1,0), (0,1)])
220 sage: lyapunov_rank(quadrant)
222 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
223 sage: lyapunov_rank(octant)
226 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
228 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
229 sage: lyapunov_rank(L31)
232 Likewise for the `L^{3}_{\infty}` cone::
234 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
235 sage: lyapunov_rank(L3infty)
238 The Lyapunov rank should be additive on a product of cones::
240 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
241 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
242 sage: K = L31.cartesian_product(octant)
243 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
246 Two isomorphic cones should have the same Lyapunov rank. The cone
247 ``K`` in the following example is isomorphic to the nonnegative
248 octant in `\mathbb{R}^{3}`::
250 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
251 sage: lyapunov_rank(K)
254 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
257 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
258 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
263 The Lyapunov rank should be additive on a product of cones::
265 sage: K1 = random_cone(0,10,0,10)
266 sage: K2 = random_cone(0,10,0,10)
267 sage: K = K1.cartesian_product(K2)
268 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
271 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
274 sage: K = random_cone(0,10,0,10)
275 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
279 V
= K
.lattice().vector_space()
281 C_of_K
= discrete_complementarity_set(K
)
283 matrices
= [x
.tensor_product(s
) for (x
,s
) in C_of_K
]
285 # Sage doesn't think matrices are vectors, so we have to convert
286 # our matrices to vectors explicitly before we can figure out how
287 # many are linearly-indepenedent.
289 # The space W has the same base ring as V, but dimension
290 # dim(V)^2. So it has the same dimension as the space of linear
291 # transformations on V. In other words, it's just the right size
292 # to create an isomorphism between it and our matrices.
293 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
297 Convert a matrix to a vector isomorphically.
301 vectors
= [phi(m
) for m
in matrices
]
303 return (W
.dimension() - W
.span(vectors
).rank())