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     """
 201     # Matrices are not vectors in Sage, so we have to convert them
 202     # to vectors explicitly before we can find a basis. We need these
 203     # two values to construct the appropriate "long vector" space.
 204     F = K.lattice().base_field()
 205     n = K.lattice_dim()
 206
 207     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 208
 209     # Convert those tensor products to long vectors.
 210     W = VectorSpace(F, n**2)
 211     vectors = [ W(tp.list()) for tp in tensor_products ]
 212
 213     # Create the *dual* cone of the positive operators, expressed as
 214     # long vectors..
 215     pi_dual = Cone(vectors, ToricLattice(W.dimension()))
 216
 217     # Now compute the desired cone from its dual...
 218     pi_cone = pi_dual.dual()
 219
 220     # And finally convert its rays back to matrix representations.
 221     M = MatrixSpace(F, n)
 222     return [ M(v.list()) for v in pi_cone.rays() ]
 223
 224
 225 def Z_transformation_gens(K):
 226     r"""
 227     Compute generators of the cone of Z-transformations on this cone.
 228
 229     OUTPUT:
 230
 231     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 232     Each matrix ``L`` in the list should have the property that
 233     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 234     discrete complementarity set of ``K``. Moreover, any nonnegative
 235     linear combination of these matrices shares the same property.
 236
 237     EXAMPLES:
 238
 239     Z-transformations on the nonnegative orthant are just Z-matrices.
 240     That is, matrices whose off-diagonal elements are nonnegative::
 241
 242         sage: K = Cone([(1,0),(0,1)])
 243         sage: Z_transformation_gens(K)
 244         [
 245         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 246         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 247         ]
 248         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 249         sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
 250         ....:                    for i in range(z.nrows())
 251         ....:                    for j in range(z.ncols())
 252         ....:                    if i != j ])
 253         True
 254
 255     The trivial cone in a trivial space has no Z-transformations::
 256
 257         sage: K = Cone([], ToricLattice(0))
 258         sage: Z_transformation_gens(K)
 259         []
 260
 261     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 262
 263         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 264         sage: K.is_full_space()
 265         True
 266         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 267         sage: zs  = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
 268         sage: zs == lls
 269         True
 270
 271     TESTS:
 272
 273     The Z-property is possessed by every Z-transformation::
 274
 275         sage: set_random_seed()
 276         sage: K = random_cone(max_ambient_dim = 6)
 277         sage: Z_of_K = Z_transformation_gens(K)
 278         sage: dcs = K.discrete_complementarity_set()
 279         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 280         ....:                                  for (x,s) in dcs])
 281         True
 282
 283     The lineality space of Z is LL::
 284
 285         sage: set_random_seed()
 286         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 287         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 288         sage: z_cone  = Cone([ z.list() for z in Z_transformation_gens(K) ])
 289         sage: z_cone.linear_subspace() == lls
 290         True
 291
 292     """
 293     # Matrices are not vectors in Sage, so we have to convert them
 294     # to vectors explicitly before we can find a basis. We need these
 295     # two values to construct the appropriate "long vector" space.
 296     F = K.lattice().base_field()
 297     n = K.lattice_dim()
 298
 299     # These tensor products contain generators for the dual cone of
 300     # the cross-positive transformations.
 301     tensor_products = [ s.tensor_product(x)
 302                         for (x,s) in K.discrete_complementarity_set() ]
 303
 304     # Turn our matrices into long vectors...
 305     W = VectorSpace(F, n**2)
 306     vectors = [ W(m.list()) for m in tensor_products ]
 307
 308     # Create the *dual* cone of the cross-positive operators,
 309     # expressed as long vectors..
 310     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
 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(F, n)
 319     return [ -M(v.list()) for v in Sigma_cone.rays() ]