]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/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 random_element(K
):
70 Return a random element of ``K`` from its ambient vector space.
74 The cone ``K`` is specified in terms of its generators, so that
75 ``K`` is equal to the convex conic combination of those generators.
76 To choose a random element of ``K``, we assign random nonnegative
77 coefficients to each generator of ``K`` and construct a new vector
80 A vector, rather than a ray, is returned so that the element may
81 have non-integer coordinates. Thus the element may have an
82 arbitrarily small norm.
86 A random element of the trivial cone is zero::
88 sage: set_random_seed()
89 sage: K = Cone([], ToricLattice(0))
90 sage: random_element(K)
92 sage: K = Cone([(0,)])
93 sage: random_element(K)
95 sage: K = Cone([(0,0)])
96 sage: random_element(K)
98 sage: K = Cone([(0,0,0)])
99 sage: random_element(K)
102 A random element of the nonnegative orthant should have all
103 components nonnegative::
105 sage: set_random_seed()
106 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
107 sage: all([ x >= 0 for x in random_element(K) ])
112 Any cone should contain a random element of itself::
114 sage: set_random_seed()
115 sage: K = random_cone(max_ambient_dim=8)
116 sage: K.contains(random_element(K))
119 A strictly convex cone contains no lines, and thus no negative
120 multiples of any of its elements besides zero::
122 sage: set_random_seed()
123 sage: K = random_cone(max_ambient_dim=8, strictly_convex=True)
124 sage: x = random_element(K)
125 sage: x.is_zero() or not K.contains(-x)
128 The sum of random elements of a cone lies in the cone::
130 sage: set_random_seed()
131 sage: K = random_cone(max_ambient_dim=8)
132 sage: K.contains(sum([random_element(K) for i in range(10)]))
136 V
= K
.lattice().vector_space()
137 scaled_gens
= [ V
.base_field().random_element().abs()*V(r
) for r
in K
]
139 # Make sure we return a vector. Without the coercion, we might
140 # return ``0`` when ``K`` has no rays.
141 return V(sum(scaled_gens
))
144 def motzkin_decomposition(K
):
146 Every convex cone is the direct sum of a pointed cone and a linear
147 subspace. Return a pair ``(P,S)`` of cones such that ``P`` is
148 pointed, ``S`` is a subspace, and ``K`` is the direct sum of ``P``
153 An ordered pair ``(P,S)`` of closed convex polyhedral cones where
154 ``P`` is pointed, ``S`` is a subspace, and ``K`` is the direct sum
159 A random point in the cone should belong to either the pointed
160 subcone ``P`` or the subspace ``S``. If the point is nonzero, it
161 should lie in one but not both of them::
163 sage: set_random_seed()
164 sage: K = random_cone(max_ambient_dim=8)
165 sage: (P,S) = motzkin_decomposition(K)
166 sage: x = random_element(K)
167 sage: P.contains(x) or S.contains(x)
169 sage: x.is_zero() or (P.contains(x) != S.contains(x))
172 linspace_gens
= [ copy(b
) for b
in K
.linear_subspace().basis() ]
173 linspace_gens
+= [ -b
for b
in linspace_gens
]
175 S
= Cone(linspace_gens
, K
.lattice())
177 # Since ``S`` is a subspace, its dual is its orthogonal complement
178 # (albeit in the wrong lattice).
179 S_perp
= Cone(S
.dual(), K
.lattice())
180 P
= K
.intersection(S_perp
)
184 def positive_operator_gens(K
):
186 Compute generators of the cone of positive operators on this cone.
190 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
191 Each matrix ``P`` in the list should have the property that ``P*x``
192 is an element of ``K`` whenever ``x`` is an element of
193 ``K``. Moreover, any nonnegative linear combination of these
194 matrices shares the same property.
198 The trivial cone in a trivial space has no positive operators::
200 sage: K = Cone([], ToricLattice(0))
201 sage: positive_operator_gens(K)
204 Positive operators on the nonnegative orthant are nonnegative matrices::
206 sage: K = Cone([(1,)])
207 sage: positive_operator_gens(K)
210 sage: K = Cone([(1,0),(0,1)])
211 sage: positive_operator_gens(K)
213 [1 0] [0 1] [0 0] [0 0]
214 [0 0], [0 0], [1 0], [0 1]
217 Every operator is positive on the ambient vector space::
219 sage: K = Cone([(1,),(-1,)])
220 sage: K.is_full_space()
222 sage: positive_operator_gens(K)
225 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
226 sage: K.is_full_space()
228 sage: positive_operator_gens(K)
230 [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0]
231 [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1]
236 A positive operator on a cone should send its generators into the cone::
238 sage: set_random_seed()
239 sage: K = random_cone(max_ambient_dim=5)
240 sage: pi_of_K = positive_operator_gens(K)
241 sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
244 The dimension of the cone of positive operators is given by the
245 corollary in my paper::
247 sage: set_random_seed()
248 sage: K = random_cone(max_ambient_dim=5)
249 sage: n = K.lattice_dim()
251 sage: l = K.lineality()
252 sage: pi_of_K = positive_operator_gens(K)
253 sage: L = ToricLattice(n**2)
254 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
255 sage: expected = n**2 - l*(m - l) - (n - m)*m
256 sage: actual == expected
259 The lineality of the cone of positive operators is given by the
260 corollary in my paper::
262 sage: set_random_seed()
263 sage: K = random_cone(max_ambient_dim=5)
264 sage: n = K.lattice_dim()
265 sage: pi_of_K = positive_operator_gens(K)
266 sage: L = ToricLattice(n**2)
267 sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
268 sage: expected = n**2 - K.dim()*K.dual().dim()
269 sage: actual == expected
272 The cone ``K`` is proper if and only if the cone of positive
273 operators on ``K`` is proper::
275 sage: set_random_seed()
276 sage: K = random_cone(max_ambient_dim=5)
277 sage: pi_of_K = positive_operator_gens(K)
278 sage: L = ToricLattice(K.lattice_dim()**2)
279 sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
280 sage: K.is_proper() == pi_cone.is_proper()
283 # Matrices are not vectors in Sage, so we have to convert them
284 # to vectors explicitly before we can find a basis. We need these
285 # two values to construct the appropriate "long vector" space.
286 F
= K
.lattice().base_field()
289 tensor_products
= [ s
.tensor_product(x
) for x
in K
for s
in K
.dual() ]
291 # Convert those tensor products to long vectors.
292 W
= VectorSpace(F
, n
**2)
293 vectors
= [ W(tp
.list()) for tp
in tensor_products
]
295 # Create the *dual* cone of the positive operators, expressed as
297 pi_dual
= Cone(vectors
, ToricLattice(W
.dimension()))
299 # Now compute the desired cone from its dual...
300 pi_cone
= pi_dual
.dual()
302 # And finally convert its rays back to matrix representations.
303 M
= MatrixSpace(F
, n
)
304 return [ M(v
.list()) for v
in pi_cone
.rays() ]
307 def Z_transformation_gens(K
):
309 Compute generators of the cone of Z-transformations on this cone.
313 A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
314 Each matrix ``L`` in the list should have the property that
315 ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
316 discrete complementarity set of ``K``. Moreover, any nonnegative
317 linear combination of these matrices shares the same property.
321 Z-transformations on the nonnegative orthant are just Z-matrices.
322 That is, matrices whose off-diagonal elements are nonnegative::
324 sage: K = Cone([(1,0),(0,1)])
325 sage: Z_transformation_gens(K)
327 [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0]
328 [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1]
330 sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
331 sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
332 ....: for i in range(z.nrows())
333 ....: for j in range(z.ncols())
337 The trivial cone in a trivial space has no Z-transformations::
339 sage: K = Cone([], ToricLattice(0))
340 sage: Z_transformation_gens(K)
343 Z-transformations on a subspace are Lyapunov-like and vice-versa::
345 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
346 sage: K.is_full_space()
348 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
349 sage: zs = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
355 The Z-property is possessed by every Z-transformation::
357 sage: set_random_seed()
358 sage: K = random_cone(max_ambient_dim=6)
359 sage: Z_of_K = Z_transformation_gens(K)
360 sage: dcs = K.discrete_complementarity_set()
361 sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
362 ....: for (x,s) in dcs])
365 The lineality space of Z is LL::
367 sage: set_random_seed()
368 sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
369 sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
370 sage: z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ])
371 sage: z_cone.linear_subspace() == lls
374 And thus, the lineality of Z is the Lyapunov rank::
376 sage: set_random_seed()
377 sage: K = random_cone(max_ambient_dim=6)
378 sage: Z_of_K = Z_transformation_gens(K)
379 sage: L = ToricLattice(K.lattice_dim()**2)
380 sage: z_cone = Cone([ z.list() for z in Z_of_K ], lattice=L)
381 sage: z_cone.lineality() == K.lyapunov_rank()
384 The lineality spaces of pi-star and Z-star are equal:
386 sage: set_random_seed()
387 sage: K = random_cone(max_ambient_dim=5)
388 sage: pi_of_K = positive_operator_gens(K)
389 sage: Z_of_K = Z_transformation_gens(K)
390 sage: L = ToricLattice(K.lattice_dim()**2)
391 sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
392 sage: z_star = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
393 sage: pi_star.linear_subspace() == z_star.linear_subspace()
396 # Matrices are not vectors in Sage, so we have to convert them
397 # to vectors explicitly before we can find a basis. We need these
398 # two values to construct the appropriate "long vector" space.
399 F
= K
.lattice().base_field()
402 # These tensor products contain generators for the dual cone of
403 # the cross-positive transformations.
404 tensor_products
= [ s
.tensor_product(x
)
405 for (x
,s
) in K
.discrete_complementarity_set() ]
407 # Turn our matrices into long vectors...
408 W
= VectorSpace(F
, n
**2)
409 vectors
= [ W(m
.list()) for m
in tensor_products
]
411 # Create the *dual* cone of the cross-positive operators,
412 # expressed as long vectors..
413 Sigma_dual
= Cone(vectors
, lattice
=ToricLattice(W
.dimension()))
415 # Now compute the desired cone from its dual...
416 Sigma_cone
= Sigma_dual
.dual()
418 # And finally convert its rays back to matrix representations.
419 # But first, make them negative, so we get Z-transformations and
420 # not cross-positive ones.
421 M
= MatrixSpace(F
, n
)
422 return [ -M(v
.list()) for v
in Sigma_cone
.rays() ]
426 gens
= Z_transformation_gens(K
)
430 return Cone([ g
.list() for g
in gens
], lattice
=L
)
433 gens
= positive_operator_gens(K
)
437 return Cone([ g
.list() for g
in gens
], lattice
=L
)