mjo/cone/cone.py

   1 from sage.all import *
   2
   3 def is_lyapunov_like(L,K):
   4     r"""
   5     Determine whether or not ``L`` is Lyapunov-like on ``K``.
   6
   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.
  12
  13     INPUT:
  14
  15     - ``L`` -- A linear transformation or matrix.
  16
  17     - ``K`` -- A polyhedral closed convex cone.
  18
  19     OUTPUT:
  20
  21     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
  22     and ``False`` otherwise.
  23
  24     .. WARNING::
  25
  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
  31         product is zero.
  32
  33     REFERENCES:
  34
  35     M. Orlitzky. The Lyapunov rank of an improper cone.
  36     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
  37
  38     EXAMPLES:
  39
  40     The identity is always Lyapunov-like in a nontrivial space::
  41
  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)
  46         True
  47
  48     As is the "zero" transformation::
  49
  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)
  54         True
  55
  56         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
  57         on ``K``::
  58
  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() ])
  61         True
  62
  63     """
  64     return all([(L*x).inner_product(s) == 0
  65                 for (x,s) in K.discrete_complementarity_set()])
  66
  67
  68 def motzkin_decomposition(K):
  69     r"""
  70     Return the pair of components in the motzkin decomposition of this cone.
  71
  72     Every convex cone is the direct sum of a strictly convex cone and a
  73     linear subspace. Return a pair ``(P,S)`` of cones such that ``P`` is
  74     strictly convex, ``S`` is a subspace, and ``K`` is the direct sum of
  75     ``P`` and ``S``.
  76
  77     OUTPUT:
  78
  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``.
  82
  83     EXAMPLES:
  84
  85     The nonnegative orthant is strictly convex, so it is its own
  86     strictly convex component and its subspace component is trivial::
  87
  88         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  89         sage: (P,S) = motzkin_decomposition(K)
  90         sage: K.is_equivalent(P)
  91         True
  92         sage: S.is_trivial()
  93         True
  94
  95     Likewise, full spaces are their own subspace components::
  96
  97         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
  98         sage: K.is_full_space()
  99         True
 100         sage: (P,S) = motzkin_decomposition(K)
 101         sage: K.is_equivalent(S)
 102         True
 103         sage: P.is_trivial()
 104         True
 105
 106     TESTS:
 107
 108     A random point in the cone should belong to either the strictly
 109     convex component or the subspace component. If the point is nonzero,
 110     it cannot be in both::
 111
 112         sage: set_random_seed()
 113         sage: K = random_cone(max_ambient_dim=8)
 114         sage: (P,S) = motzkin_decomposition(K)
 115         sage: x = K.random_element()
 116         sage: P.contains(x) or S.contains(x)
 117         True
 118         sage: x.is_zero() or (P.contains(x) != S.contains(x))
 119         True
 120
 121     The strictly convex component should always be strictly convex, and
 122     the subspace component should always be a subspace::
 123
 124         sage: set_random_seed()
 125         sage: K = random_cone(max_ambient_dim=8)
 126         sage: (P,S) = motzkin_decomposition(K)
 127         sage: P.is_strictly_convex()
 128         True
 129         sage: S.lineality() == S.dim()
 130         True
 131
 132     The generators of the strictly convex component are obtained from
 133     the orthogonal projections of the original generators onto the
 134     orthogonal complement of the subspace component::
 135
 136         sage: set_random_seed()
 137         sage: K = random_cone(max_ambient_dim=8)
 138         sage: (P,S) = motzkin_decomposition(K)
 139         sage: S_perp = S.linear_subspace().complement()
 140         sage: A = S_perp.matrix().transpose()
 141         sage: proj = A * (A.transpose()*A).inverse() * A.transpose()
 142         sage: expected = Cone([ proj*g for g in K ], K.lattice())
 143         sage: P.is_equivalent(expected)
 144         True
 145     """
 146     linspace_gens  = [ copy(b) for b in K.linear_subspace().basis() ]
 147     linspace_gens += [ -b for b in linspace_gens ]
 148
 149     S = Cone(linspace_gens, K.lattice())
 150
 151     # Since ``S`` is a subspace, its dual is its orthogonal complement
 152     # (albeit in the wrong lattice).
 153     S_perp = Cone(S.dual(), K.lattice())
 154     P = K.intersection(S_perp)
 155
 156     return (P,S)
 157
 158 def positive_operator_gens(K):
 159     r"""
 160     Compute generators of the cone of positive operators on this cone.
 161
 162     OUTPUT:
 163
 164     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 165     Each matrix ``P`` in the list should have the property that ``P*x``
 166     is an element of ``K`` whenever ``x`` is an element of
 167     ``K``. Moreover, any nonnegative linear combination of these
 168     matrices shares the same property.
 169
 170     EXAMPLES:
 171
 172     The trivial cone in a trivial space has no positive operators::
 173
 174         sage: K = Cone([], ToricLattice(0))
 175         sage: positive_operator_gens(K)
 176         []
 177
 178     Positive operators on the nonnegative orthant are nonnegative matrices::
 179
 180         sage: K = Cone([(1,)])
 181         sage: positive_operator_gens(K)
 182         [[1]]
 183
 184         sage: K = Cone([(1,0),(0,1)])
 185         sage: positive_operator_gens(K)
 186         [
 187         [1 0]  [0 1]  [0 0]  [0 0]
 188         [0 0], [0 0], [1 0], [0 1]
 189         ]
 190
 191     Every operator is positive on the ambient vector space::
 192
 193         sage: K = Cone([(1,),(-1,)])
 194         sage: K.is_full_space()
 195         True
 196         sage: positive_operator_gens(K)
 197         [[1], [-1]]
 198
 199         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 200         sage: K.is_full_space()
 201         True
 202         sage: positive_operator_gens(K)
 203         [
 204         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 205         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 206         ]
 207
 208     TESTS:
 209
 210     A positive operator on a cone should send its generators into the cone::
 211
 212         sage: set_random_seed()
 213         sage: K = random_cone(max_ambient_dim=5)
 214         sage: pi_of_K = positive_operator_gens(K)
 215         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 216         True
 217
 218     The dimension of the cone of positive operators is given by the
 219     corollary in my paper::
 220
 221         sage: set_random_seed()
 222         sage: K = random_cone(max_ambient_dim=5)
 223         sage: n = K.lattice_dim()
 224         sage: m = K.dim()
 225         sage: l = K.lineality()
 226         sage: pi_of_K = positive_operator_gens(K)
 227         sage: L = ToricLattice(n**2)
 228         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
 229         sage: expected = n**2 - l*(m - l) - (n - m)*m
 230         sage: actual == expected
 231         True
 232
 233     The lineality of the cone of positive operators is given by the
 234     corollary in my paper::
 235
 236         sage: set_random_seed()
 237         sage: K = random_cone(max_ambient_dim=5)
 238         sage: n = K.lattice_dim()
 239         sage: pi_of_K = positive_operator_gens(K)
 240         sage: L = ToricLattice(n**2)
 241         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
 242         sage: expected = n**2 - K.dim()*K.dual().dim()
 243         sage: actual == expected
 244         True
 245
 246     The cone ``K`` is proper if and only if the cone of positive
 247     operators on ``K`` is proper::
 248
 249         sage: set_random_seed()
 250         sage: K = random_cone(max_ambient_dim=5)
 251         sage: pi_of_K = positive_operator_gens(K)
 252         sage: L = ToricLattice(K.lattice_dim()**2)
 253         sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
 254         sage: K.is_proper() == pi_cone.is_proper()
 255         True
 256     """
 257     # Matrices are not vectors in Sage, so we have to convert them
 258     # to vectors explicitly before we can find a basis. We need these
 259     # two values to construct the appropriate "long vector" space.
 260     F = K.lattice().base_field()
 261     n = K.lattice_dim()
 262
 263     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 264
 265     # Convert those tensor products to long vectors.
 266     W = VectorSpace(F, n**2)
 267     vectors = [ W(tp.list()) for tp in tensor_products ]
 268
 269     # Create the *dual* cone of the positive operators, expressed as
 270     # long vectors..
 271     pi_dual = Cone(vectors, ToricLattice(W.dimension()))
 272
 273     # Now compute the desired cone from its dual...
 274     pi_cone = pi_dual.dual()
 275
 276     # And finally convert its rays back to matrix representations.
 277     M = MatrixSpace(F, n)
 278     return [ M(v.list()) for v in pi_cone.rays() ]
 279
 280
 281 def Z_transformation_gens(K):
 282     r"""
 283     Compute generators of the cone of Z-transformations on this cone.
 284
 285     OUTPUT:
 286
 287     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 288     Each matrix ``L`` in the list should have the property that
 289     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 290     discrete complementarity set of ``K``. Moreover, any nonnegative
 291     linear combination of these matrices shares the same property.
 292
 293     EXAMPLES:
 294
 295     Z-transformations on the nonnegative orthant are just Z-matrices.
 296     That is, matrices whose off-diagonal elements are nonnegative::
 297
 298         sage: K = Cone([(1,0),(0,1)])
 299         sage: Z_transformation_gens(K)
 300         [
 301         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 302         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 303         ]
 304         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 305         sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
 306         ....:                    for i in range(z.nrows())
 307         ....:                    for j in range(z.ncols())
 308         ....:                    if i != j ])
 309         True
 310
 311     The trivial cone in a trivial space has no Z-transformations::
 312
 313         sage: K = Cone([], ToricLattice(0))
 314         sage: Z_transformation_gens(K)
 315         []
 316
 317     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 318
 319         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 320         sage: K.is_full_space()
 321         True
 322         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 323         sage: zs  = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
 324         sage: zs == lls
 325         True
 326
 327     TESTS:
 328
 329     The Z-property is possessed by every Z-transformation::
 330
 331         sage: set_random_seed()
 332         sage: K = random_cone(max_ambient_dim=6)
 333         sage: Z_of_K = Z_transformation_gens(K)
 334         sage: dcs = K.discrete_complementarity_set()
 335         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 336         ....:                                  for (x,s) in dcs])
 337         True
 338
 339     The lineality space of Z is LL::
 340
 341         sage: set_random_seed()
 342         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 343         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 344         sage: z_cone  = Cone([ z.list() for z in Z_transformation_gens(K) ])
 345         sage: z_cone.linear_subspace() == lls
 346         True
 347
 348     And thus, the lineality of Z is the Lyapunov rank::
 349
 350         sage: set_random_seed()
 351         sage: K = random_cone(max_ambient_dim=6)
 352         sage: Z_of_K = Z_transformation_gens(K)
 353         sage: L = ToricLattice(K.lattice_dim()**2)
 354         sage: z_cone  = Cone([ z.list() for z in Z_of_K ], lattice=L)
 355         sage: z_cone.lineality() == K.lyapunov_rank()
 356         True
 357
 358     The lineality spaces of pi-star and Z-star are equal:
 359
 360         sage: set_random_seed()
 361         sage: K = random_cone(max_ambient_dim=5)
 362         sage: pi_of_K = positive_operator_gens(K)
 363         sage: Z_of_K = Z_transformation_gens(K)
 364         sage: L = ToricLattice(K.lattice_dim()**2)
 365         sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
 366         sage: z_star  = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
 367         sage: pi_star.linear_subspace() == z_star.linear_subspace()
 368         True
 369     """
 370     # Matrices are not vectors in Sage, so we have to convert them
 371     # to vectors explicitly before we can find a basis. We need these
 372     # two values to construct the appropriate "long vector" space.
 373     F = K.lattice().base_field()
 374     n = K.lattice_dim()
 375
 376     # These tensor products contain generators for the dual cone of
 377     # the cross-positive transformations.
 378     tensor_products = [ s.tensor_product(x)
 379                         for (x,s) in K.discrete_complementarity_set() ]
 380
 381     # Turn our matrices into long vectors...
 382     W = VectorSpace(F, n**2)
 383     vectors = [ W(m.list()) for m in tensor_products ]
 384
 385     # Create the *dual* cone of the cross-positive operators,
 386     # expressed as long vectors..
 387     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
 388
 389     # Now compute the desired cone from its dual...
 390     Sigma_cone = Sigma_dual.dual()
 391
 392     # And finally convert its rays back to matrix representations.
 393     # But first, make them negative, so we get Z-transformations and
 394     # not cross-positive ones.
 395     M = MatrixSpace(F, n)
 396     return [ -M(v.list()) for v in Sigma_cone.rays() ]
 397
 398
 399 def Z_cone(K):
 400     gens = Z_transformation_gens(K)
 401     L = None
 402     if len(gens) == 0:
 403         L = ToricLattice(0)
 404     return Cone([ g.list() for g in gens ], lattice=L)
 405
 406 def pi_cone(K):
 407     gens = positive_operator_gens(K)
 408     L = None
 409     if len(gens) == 0:
 410         L = ToricLattice(0)
 411     return Cone([ g.list() for g in gens ], lattice=L)