]>
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 discrete_complementarity_set(K
):
13 Compute the discrete complementarity set of this cone.
15 The complementarity set of this cone is the set of all orthogonal
16 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
17 dual. The discrete complementarity set restricts `x` and `s` to be
18 generators of their respective cones.
22 A list of pairs `(x,s)` such that,
24 * `x` is in this cone.
25 * `x` is a generator of this cone.
26 * `s` is in this cone's dual.
27 * `s` is a generator of this cone's dual.
28 * `x` and `s` are orthogonal.
32 The discrete complementarity set of the nonnegative orthant consists
33 of pairs of standard basis vectors::
35 sage: K = Cone([(1,0),(0,1)])
36 sage: discrete_complementarity_set(K)
37 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
39 If the cone consists of a single ray, the second components of the
40 discrete complementarity set should generate the orthogonal
41 complement of that ray::
43 sage: K = Cone([(1,0)])
44 sage: discrete_complementarity_set(K)
45 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
46 sage: K = Cone([(1,0,0)])
47 sage: discrete_complementarity_set(K)
48 [((1, 0, 0), (0, 1, 0)),
49 ((1, 0, 0), (0, -1, 0)),
50 ((1, 0, 0), (0, 0, 1)),
51 ((1, 0, 0), (0, 0, -1))]
53 When the cone is the entire space, its dual is the trivial cone, so
54 the discrete complementarity set is empty::
56 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
57 sage: discrete_complementarity_set(K)
62 The complementarity set of the dual can be obtained by switching the
63 components of the complementarity set of the original cone::
65 sage: K1 = random_cone(max_dim=10, max_rays=10)
67 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
68 sage: actual = discrete_complementarity_set(K1)
69 sage: actual == expected
73 V
= K
.lattice().vector_space()
75 # Convert the rays to vectors so that we can compute inner
77 xs
= [V(x
) for x
in K
.rays()]
78 ss
= [V(s
) for s
in K
.dual().rays()]
80 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
85 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
90 A ``MatrixSpace`` object `M` such that every matrix `L \in M` is
91 Lyapunov-like on this cone.
94 V
= K
.lattice().vector_space()
96 C_of_K
= discrete_complementarity_set(K
)
98 tensor_products
= [s
.tensor_product(x
) for (x
,s
) in C_of_K
]
100 # Sage doesn't think matrices are vectors, so we have to convert
101 # our matrices to vectors explicitly before we can figure out how
102 # many are linearly-indepenedent.
104 # The space W has the same base ring as V, but dimension
105 # dim(V)^2. So it has the same dimension as the space of linear
106 # transformations on V. In other words, it's just the right size
107 # to create an isomorphism between it and our matrices.
108 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
110 # Turn our matrices into long vectors...
111 vectors
= [ W(m
.list()) for m
in tensor_products
]
113 # Vector space representation of Lyapunov-like matrices
114 # (i.e. vec(L) where L is Luapunov-like).
115 LL_vector
= W
.span(vectors
).complement()
117 # Now construct an ambient MatrixSpace in which to stick our
119 M
= MatrixSpace(V
.base_ring(), V
.dimension())
121 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
127 def lyapunov_rank(K
):
129 Compute the Lyapunov (or bilinearity) rank of this cone.
131 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
133 1. The dimension of the Lie algebra of the automorphism group of the
136 2. The dimension of the linear space of all Lyapunov-like
137 transformations on the cone.
141 A closed, convex polyhedral cone.
145 An integer representing the Lyapunov rank of the cone. If the
146 dimension of the ambient vector space is `n`, then the Lyapunov rank
147 will be between `1` and `n` inclusive; however a rank of `n-1` is
148 not possible (see the first reference).
152 In the references, the cones are always assumed to be proper. We
153 do not impose this restriction.
161 The codimension formula from the second reference is used. We find
162 all pairs `(x,s)` in the complementarity set of `K` such that `x`
163 and `s` are rays of our cone. It is known that these vectors are
164 sufficient to apply the codimension formula. Once we have all such
165 pairs, we "brute force" the codimension formula by finding all
166 linearly-independent `xs^{T}`.
170 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
171 cone and Lyapunov-like transformations, Mathematical Programming, 147
174 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
175 optimality constraints for the cone of positive polynomials,
176 Mathematical Programming, Series B, 129 (2011) 5-31.
180 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
183 sage: positives = Cone([(1,)])
184 sage: lyapunov_rank(positives)
186 sage: quadrant = Cone([(1,0), (0,1)])
187 sage: lyapunov_rank(quadrant)
189 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
190 sage: lyapunov_rank(octant)
193 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
196 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
197 sage: lyapunov_rank(L31)
200 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
202 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
203 sage: lyapunov_rank(L3infty)
206 The Lyapunov rank should be additive on a product of cones
209 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
210 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
211 sage: K = L31.cartesian_product(octant)
212 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
215 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
216 The cone ``K`` in the following example is isomorphic to the nonnegative
217 octant in `\mathbb{R}^{3}`::
219 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
220 sage: lyapunov_rank(K)
223 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
224 itself [Rudolf et al.]_::
226 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
227 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
232 The Lyapunov rank should be additive on a product of cones
235 sage: K1 = random_cone(max_dim=10, max_rays=10)
236 sage: K2 = random_cone(max_dim=10, max_rays=10)
237 sage: K = K1.cartesian_product(K2)
238 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
241 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
242 itself [Rudolf et al.]_::
244 sage: K = random_cone(max_dim=10, max_rays=10)
245 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
248 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
249 be any number between `1` and `n` inclusive, excluding `n-1`
250 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
251 trivial cone in a trivial space as well. However, in zero dimensions,
252 the Lyapunov rank of the trivial cone will be zero::
254 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
255 sage: b = lyapunov_rank(K)
256 sage: n = K.lattice_dim()
257 sage: (n == 0 or 1 <= b) and b <= n