]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
eb5316330f6afefdc3b452482bfbd6d59deabbd2
3 def is_lyapunov_like(L
,K
):
5 Determine whether or not ``L`` is Lyapunov-like on ``K``.
7 We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
8 L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
9 `\left\langle x,s \right\rangle` in the complementarity set of
10 ``K``. It is known [Orlitzky]_ that this property need only be
11 checked for generators of ``K`` and its dual.
15 - ``L`` -- A linear transformation or matrix.
17 - ``K`` -- A polyhedral closed convex cone.
21 ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
22 and ``False`` otherwise.
26 If this function returns ``True``, then ``L`` is Lyapunov-like
27 on ``K``. However, if ``False`` is returned, that could mean one
28 of two things. The first is that ``L`` is definitely not
29 Lyapunov-like on ``K``. The second is more of an "I don't know"
30 answer, returned (for example) if we cannot prove that an inner
35 M. Orlitzky. The Lyapunov rank of an improper cone.
36 http://www.optimization-online.org/DB_HTML/2015/10/5135.html
40 The identity is always Lyapunov-like in a nontrivial space::
42 sage: set_random_seed()
43 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
44 sage: L = identity_matrix(K.lattice_dim())
45 sage: is_lyapunov_like(L,K)
48 As is the "zero" transformation::
50 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
51 sage: R = K.lattice().vector_space().base_ring()
52 sage: L = zero_matrix(R, K.lattice_dim())
53 sage: is_lyapunov_like(L,K)
56 Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
59 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
60 sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
64 return all([(L
*x
).inner_product(s
) == 0
65 for (x
,s
) in K
.discrete_complementarity_set()])
68 def motzkin_decomposition(K
):
70 Return the pair of components in the Motzkin decomposition of this cone.
72 Every convex cone is the direct sum of a strictly convex cone and a
73 linear subspace [Stoer-Witzgall]_. Return a pair ``(P,S)`` of cones
74 such that ``P`` is strictly convex, ``S`` is a subspace, and ``K``
75 is the direct sum of ``P`` and ``S``.
79 An ordered pair ``(P,S)`` of closed convex polyhedral cones where
80 ``P`` is strictly convex, ``S`` is a subspace, and ``K`` is the
81 direct sum of ``P`` and ``S``.
85 .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and
86 Optimization in Finite Dimensions I. Springer-Verlag, New
91 The nonnegative orthant is strictly convex, so it is its own
92 strictly convex component and its subspace component is trivial::
94 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
95 sage: (P,S) = motzkin_decomposition(K)
96 sage: K.is_equivalent(P)
101 Likewise, full spaces are their own subspace components::
103 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
104 sage: K.is_full_space()
106 sage: (P,S) = motzkin_decomposition(K)
107 sage: K.is_equivalent(S)
114 A random point in the cone should belong to either the strictly
115 convex component or the subspace component. If the point is nonzero,
116 it cannot be in both::
118 sage: set_random_seed()
119 sage: K = random_cone(max_ambient_dim=8)
120 sage: (P,S) = motzkin_decomposition(K)
121 sage: x = K.random_element()
122 sage: P.contains(x) or S.contains(x)
124 sage: x.is_zero() or (P.contains(x) != S.contains(x))
127 The strictly convex component should always be strictly convex, and
128 the subspace component should always be a subspace::
130 sage: set_random_seed()
131 sage: K = random_cone(max_ambient_dim=8)
132 sage: (P,S) = motzkin_decomposition(K)
133 sage: P.is_strictly_convex()
135 sage: S.lineality() == S.dim()
138 The generators of the components are obtained from orthogonal
139 projections of the original generators [Stoer-Witzgall]_::
141 sage: set_random_seed()
142 sage: K = random_cone(max_ambient_dim=8)
143 sage: (P,S) = motzkin_decomposition(K)
144 sage: A = S.linear_subspace().complement().matrix()
145 sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A
146 sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice())
147 sage: P.is_equivalent(expected_P)
149 sage: A = S.linear_subspace().matrix()
150 sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A
151 sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice())
152 sage: S.is_equivalent(expected_S)
155 # The lines() method only gives us one generator for each line,
156 # so we negate the result and combine everything for the full set.
157 S
= Cone([p
*l
for p
in [1,-1] for l
in K
.lines()], K
.lattice())
159 # Since ``S`` is a subspace, the rays of its dual generate its
160 # orthogonal complement.
161 P
= K
.intersection( Cone(S
.dual(), K
.lattice()) )
166 def positive_operator_gens(K
):
168 Compute generators of the cone of positive operators on this cone.
172 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
173 Each matrix ``P`` in the list should have the property that ``P*x``
174 is an element of ``K`` whenever ``x`` is an element of
175 ``K``. Moreover, any nonnegative linear combination of these
176 matrices shares the same property.
180 The trivial cone in a trivial space has no positive operators::
182 sage: K = Cone([], ToricLattice(0))
183 sage: positive_operator_gens(K)
186 Positive operators on the nonnegative orthant are nonnegative matrices::
188 sage: K = Cone([(1,)])
189 sage: positive_operator_gens(K)
192 sage: K = Cone([(1,0),(0,1)])
193 sage: positive_operator_gens(K)
195 [1 0] [0 1] [0 0] [0 0]
196 [0 0], [0 0], [1 0], [0 1]
199 Every operator is positive on the ambient vector space::
201 sage: K = Cone([(1,),(-1,)])
202 sage: K.is_full_space()
204 sage: positive_operator_gens(K)
207 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
208 sage: K.is_full_space()
210 sage: positive_operator_gens(K)
212 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
213 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
218 A positive operator on a cone should send its generators into the cone::
220 sage: set_random_seed()
221 sage: K = random_cone(max_ambient_dim=5)
222 sage: pi_of_K = positive_operator_gens(K)
223 sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
226 The dimension of the cone of positive operators is given by the
227 corollary in my paper::
229 sage: set_random_seed()
230 sage: K = random_cone(max_ambient_dim=5)
231 sage: n = K.lattice_dim()
233 sage: l = K.lineality()
234 sage: pi_of_K = positive_operator_gens(K)
235 sage: L = ToricLattice(n**2)
236 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
237 sage: expected = n**2 - l*(m - l) - (n - m)*m
238 sage: actual == expected
241 The lineality of the cone of positive operators is given by the
242 corollary in my paper::
244 sage: set_random_seed()
245 sage: K = random_cone(max_ambient_dim=5)
246 sage: n = K.lattice_dim()
247 sage: pi_of_K = positive_operator_gens(K)
248 sage: L = ToricLattice(n**2)
249 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
250 sage: expected = n**2 - K.dim()*K.dual().dim()
251 sage: actual == expected
254 The cone ``K`` is proper if and only if the cone of positive
255 operators on ``K`` is proper::
257 sage: set_random_seed()
258 sage: K = random_cone(max_ambient_dim=5)
259 sage: pi_of_K = positive_operator_gens(K)
260 sage: L = ToricLattice(K.lattice_dim()**2)
261 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
262 sage: K.is_proper() == pi_cone.is_proper()
265 # Matrices are not vectors in Sage, so we have to convert them
266 # to vectors explicitly before we can find a basis. We need these
267 # two values to construct the appropriate "long vector" space.
268 F
= K
.lattice().base_field()
271 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
273 # Convert those tensor products to long vectors.
274 W
= VectorSpace(F
, n
**2)
275 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
277 # Create the *dual* cone of the positive operators, expressed as
279 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()))
281 # Now compute the desired cone from its dual...
282 pi_cone
= pi_dual
.dual()
284 # And finally convert its rays back to matrix representations.
285 M
= MatrixSpace(F
, n
)
286 return [ M(v
.list()) for v
in pi_cone
.rays() ]
289 def Z_transformation_gens(K
):
291 Compute generators of the cone of Z-transformations on this cone.
295 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
296 Each matrix ``L`` in the list should have the property that
297 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
298 discrete complementarity set of ``K``. Moreover, any nonnegative
299 linear combination of these matrices shares the same property.
303 Z-transformations on the nonnegative orthant are just Z-matrices.
304 That is, matrices whose off-diagonal elements are nonnegative::
306 sage: K = Cone([(1,0),(0,1)])
307 sage: Z_transformation_gens(K)
309 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
310 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
312 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
313 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
314 ....: for i in range(z.nrows())
315 ....: for j in range(z.ncols())
319 The trivial cone in a trivial space has no Z-transformations::
321 sage: K = Cone([], ToricLattice(0))
322 sage: Z_transformation_gens(K)
325 Z-transformations on a subspace are Lyapunov-like and vice-versa::
327 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
328 sage: K.is_full_space()
330 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
331 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
337 The Z-property is possessed by every Z-transformation::
339 sage: set_random_seed()
340 sage: K = random_cone(max_ambient_dim=6)
341 sage: Z_of_K = Z_transformation_gens(K)
342 sage: dcs = K.discrete_complementarity_set()
343 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
344 ....: for (x,s) in dcs])
347 The lineality space of Z is LL::
349 sage: set_random_seed()
350 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
351 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
352 sage: z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ])
353 sage: z_cone.linear_subspace() == lls
356 And thus, the lineality of Z is the Lyapunov rank::
358 sage: set_random_seed()
359 sage: K = random_cone(max_ambient_dim=6)
360 sage: Z_of_K = Z_transformation_gens(K)
361 sage: L = ToricLattice(K.lattice_dim()**2)
362 sage: z_cone = Cone([ z.list() for z in Z_of_K ], lattice=L)
363 sage: z_cone.lineality() == K.lyapunov_rank()
366 The lineality spaces of pi-star and Z-star are equal:
368 sage: set_random_seed()
369 sage: K = random_cone(max_ambient_dim=5)
370 sage: pi_of_K = positive_operator_gens(K)
371 sage: Z_of_K = Z_transformation_gens(K)
372 sage: L = ToricLattice(K.lattice_dim()**2)
373 sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
374 sage: z_star = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
375 sage: pi_star.linear_subspace() == z_star.linear_subspace()
378 # Matrices are not vectors in Sage, so we have to convert them
379 # to vectors explicitly before we can find a basis. We need these
380 # two values to construct the appropriate "long vector" space.
381 F
= K
.lattice().base_field()
384 # These tensor products contain generators for the dual cone of
385 # the cross-positive transformations.
386 tensor_products
= [ s
.tensor_product(x
)
387 for (x
,s
) in K
.discrete_complementarity_set() ]
389 # Turn our matrices into long vectors...
390 W
= VectorSpace(F
, n
**2)
391 vectors
= [ W(m
.list()) for m
in tensor_products
]
393 # Create the *dual* cone of the cross-positive operators,
394 # expressed as long vectors..
395 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()))
397 # Now compute the desired cone from its dual...
398 Sigma_cone
= Sigma_dual
.dual()
400 # And finally convert its rays back to matrix representations.
401 # But first, make them negative, so we get Z-transformations and
402 # not cross-positive ones.
403 M
= MatrixSpace(F
, n
)
404 return [ -M(v
.list()) for v
in Sigma_cone
.rays() ]
408 gens
= Z_transformation_gens(K
)
412 return Cone([ g
.list() for g
in gens
], lattice
=L
)
415 gens
= positive_operator_gens(K
)
419 return Cone([ g
.list() for g
in gens
], lattice
=L
)