mjo/cone/cone.py

   1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
   2 # have to explicitly mangle our sitedir here so that "mjo.cone"
   3 # resolves.
   4 from os.path import abspath
   5 from site import addsitedir
   6 addsitedir(abspath('../../'))
   7
   8 from sage.all import *
   9
  10 def is_lyapunov_like(L,K):
  11     r"""
  12     Determine whether or not ``L`` is Lyapunov-like on ``K``.
  13
  14     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
  15     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
  16     `\left\langle x,s \right\rangle` in the complementarity set of
  17     ``K``. It is known [Orlitzky]_ that this property need only be
  18     checked for generators of ``K`` and its dual.
  19
  20     INPUT:
  21
  22     - ``L`` -- A linear transformation or matrix.
  23
  24     - ``K`` -- A polyhedral closed convex cone.
  25
  26     OUTPUT:
  27
  28     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
  29     and ``False`` otherwise.
  30
  31     .. WARNING::
  32
  33         If this function returns ``True``, then ``L`` is Lyapunov-like
  34         on ``K``. However, if ``False`` is returned, that could mean one
  35         of two things. The first is that ``L`` is definitely not
  36         Lyapunov-like on ``K``. The second is more of an "I don't know"
  37         answer, returned (for example) if we cannot prove that an inner
  38         product is zero.
  39
  40     REFERENCES:
  41
  42     M. Orlitzky. The Lyapunov rank of an improper cone.
  43     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
  44
  45     EXAMPLES:
  46
  47     The identity is always Lyapunov-like in a nontrivial space::
  48
  49         sage: set_random_seed()
  50         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
  51         sage: L = identity_matrix(K.lattice_dim())
  52         sage: is_lyapunov_like(L,K)
  53         True
  54
  55     As is the "zero" transformation::
  56
  57         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
  58         sage: R = K.lattice().vector_space().base_ring()
  59         sage: L = zero_matrix(R, K.lattice_dim())
  60         sage: is_lyapunov_like(L,K)
  61         True
  62
  63         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
  64         on ``K``::
  65
  66         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
  67         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
  68         True
  69
  70     """
  71     return all([(L*x).inner_product(s) == 0
  72                 for (x,s) in K.discrete_complementarity_set()])
  73
  74
  75 def random_element(K):
  76     r"""
  77     Return a random element of ``K`` from its ambient vector space.
  78
  79     ALGORITHM:
  80
  81     The cone ``K`` is specified in terms of its generators, so that
  82     ``K`` is equal to the convex conic combination of those generators.
  83     To choose a random element of ``K``, we assign random nonnegative
  84     coefficients to each generator of ``K`` and construct a new vector
  85     from the scaled rays.
  86
  87     A vector, rather than a ray, is returned so that the element may
  88     have non-integer coordinates. Thus the element may have an
  89     arbitrarily small norm.
  90
  91     EXAMPLES:
  92
  93     A random element of the trivial cone is zero::
  94
  95         sage: set_random_seed()
  96         sage: K = Cone([], ToricLattice(0))
  97         sage: random_element(K)
  98         ()
  99         sage: K = Cone([(0,)])
 100         sage: random_element(K)
 101         (0)
 102         sage: K = Cone([(0,0)])
 103         sage: random_element(K)
 104         (0, 0)
 105         sage: K = Cone([(0,0,0)])
 106         sage: random_element(K)
 107         (0, 0, 0)
 108
 109     TESTS:
 110
 111     Any cone should contain an element of itself::
 112
 113         sage: set_random_seed()
 114         sage: K = random_cone(max_rays = 8)
 115         sage: K.contains(random_element(K))
 116         True
 117
 118     """
 119     V = K.lattice().vector_space()
 120     F = V.base_ring()
 121     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 122     vector_gens  = map(V, K.rays())
 123     scaled_gens  = [ coefficients[i]*vector_gens[i]
 124                          for i in range(len(vector_gens)) ]
 125
 126     # Make sure we return a vector. Without the coercion, we might
 127     # return ``0`` when ``K`` has no rays.
 128     v = V(sum(scaled_gens))
 129     return v
 130
 131
 132 def positive_operators(K):
 133     r"""
 134     Compute generators of the cone of positive operators on this cone.
 135
 136     OUTPUT:
 137
 138     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 139     Each matrix ``P`` in the list should have the property that ``P*x``
 140     is an element of ``K`` whenever ``x`` is an element of
 141     ``K``. Moreover, any nonnegative linear combination of these
 142     matrices shares the same property.
 143
 144     EXAMPLES:
 145
 146     The trivial cone in a trivial space has no positive operators::
 147
 148         sage: K = Cone([], ToricLattice(0))
 149         sage: positive_operators(K)
 150         []
 151
 152     Positive operators on the nonnegative orthant are nonnegative matrices::
 153
 154         sage: K = Cone([(1,)])
 155         sage: positive_operators(K)
 156         [[1]]
 157
 158         sage: K = Cone([(1,0),(0,1)])
 159         sage: positive_operators(K)
 160         [
 161         [1 0]  [0 1]  [0 0]  [0 0]
 162         [0 0], [0 0], [1 0], [0 1]
 163         ]
 164
 165     Every operator is positive on the ambient vector space::
 166
 167         sage: K = Cone([(1,),(-1,)])
 168         sage: K.is_full_space()
 169         True
 170         sage: positive_operators(K)
 171         [[1], [-1]]
 172
 173         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 174         sage: K.is_full_space()
 175         True
 176         sage: positive_operators(K)
 177         [
 178         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 179         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 180         ]
 181
 182     TESTS:
 183
 184     A positive operator on a cone should send its generators into the cone::
 185
 186         sage: K = random_cone(max_ambient_dim = 6)
 187         sage: pi_of_K = positive_operators(K)
 188         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 189         True
 190
 191     """
 192     # Sage doesn't think matrices are vectors, so we have to convert
 193     # our matrices to vectors explicitly before we can figure out how
 194     # many are linearly-indepenedent.
 195     #
 196     # The space W has the same base ring as V, but dimension
 197     # dim(V)^2. So it has the same dimension as the space of linear
 198     # transformations on V. In other words, it's just the right size
 199     # to create an isomorphism between it and our matrices.
 200     V = K.lattice().vector_space()
 201     W = VectorSpace(V.base_ring(), V.dimension()**2)
 202
 203     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 204
 205     # Turn our matrices into long vectors...
 206     vectors = [ W(m.list()) for m in tensor_products ]
 207
 208     # Create the *dual* cone of the positive operators, expressed as
 209     # long vectors..
 210     L = ToricLattice(W.dimension())
 211     pi_dual = Cone(vectors, lattice=L)
 212
 213     # Now compute the desired cone from its dual...
 214     pi_cone = pi_dual.dual()
 215
 216     # And finally convert its rays back to matrix representations.
 217     M = MatrixSpace(V.base_ring(), V.dimension())
 218
 219     return [ M(v.list()) for v in pi_cone.rays() ]
 220
 221
 222 def Z_transformations(K):
 223     r"""
 224     Compute generators of the cone of Z-transformations on this cone.
 225
 226     OUTPUT:
 227
 228     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 229     Each matrix ``L`` in the list should have the property that
 230     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 231     discrete complementarity set of ``K``. Moreover, any nonnegative
 232     linear combination of these matrices shares the same property.
 233
 234     EXAMPLES:
 235
 236     Z-transformations on the nonnegative orthant are just Z-matrices.
 237     That is, matrices whose off-diagonal elements are nonnegative::
 238
 239         sage: K = Cone([(1,0),(0,1)])
 240         sage: Z_transformations(K)
 241         [
 242         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 243         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 244         ]
 245         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 246         sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
 247         ....:                    for i in range(z.nrows())
 248         ....:                    for j in range(z.ncols())
 249         ....:                    if i != j ])
 250         True
 251
 252     The trivial cone in a trivial space has no Z-transformations::
 253
 254         sage: K = Cone([], ToricLattice(0))
 255         sage: Z_transformations(K)
 256         []
 257
 258     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 259
 260         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 261         sage: K.is_full_space()
 262         True
 263         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 264         sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
 265         sage: zs == lls
 266         True
 267
 268     TESTS:
 269
 270     The Z-property is possessed by every Z-transformation::
 271
 272         sage: set_random_seed()
 273         sage: K = random_cone(max_ambient_dim = 6)
 274         sage: Z_of_K = Z_transformations(K)
 275         sage: dcs = K.discrete_complementarity_set()
 276         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 277         ....:                                  for (x,s) in dcs])
 278         True
 279
 280     The lineality space of Z is LL::
 281
 282         sage: set_random_seed()
 283         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 284         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 285         sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
 286         sage: z_cone.linear_subspace() == lls
 287         True
 288
 289     """
 290     # Sage doesn't think matrices are vectors, so we have to convert
 291     # our matrices to vectors explicitly before we can figure out how
 292     # many are linearly-indepenedent.
 293     #
 294     # The space W has the same base ring as V, but dimension
 295     # dim(V)^2. So it has the same dimension as the space of linear
 296     # transformations on V. In other words, it's just the right size
 297     # to create an isomorphism between it and our matrices.
 298     V = K.lattice().vector_space()
 299     W = VectorSpace(V.base_ring(), V.dimension()**2)
 300
 301     C_of_K = K.discrete_complementarity_set()
 302     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 303
 304     # Turn our matrices into long vectors...
 305     vectors = [ W(m.list()) for m in tensor_products ]
 306
 307     # Create the *dual* cone of the cross-positive operators,
 308     # expressed as long vectors..
 309     L = ToricLattice(W.dimension())
 310     Sigma_dual = Cone(vectors, lattice=L)
 311
 312     # Now compute the desired cone from its dual...
 313     Sigma_cone = Sigma_dual.dual()
 314
 315     # And finally convert its rays back to matrix representations.
 316     # But first, make them negative, so we get Z-transformations and
 317     # not cross-positive ones.
 318     M = MatrixSpace(V.base_ring(), V.dimension())
 319
 320     return [ -M(v.list()) for v in Sigma_cone.rays() ]