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: K = random_cone(max_ambient_dim = 6)
 180         sage: pi_of_K = positive_operator_gens(K)
 181         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 182         True
 183
 184     The dimension of the cone of positive operators is given by the
 185     corollary in my paper::
 186
 187         sage: K = random_cone(max_ambient_dim = 6)
 188         sage: n = K.lattice_dim()
 189         sage: m = K.dim()
 190         sage: l = K.lineality()
 191         sage: pi_of_K = positive_operator_gens(K)
 192         sage: actual = Cone([p.list() for p in pi_of_K]).dim()
 193         sage: expected = n**2 - l*(n - l) - (n - m)*m
 194         sage: actual == expected
 195         True
 196
 197     """
 198     # Matrices are not vectors in Sage, so we have to convert them
 199     # to vectors explicitly before we can find a basis. We need these
 200     # two values to construct the appropriate "long vector" space.
 201     F = K.lattice().base_field()
 202     n = K.lattice_dim()
 203
 204     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 205
 206     # Convert those tensor products to long vectors.
 207     W = VectorSpace(F, n**2)
 208     vectors = [ W(tp.list()) for tp in tensor_products ]
 209
 210     # Create the *dual* cone of the positive operators, expressed as
 211     # long vectors..
 212     pi_dual = Cone(vectors, ToricLattice(W.dimension()))
 213
 214     # Now compute the desired cone from its dual...
 215     pi_cone = pi_dual.dual()
 216
 217     # And finally convert its rays back to matrix representations.
 218     M = MatrixSpace(F, n)
 219     return [ M(v.list()) for v in pi_cone.rays() ]
 220
 221
 222 def Z_transformation_gens(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_transformation_gens(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_transformation_gens(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_transformation_gens(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_transformation_gens(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_transformation_gens(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_transformation_gens(K) ])
 286         sage: z_cone.linear_subspace() == lls
 287         True
 288
 289     """
 290     # Matrices are not vectors in Sage, so we have to convert them
 291     # to vectors explicitly before we can find a basis. We need these
 292     # two values to construct the appropriate "long vector" space.
 293     F = K.lattice().base_field()
 294     n = K.lattice_dim()
 295
 296     # These tensor products contain generators for the dual cone of
 297     # the cross-positive transformations.
 298     tensor_products = [ s.tensor_product(x)
 299                         for (x,s) in K.discrete_complementarity_set() ]
 300
 301     # Turn our matrices into long vectors...
 302     W = VectorSpace(F, n**2)
 303     vectors = [ W(m.list()) for m in tensor_products ]
 304
 305     # Create the *dual* cone of the cross-positive operators,
 306     # expressed as long vectors..
 307     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
 308
 309     # Now compute the desired cone from its dual...
 310     Sigma_cone = Sigma_dual.dual()
 311
 312     # And finally convert its rays back to matrix representations.
 313     # But first, make them negative, so we get Z-transformations and
 314     # not cross-positive ones.
 315     M = MatrixSpace(F, n)
 316     return [ -M(v.list()) for v in Sigma_cone.rays() ]