]>
gitweb.michael.orlitzky.com - sage.d.git/blob - cone/cone.py
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 returns one generator per line. For a true
156 # line, we also need a generator pointing in the opposite direction.
157 S_gens
= [ direction
*gen
for direction
in [1,-1] for gen
in K
.lines() ]
158 S
= Cone(S_gens
, K
.lattice())
160 # Since ``S`` is a subspace, the rays of its dual generate its
161 # orthogonal complement.
162 S_perp
= Cone(S
.dual(), K
.lattice())
163 P
= K
.intersection(S_perp
)
168 def positive_operator_gens(K
):
170 Compute generators of the cone of positive operators on this cone.
174 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
175 Each matrix ``P`` in the list should have the property that ``P*x``
176 is an element of ``K`` whenever ``x`` is an element of
177 ``K``. Moreover, any nonnegative linear combination of these
178 matrices shares the same property.
182 The trivial cone in a trivial space has no positive operators::
184 sage: K = Cone([], ToricLattice(0))
185 sage: positive_operator_gens(K)
188 Positive operators on the nonnegative orthant are nonnegative matrices::
190 sage: K = Cone([(1,)])
191 sage: positive_operator_gens(K)
194 sage: K = Cone([(1,0),(0,1)])
195 sage: positive_operator_gens(K)
197 [1 0] [0 1] [0 0] [0 0]
198 [0 0], [0 0], [1 0], [0 1]
201 Every operator is positive on the ambient vector space::
203 sage: K = Cone([(1,),(-1,)])
204 sage: K.is_full_space()
206 sage: positive_operator_gens(K)
209 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
210 sage: K.is_full_space()
212 sage: positive_operator_gens(K)
214 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
215 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
220 Each positive operator generator should send the generators of the
223 sage: set_random_seed()
224 sage: K = random_cone(max_ambient_dim=5)
225 sage: pi_of_K = positive_operator_gens(K)
226 sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
229 Each positive operator generator should send a random element of the
232 sage: set_random_seed()
233 sage: K = random_cone(max_ambient_dim=5)
234 sage: pi_of_K = positive_operator_gens(K)
235 sage: all([ K.contains(P*K.random_element()) for P in pi_of_K ])
238 A random element of the positive operator cone should send the
239 generators of the cone into the cone::
241 sage: set_random_seed()
242 sage: K = random_cone(max_ambient_dim=5)
243 sage: pi_of_K = positive_operator_gens(K)
244 sage: L = ToricLattice(K.lattice_dim()**2)
245 sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
246 sage: P = matrix(K.lattice_dim(), pi_cone.random_element().list())
247 sage: all([ K.contains(P*x) for x in K ])
250 A random element of the positive operator cone should send a random
251 element of the cone into the cone::
253 sage: set_random_seed()
254 sage: K = random_cone(max_ambient_dim=5)
255 sage: pi_of_K = positive_operator_gens(K)
256 sage: L = ToricLattice(K.lattice_dim()**2)
257 sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
258 sage: P = matrix(K.lattice_dim(), pi_cone.random_element().list())
259 sage: K.contains(P*K.random_element())
262 The dimension of the cone of positive operators is given by the
263 corollary in my paper::
265 sage: set_random_seed()
266 sage: K = random_cone(max_ambient_dim=5)
267 sage: n = K.lattice_dim()
269 sage: l = K.lineality()
270 sage: pi_of_K = positive_operator_gens(K)
271 sage: L = ToricLattice(n**2)
272 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
273 sage: expected = n**2 - l*(m - l) - (n - m)*m
274 sage: actual == expected
277 The lineality of the cone of positive operators is given by the
278 corollary in my paper::
280 sage: set_random_seed()
281 sage: K = random_cone(max_ambient_dim=5)
282 sage: n = K.lattice_dim()
283 sage: pi_of_K = positive_operator_gens(K)
284 sage: L = ToricLattice(n**2)
285 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
286 sage: expected = n**2 - K.dim()*K.dual().dim()
287 sage: actual == expected
290 The cone ``K`` is proper if and only if the cone of positive
291 operators on ``K`` is proper::
293 sage: set_random_seed()
294 sage: K = random_cone(max_ambient_dim=5)
295 sage: pi_of_K = positive_operator_gens(K)
296 sage: L = ToricLattice(K.lattice_dim()**2)
297 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
298 sage: K.is_proper() == pi_cone.is_proper()
301 # Matrices are not vectors in Sage, so we have to convert them
302 # to vectors explicitly before we can find a basis. We need these
303 # two values to construct the appropriate "long vector" space.
304 F
= K
.lattice().base_field()
307 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
309 # Convert those tensor products to long vectors.
310 W
= VectorSpace(F
, n
**2)
311 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
313 # Create the *dual* cone of the positive operators, expressed as
315 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()))
317 # Now compute the desired cone from its dual...
318 pi_cone
= pi_dual
.dual()
320 # And finally convert its rays back to matrix representations.
321 M
= MatrixSpace(F
, n
)
322 return [ M(v
.list()) for v
in pi_cone
.rays() ]
325 def Z_transformation_gens(K
):
327 Compute generators of the cone of Z-transformations on this cone.
331 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
332 Each matrix ``L`` in the list should have the property that
333 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
334 discrete complementarity set of ``K``. Moreover, any nonnegative
335 linear combination of these matrices shares the same property.
339 Z-transformations on the nonnegative orthant are just Z-matrices.
340 That is, matrices whose off-diagonal elements are nonnegative::
342 sage: K = Cone([(1,0),(0,1)])
343 sage: Z_transformation_gens(K)
345 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
346 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
348 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
349 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
350 ....: for i in range(z.nrows())
351 ....: for j in range(z.ncols())
355 The trivial cone in a trivial space has no Z-transformations::
357 sage: K = Cone([], ToricLattice(0))
358 sage: Z_transformation_gens(K)
361 Z-transformations on a subspace are Lyapunov-like and vice-versa::
363 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
364 sage: K.is_full_space()
366 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
367 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
373 The Z-property is possessed by every Z-transformation::
375 sage: set_random_seed()
376 sage: K = random_cone(max_ambient_dim=6)
377 sage: Z_of_K = Z_transformation_gens(K)
378 sage: dcs = K.discrete_complementarity_set()
379 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
380 ....: for (x,s) in dcs])
383 The lineality space of Z is LL::
385 sage: set_random_seed()
386 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
387 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
388 sage: z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ])
389 sage: z_cone.linear_subspace() == lls
392 And thus, the lineality of Z is the Lyapunov rank::
394 sage: set_random_seed()
395 sage: K = random_cone(max_ambient_dim=6)
396 sage: Z_of_K = Z_transformation_gens(K)
397 sage: L = ToricLattice(K.lattice_dim()**2)
398 sage: z_cone = Cone([ z.list() for z in Z_of_K ], lattice=L)
399 sage: z_cone.lineality() == K.lyapunov_rank()
402 The lineality spaces of pi-star and Z-star are equal:
404 sage: set_random_seed()
405 sage: K = random_cone(max_ambient_dim=5)
406 sage: pi_of_K = positive_operator_gens(K)
407 sage: Z_of_K = Z_transformation_gens(K)
408 sage: L = ToricLattice(K.lattice_dim()**2)
409 sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
410 sage: z_star = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
411 sage: pi_star.linear_subspace() == z_star.linear_subspace()
414 # Matrices are not vectors in Sage, so we have to convert them
415 # to vectors explicitly before we can find a basis. We need these
416 # two values to construct the appropriate "long vector" space.
417 F
= K
.lattice().base_field()
420 # These tensor products contain generators for the dual cone of
421 # the cross-positive transformations.
422 tensor_products
= [ s
.tensor_product(x
)
423 for (x
,s
) in K
.discrete_complementarity_set() ]
425 # Turn our matrices into long vectors...
426 W
= VectorSpace(F
, n
**2)
427 vectors
= [ W(m
.list()) for m
in tensor_products
]
429 # Create the *dual* cone of the cross-positive operators,
430 # expressed as long vectors..
431 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()))
433 # Now compute the desired cone from its dual...
434 Sigma_cone
= Sigma_dual
.dual()
436 # And finally convert its rays back to matrix representations.
437 # But first, make them negative, so we get Z-transformations and
438 # not cross-positive ones.
439 M
= MatrixSpace(F
, n
)
440 return [ -M(v
.list()) for v
in Sigma_cone
.rays() ]
444 gens
= Z_transformation_gens(K
)
448 return Cone([ g
.list() for g
in gens
], lattice
=L
)
451 gens
= positive_operator_gens(K
)
455 return Cone([ g
.list() for g
in gens
], lattice
=L
)