]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
dbe369a8c2feabaee32bade898106641c3a3c3f0
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 pass # implement me lol
99 Compute the Lyapunov (or bilinearity) rank of this cone.
101 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
103 1. The dimension of the Lie algebra of the automorphism group of the
106 2. The dimension of the linear space of all Lyapunov-like
107 transformations on the cone.
111 A closed, convex polyhedral cone.
115 An integer representing the Lyapunov rank of the cone. If the
116 dimension of the ambient vector space is `n`, then the Lyapunov rank
117 will be between `1` and `n` inclusive; however a rank of `n-1` is
118 not possible (see the first reference).
122 In the references, the cones are always assumed to be proper. We
123 do not impose this restriction.
131 The codimension formula from the second reference is used. We find
132 all pairs `(x,s)` in the complementarity set of `K` such that `x`
133 and `s` are rays of our cone. It is known that these vectors are
134 sufficient to apply the codimension formula. Once we have all such
135 pairs, we "brute force" the codimension formula by finding all
136 linearly-independent `xs^{T}`.
140 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
141 cone and Lyapunov-like transformations, Mathematical Programming, 147
144 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
145 optimality constraints for the cone of positive polynomials,
146 Mathematical Programming, Series B, 129 (2011) 5-31.
150 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
153 sage: positives = Cone([(1,)])
154 sage: lyapunov_rank(positives)
156 sage: quadrant = Cone([(1,0), (0,1)])
157 sage: lyapunov_rank(quadrant)
159 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
160 sage: lyapunov_rank(octant)
163 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
166 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
167 sage: lyapunov_rank(L31)
170 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
172 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
173 sage: lyapunov_rank(L3infty)
176 The Lyapunov rank should be additive on a product of cones
179 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
180 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
181 sage: K = L31.cartesian_product(octant)
182 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
185 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
186 The cone ``K`` in the following example is isomorphic to the nonnegative
187 octant in `\mathbb{R}^{3}`::
189 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
190 sage: lyapunov_rank(K)
193 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
194 itself [Rudolf et al.]_::
196 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
197 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
202 The Lyapunov rank should be additive on a product of cones
205 sage: K1 = random_cone(max_dim=10, max_rays=10)
206 sage: K2 = random_cone(max_dim=10, max_rays=10)
207 sage: K = K1.cartesian_product(K2)
208 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
211 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
212 itself [Rudolf et al.]_::
214 sage: K = random_cone(max_dim=10, max_rays=10)
215 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
218 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
219 be any number between `1` and `n` inclusive, excluding `n-1`
220 [Gowda/Tao]_ (by accident, this holds for the trivial cone in a
221 trivial space as well)::
223 sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
224 sage: b = lyapunov_rank(K)
225 sage: n = K.lattice_dim()
226 sage: 1 <= b and b <= n
232 V
= K
.lattice().vector_space()
234 C_of_K
= discrete_complementarity_set(K
)
236 matrices
= [x
.tensor_product(s
) for (x
,s
) in C_of_K
]
238 # Sage doesn't think matrices are vectors, so we have to convert
239 # our matrices to vectors explicitly before we can figure out how
240 # many are linearly-indepenedent.
242 # The space W has the same base ring as V, but dimension
243 # dim(V)^2. So it has the same dimension as the space of linear
244 # transformations on V. In other words, it's just the right size
245 # to create an isomorphism between it and our matrices.
246 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
250 Convert a matrix to a vector isomorphically.
254 vectors
= [phi(m
) for m
in matrices
]
256 return (W
.dimension() - W
.span(vectors
).rank())