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_rays = 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_rays = 5)
  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_rays = 5)
  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_operators(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_operators(K)
 143         []
 144
 145     Positive operators on the nonnegative orthant are nonnegative matrices::
 146
 147         sage: K = Cone([(1,)])
 148         sage: positive_operators(K)
 149         [[1]]
 150
 151         sage: K = Cone([(1,0),(0,1)])
 152         sage: positive_operators(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_operators(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_operators(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: K = random_cone(max_ambient_dim = 6)
 180         sage: pi_of_K = positive_operators(K)
 181         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 182         True
 183
 184     """
 185     # Sage doesn't think matrices are vectors, so we have to convert
 186     # our matrices to vectors explicitly before we can figure out how
 187     # many are linearly-indepenedent.
 188     #
 189     # The space W has the same base ring as V, but dimension
 190     # dim(V)^2. So it has the same dimension as the space of linear
 191     # transformations on V. In other words, it's just the right size
 192     # to create an isomorphism between it and our matrices.
 193     V = K.lattice().vector_space()
 194     W = VectorSpace(V.base_ring(), V.dimension()**2)
 195
 196     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 197
 198     # Turn our matrices into long vectors...
 199     vectors = [ W(m.list()) for m in tensor_products ]
 200
 201     # Create the *dual* cone of the positive operators, expressed as
 202     # long vectors..
 203     L = ToricLattice(W.dimension())
 204     pi_dual = Cone(vectors, lattice=L)
 205
 206     # Now compute the desired cone from its dual...
 207     pi_cone = pi_dual.dual()
 208
 209     # And finally convert its rays back to matrix representations.
 210     M = MatrixSpace(V.base_ring(), V.dimension())
 211
 212     return [ M(v.list()) for v in pi_cone.rays() ]
 213
 214
 215 def Z_transformations(K):
 216     r"""
 217     Compute generators of the cone of Z-transformations on this cone.
 218
 219     OUTPUT:
 220
 221     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 222     Each matrix ``L`` in the list should have the property that
 223     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 224     discrete complementarity set of ``K``. Moreover, any nonnegative
 225     linear combination of these matrices shares the same property.
 226
 227     EXAMPLES:
 228
 229     Z-transformations on the nonnegative orthant are just Z-matrices.
 230     That is, matrices whose off-diagonal elements are nonnegative::
 231
 232         sage: K = Cone([(1,0),(0,1)])
 233         sage: Z_transformations(K)
 234         [
 235         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 236         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 237         ]
 238         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 239         sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
 240         ....:                    for i in range(z.nrows())
 241         ....:                    for j in range(z.ncols())
 242         ....:                    if i != j ])
 243         True
 244
 245     The trivial cone in a trivial space has no Z-transformations::
 246
 247         sage: K = Cone([], ToricLattice(0))
 248         sage: Z_transformations(K)
 249         []
 250
 251     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 252
 253         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 254         sage: K.is_full_space()
 255         True
 256         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 257         sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
 258         sage: zs == lls
 259         True
 260
 261     TESTS:
 262
 263     The Z-property is possessed by every Z-transformation::
 264
 265         sage: set_random_seed()
 266         sage: K = random_cone(max_ambient_dim = 6)
 267         sage: Z_of_K = Z_transformations(K)
 268         sage: dcs = K.discrete_complementarity_set()
 269         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 270         ....:                                  for (x,s) in dcs])
 271         True
 272
 273     The lineality space of Z is LL::
 274
 275         sage: set_random_seed()
 276         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 277         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 278         sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
 279         sage: z_cone.linear_subspace() == lls
 280         True
 281
 282     """
 283     # Sage doesn't think matrices are vectors, so we have to convert
 284     # our matrices to vectors explicitly before we can figure out how
 285     # many are linearly-indepenedent.
 286     #
 287     # The space W has the same base ring as V, but dimension
 288     # dim(V)^2. So it has the same dimension as the space of linear
 289     # transformations on V. In other words, it's just the right size
 290     # to create an isomorphism between it and our matrices.
 291     V = K.lattice().vector_space()
 292     W = VectorSpace(V.base_ring(), V.dimension()**2)
 293
 294     C_of_K = K.discrete_complementarity_set()
 295     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 296
 297     # Turn our matrices into long vectors...
 298     vectors = [ W(m.list()) for m in tensor_products ]
 299
 300     # Create the *dual* cone of the cross-positive operators,
 301     # expressed as long vectors..
 302     L = ToricLattice(W.dimension())
 303     Sigma_dual = Cone(vectors, lattice=L)
 304
 305     # Now compute the desired cone from its dual...
 306     Sigma_cone = Sigma_dual.dual()
 307
 308     # And finally convert its rays back to matrix representations.
 309     # But first, make them negative, so we get Z-transformations and
 310     # not cross-positive ones.
 311     M = MatrixSpace(V.base_ring(), V.dimension())
 312
 313     return [ -M(v.list()) for v in Sigma_cone.rays() ]