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 random_element(K):
  69     r"""
  70     Return a random element of ``K`` from its ambient vector space.
  71
  72     ALGORITHM:
  73
  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
  78     from the scaled rays.
  79
  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.
  83
  84     EXAMPLES:
  85
  86     A random element of the trivial cone is zero::
  87
  88         sage: set_random_seed()
  89         sage: K = Cone([], ToricLattice(0))
  90         sage: random_element(K)
  91         ()
  92         sage: K = Cone([(0,)])
  93         sage: random_element(K)
  94         (0)
  95         sage: K = Cone([(0,0)])
  96         sage: random_element(K)
  97         (0, 0)
  98         sage: K = Cone([(0,0,0)])
  99         sage: random_element(K)
 100         (0, 0, 0)
 101
 102     TESTS:
 103
 104     Any cone should contain an element of itself::
 105
 106         sage: set_random_seed()
 107         sage: K = random_cone(max_rays = 8)
 108         sage: K.contains(random_element(K))
 109         True
 110
 111     """
 112     V = K.lattice().vector_space()
 113     F = V.base_ring()
 114     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 115     vector_gens  = map(V, K.rays())
 116     scaled_gens  = [ coefficients[i]*vector_gens[i]
 117                          for i in range(len(vector_gens)) ]
 118
 119     # Make sure we return a vector. Without the coercion, we might
 120     # return ``0`` when ``K`` has no rays.
 121     v = V(sum(scaled_gens))
 122     return v
 123
 124
 125 def positive_operator_gens(K):
 126     r"""
 127     Compute generators of the cone of positive operators on this cone.
 128
 129     OUTPUT:
 130
 131     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 132     Each matrix ``P`` in the list should have the property that ``P*x``
 133     is an element of ``K`` whenever ``x`` is an element of
 134     ``K``. Moreover, any nonnegative linear combination of these
 135     matrices shares the same property.
 136
 137     EXAMPLES:
 138
 139     The trivial cone in a trivial space has no positive operators::
 140
 141         sage: K = Cone([], ToricLattice(0))
 142         sage: positive_operator_gens(K)
 143         []
 144
 145     Positive operators on the nonnegative orthant are nonnegative matrices::
 146
 147         sage: K = Cone([(1,)])
 148         sage: positive_operator_gens(K)
 149         [[1]]
 150
 151         sage: K = Cone([(1,0),(0,1)])
 152         sage: positive_operator_gens(K)
 153         [
 154         [1 0]  [0 1]  [0 0]  [0 0]
 155         [0 0], [0 0], [1 0], [0 1]
 156         ]
 157
 158     Every operator is positive on the ambient vector space::
 159
 160         sage: K = Cone([(1,),(-1,)])
 161         sage: K.is_full_space()
 162         True
 163         sage: positive_operator_gens(K)
 164         [[1], [-1]]
 165
 166         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 167         sage: K.is_full_space()
 168         True
 169         sage: positive_operator_gens(K)
 170         [
 171         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 172         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 173         ]
 174
 175     TESTS:
 176
 177     A positive operator on a cone should send its generators into the cone::
 178
 179         sage: set_random_seed()
 180         sage: K = random_cone(max_ambient_dim=5)
 181         sage: pi_of_K = positive_operator_gens(K)
 182         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 183         True
 184
 185     The dimension of the cone of positive operators is given by the
 186     corollary in my paper::
 187
 188         sage: set_random_seed()
 189         sage: K = random_cone(max_ambient_dim = 5)
 190         sage: n = K.lattice_dim()
 191         sage: m = K.dim()
 192         sage: l = K.lineality()
 193         sage: pi_of_K = positive_operator_gens(K)
 194         sage: L = ToricLattice(n**2)
 195         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
 196         sage: expected = n**2 - l*(m - l) - (n - m)*m
 197         sage: actual == expected
 198         True
 199
 200     The lineality of the cone of positive operators is given by the
 201     corollary in my paper::
 202
 203         sage: set_random_seed()
 204         sage: K = random_cone(max_ambient_dim=5)
 205         sage: n = K.lattice_dim()
 206         sage: pi_of_K = positive_operator_gens(K)
 207         sage: L = ToricLattice(n**2)
 208         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
 209         sage: expected = n**2 - K.dim()*K.dual().dim()
 210         sage: actual == expected
 211         True
 212
 213     The cone ``K`` is proper if and only if the cone of positive
 214     operators on ``K`` is proper::
 215
 216         sage: set_random_seed()
 217         sage: K = random_cone(max_ambient_dim=5)
 218         sage: pi_of_K = positive_operator_gens(K)
 219         sage: L = ToricLattice(K.lattice_dim()**2)
 220         sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
 221         sage: K.is_proper() == pi_cone.is_proper()
 222         True
 223     """
 224     # Matrices are not vectors in Sage, so we have to convert them
 225     # to vectors explicitly before we can find a basis. We need these
 226     # two values to construct the appropriate "long vector" space.
 227     F = K.lattice().base_field()
 228     n = K.lattice_dim()
 229
 230     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 231
 232     # Convert those tensor products to long vectors.
 233     W = VectorSpace(F, n**2)
 234     vectors = [ W(tp.list()) for tp in tensor_products ]
 235
 236     # Create the *dual* cone of the positive operators, expressed as
 237     # long vectors..
 238     pi_dual = Cone(vectors, ToricLattice(W.dimension()))
 239
 240     # Now compute the desired cone from its dual...
 241     pi_cone = pi_dual.dual()
 242
 243     # And finally convert its rays back to matrix representations.
 244     M = MatrixSpace(F, n)
 245     return [ M(v.list()) for v in pi_cone.rays() ]
 246
 247
 248 def Z_transformation_gens(K):
 249     r"""
 250     Compute generators of the cone of Z-transformations on this cone.
 251
 252     OUTPUT:
 253
 254     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 255     Each matrix ``L`` in the list should have the property that
 256     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 257     discrete complementarity set of ``K``. Moreover, any nonnegative
 258     linear combination of these matrices shares the same property.
 259
 260     EXAMPLES:
 261
 262     Z-transformations on the nonnegative orthant are just Z-matrices.
 263     That is, matrices whose off-diagonal elements are nonnegative::
 264
 265         sage: K = Cone([(1,0),(0,1)])
 266         sage: Z_transformation_gens(K)
 267         [
 268         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 269         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 270         ]
 271         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 272         sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
 273         ....:                    for i in range(z.nrows())
 274         ....:                    for j in range(z.ncols())
 275         ....:                    if i != j ])
 276         True
 277
 278     The trivial cone in a trivial space has no Z-transformations::
 279
 280         sage: K = Cone([], ToricLattice(0))
 281         sage: Z_transformation_gens(K)
 282         []
 283
 284     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 285
 286         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 287         sage: K.is_full_space()
 288         True
 289         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 290         sage: zs  = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
 291         sage: zs == lls
 292         True
 293
 294     TESTS:
 295
 296     The Z-property is possessed by every Z-transformation::
 297
 298         sage: set_random_seed()
 299         sage: K = random_cone(max_ambient_dim = 6)
 300         sage: Z_of_K = Z_transformation_gens(K)
 301         sage: dcs = K.discrete_complementarity_set()
 302         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 303         ....:                                  for (x,s) in dcs])
 304         True
 305
 306     The lineality space of Z is LL::
 307
 308         sage: set_random_seed()
 309         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 310         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 311         sage: z_cone  = Cone([ z.list() for z in Z_transformation_gens(K) ])
 312         sage: z_cone.linear_subspace() == lls
 313         True
 314
 315     And thus, the lineality of Z is the Lyapunov rank::
 316
 317         sage: set_random_seed()
 318         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 319         sage: z_cone  = Cone([ z.list() for z in Z_transformation_gens(K) ])
 320         sage: z_cone.lineality() == K.lyapunov_rank()
 321         True
 322     """
 323     # Matrices are not vectors in Sage, so we have to convert them
 324     # to vectors explicitly before we can find a basis. We need these
 325     # two values to construct the appropriate "long vector" space.
 326     F = K.lattice().base_field()
 327     n = K.lattice_dim()
 328
 329     # These tensor products contain generators for the dual cone of
 330     # the cross-positive transformations.
 331     tensor_products = [ s.tensor_product(x)
 332                         for (x,s) in K.discrete_complementarity_set() ]
 333
 334     # Turn our matrices into long vectors...
 335     W = VectorSpace(F, n**2)
 336     vectors = [ W(m.list()) for m in tensor_products ]
 337
 338     # Create the *dual* cone of the cross-positive operators,
 339     # expressed as long vectors..
 340     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
 341
 342     # Now compute the desired cone from its dual...
 343     Sigma_cone = Sigma_dual.dual()
 344
 345     # And finally convert its rays back to matrix representations.
 346     # But first, make them negative, so we get Z-transformations and
 347     # not cross-positive ones.
 348     M = MatrixSpace(F, n)
 349     return [ -M(v.list()) for v in Sigma_cone.rays() ]
 350
 351
 352 def Z_cone(K):
 353     gens = Z_transformation_gens(K)
 354     L = None
 355     if len(gens) == 0:
 356         L = ToricLattice(0)
 357     return Cone([ g.list() for g in gens ], lattice=L)
 358
 359 def pi_cone(K):
 360     gens = positive_operator_gens(K)
 361     L = None
 362     if len(gens) == 0:
 363         L = ToricLattice(0)
 364     return Cone([ g.list() for g in gens ], lattice=L)