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 rename_lattice(L,s):
  11     r"""
  12     Change all names of the given lattice to ``s``.
  13     """
  14     L._name = s
  15     L._dual_name = s
  16     L._latex_name = s
  17     L._latex_dual_name = s
  18
  19 def span_iso(K):
  20     r"""
  21     Return an isomorphism (and its inverse) that will send ``K`` into a
  22     lower-dimensional space isomorphic to its span (and back).
  23
  24     EXAMPLES:
  25
  26     The inverse composed with the isomorphism should be the identity::
  27
  28         sage: K = random_cone(max_dim=10)
  29         sage: (phi, phi_inv) = span_iso(K)
  30         sage: phi_inv(phi(K)) == K
  31         True
  32
  33     The image of ``K`` under the isomorphism should have full dimension::
  34
  35         sage: K = random_cone(max_dim=10)
  36         sage: (phi, phi_inv) = span_iso(K)
  37         sage: phi(K).dim() == phi(K).lattice_dim()
  38         True
  39
  40     """
  41     phi_domain = K.sublattice().vector_space()
  42     phi_codo = VectorSpace(phi_domain.base_field(), phi_domain.dimension())
  43
  44     # S goes from the new space to the cone space.
  45     S = linear_transformation(phi_codo, phi_domain, phi_domain.basis())
  46
  47     # phi goes from the cone space to the new space.
  48     def phi(J_orig):
  49         r"""
  50         Takes a cone ``J`` and sends it into the new space.
  51         """
  52         newrays = map(S.inverse(), J_orig.rays())
  53         L = None
  54         if len(newrays) == 0:
  55             L = ToricLattice(0)
  56
  57         return Cone(newrays, lattice=L)
  58
  59     def phi_inverse(J_sub):
  60         r"""
  61         The inverse to phi which goes from the new space to the cone space.
  62         """
  63         newrays = map(S, J_sub.rays())
  64         return Cone(newrays, lattice=K.lattice())
  65
  66
  67     return (phi, phi_inverse)
  68
  69
  70
  71 def discrete_complementarity_set(K):
  72     r"""
  73     Compute the discrete complementarity set of this cone.
  74
  75     The complementarity set of this cone is the set of all orthogonal
  76     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
  77     dual. The discrete complementarity set restricts `x` and `s` to be
  78     generators of their respective cones.
  79
  80     OUTPUT:
  81
  82     A list of pairs `(x,s)` such that,
  83
  84       * `x` is in this cone.
  85       * `x` is a generator of this cone.
  86       * `s` is in this cone's dual.
  87       * `s` is a generator of this cone's dual.
  88       * `x` and `s` are orthogonal.
  89
  90     EXAMPLES:
  91
  92     The discrete complementarity set of the nonnegative orthant consists
  93     of pairs of standard basis vectors::
  94
  95         sage: K = Cone([(1,0),(0,1)])
  96         sage: discrete_complementarity_set(K)
  97         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
  98
  99     If the cone consists of a single ray, the second components of the
 100     discrete complementarity set should generate the orthogonal
 101     complement of that ray::
 102
 103         sage: K = Cone([(1,0)])
 104         sage: discrete_complementarity_set(K)
 105         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 106         sage: K = Cone([(1,0,0)])
 107         sage: discrete_complementarity_set(K)
 108         [((1, 0, 0), (0, 1, 0)),
 109           ((1, 0, 0), (0, -1, 0)),
 110           ((1, 0, 0), (0, 0, 1)),
 111           ((1, 0, 0), (0, 0, -1))]
 112
 113     When the cone is the entire space, its dual is the trivial cone, so
 114     the discrete complementarity set is empty::
 115
 116         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 117         sage: discrete_complementarity_set(K)
 118         []
 119
 120     TESTS:
 121
 122     The complementarity set of the dual can be obtained by switching the
 123     components of the complementarity set of the original cone::
 124
 125         sage: K1 = random_cone(max_dim=10, max_rays=10)
 126         sage: K2 = K1.dual()
 127         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 128         sage: actual = discrete_complementarity_set(K1)
 129         sage: actual == expected
 130         True
 131
 132     """
 133     V = K.lattice().vector_space()
 134
 135     # Convert the rays to vectors so that we can compute inner
 136     # products.
 137     xs = [V(x) for x in K.rays()]
 138     ss = [V(s) for s in K.dual().rays()]
 139
 140     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 141
 142
 143 def LL(K):
 144     r"""
 145     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 146     on this cone.
 147
 148     OUTPUT:
 149
 150     A list of matrices forming a basis for the space of all
 151     Lyapunov-like transformations on the given cone.
 152
 153     EXAMPLES:
 154
 155     The trivial cone has no Lyapunov-like transformations::
 156
 157         sage: L = ToricLattice(0)
 158         sage: K = Cone([], lattice=L)
 159         sage: LL(K)
 160         []
 161
 162     The Lyapunov-like transformations on the nonnegative orthant are
 163     simply diagonal matrices::
 164
 165         sage: K = Cone([(1,)])
 166         sage: LL(K)
 167         [[1]]
 168
 169         sage: K = Cone([(1,0),(0,1)])
 170         sage: LL(K)
 171         [
 172         [1 0]  [0 0]
 173         [0 0], [0 1]
 174         ]
 175
 176         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 177         sage: LL(K)
 178         [
 179         [1 0 0]  [0 0 0]  [0 0 0]
 180         [0 0 0]  [0 1 0]  [0 0 0]
 181         [0 0 0], [0 0 0], [0 0 1]
 182         ]
 183
 184     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 185     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 186
 187         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 188         sage: LL(L31)
 189         [
 190         [1 0 0]
 191         [0 1 0]
 192         [0 0 1]
 193         ]
 194
 195         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 196         sage: LL(L3infty)
 197         [
 198         [1 0 0]
 199         [0 1 0]
 200         [0 0 1]
 201         ]
 202
 203     TESTS:
 204
 205     The inner product `\left< L\left(x\right), s \right>` is zero for
 206     every pair `\left( x,s \right)` in the discrete complementarity set
 207     of the cone::
 208
 209         sage: K = random_cone(max_dim=8, max_rays=10)
 210         sage: C_of_K = discrete_complementarity_set(K)
 211         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 212         sage: sum(map(abs, l))
 213         0
 214
 215     """
 216     V = K.lattice().vector_space()
 217
 218     C_of_K = discrete_complementarity_set(K)
 219
 220     tensor_products = [s.tensor_product(x) for (x,s) in C_of_K]
 221
 222     # Sage doesn't think matrices are vectors, so we have to convert
 223     # our matrices to vectors explicitly before we can figure out how
 224     # many are linearly-indepenedent.
 225     #
 226     # The space W has the same base ring as V, but dimension
 227     # dim(V)^2. So it has the same dimension as the space of linear
 228     # transformations on V. In other words, it's just the right size
 229     # to create an isomorphism between it and our matrices.
 230     W = VectorSpace(V.base_ring(), V.dimension()**2)
 231
 232     # Turn our matrices into long vectors...
 233     vectors = [ W(m.list()) for m in tensor_products ]
 234
 235     # Vector space representation of Lyapunov-like matrices
 236     # (i.e. vec(L) where L is Luapunov-like).
 237     LL_vector = W.span(vectors).complement()
 238
 239     # Now construct an ambient MatrixSpace in which to stick our
 240     # transformations.
 241     M = MatrixSpace(V.base_ring(), V.dimension())
 242
 243     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 244
 245     return matrix_basis
 246
 247
 248
 249 def lyapunov_rank(K):
 250     r"""
 251     Compute the Lyapunov (or bilinearity) rank of this cone.
 252
 253     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 254
 255     1. The dimension of the Lie algebra of the automorphism group of the
 256        cone.
 257
 258     2. The dimension of the linear space of all Lyapunov-like
 259        transformations on the cone.
 260
 261     INPUT:
 262
 263     A closed, convex polyhedral cone.
 264
 265     OUTPUT:
 266
 267     An integer representing the Lyapunov rank of the cone. If the
 268     dimension of the ambient vector space is `n`, then the Lyapunov rank
 269     will be between `1` and `n` inclusive; however a rank of `n-1` is
 270     not possible (see the first reference).
 271
 272     .. note::
 273
 274         In the references, the cones are always assumed to be proper. We
 275         do not impose this restriction.
 276
 277     .. seealso::
 278
 279         :meth:`is_proper`
 280
 281     ALGORITHM:
 282
 283     The codimension formula from the second reference is used. We find
 284     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 285     and `s` are rays of our cone. It is known that these vectors are
 286     sufficient to apply the codimension formula. Once we have all such
 287     pairs, we "brute force" the codimension formula by finding all
 288     linearly-independent `xs^{T}`.
 289
 290     REFERENCES:
 291
 292     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 293        cone and Lyapunov-like transformations, Mathematical Programming, 147
 294        (2014) 155-170.
 295
 296     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 297        Improper Cone. Work in-progress.
 298
 299     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 300        optimality constraints for the cone of positive polynomials,
 301        Mathematical Programming, Series B, 129 (2011) 5-31.
 302
 303     EXAMPLES:
 304
 305     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 306     [Rudolf et al.]_::
 307
 308         sage: positives = Cone([(1,)])
 309         sage: lyapunov_rank(positives)
 310         1
 311         sage: quadrant = Cone([(1,0), (0,1)])
 312         sage: lyapunov_rank(quadrant)
 313         2
 314        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 315         sage: lyapunov_rank(octant)
 316         3
 317
 318     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 319     [Rudolf et al.]_::
 320
 321         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 322         sage: lyapunov_rank(L31)
 323         1
 324
 325     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 326
 327         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 328         sage: lyapunov_rank(L3infty)
 329         1
 330
 331     The Lyapunov rank should be additive on a product of cones
 332     [Rudolf et al.]_::
 333
 334         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 335         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 336         sage: K = L31.cartesian_product(octant)
 337         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 338         True
 339
 340     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 341     The cone ``K`` in the following example is isomorphic to the nonnegative
 342     octant in `\mathbb{R}^{3}`::
 343
 344         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 345         sage: lyapunov_rank(K)
 346         3
 347
 348     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 349     itself [Rudolf et al.]_::
 350
 351         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 352         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 353         True
 354
 355     TESTS:
 356
 357     The Lyapunov rank should be additive on a product of cones
 358     [Rudolf et al.]_::
 359
 360         sage: K1 = random_cone(max_dim=10, max_rays=10)
 361         sage: K2 = random_cone(max_dim=10, max_rays=10)
 362         sage: K = K1.cartesian_product(K2)
 363         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 364         True
 365
 366     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 367     itself [Rudolf et al.]_::
 368
 369         sage: K = random_cone(max_dim=10, max_rays=10)
 370         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 371         True
 372
 373     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 374     be any number between `1` and `n` inclusive, excluding `n-1`
 375     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 376     trivial cone in a trivial space as well. However, in zero dimensions,
 377     the Lyapunov rank of the trivial cone will be zero::
 378
 379         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 380         sage: b = lyapunov_rank(K)
 381         sage: n = K.lattice_dim()
 382         sage: (n == 0 or 1 <= b) and b <= n
 383         True
 384         sage: b == n-1
 385         False
 386
 387     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 388     Lyapunov rank `n-1` in `n` dimensions::
 389
 390         sage: K = random_cone(max_dim=10)
 391         sage: b = lyapunov_rank(K)
 392         sage: n = K.lattice_dim()
 393         sage: b == n-1
 394         False
 395
 396     The calculation of the Lyapunov rank of an improper cone can be
 397     reduced to that of a proper cone [Orlitzky/Gowda]_::
 398
 399         sage: K = random_cone(max_dim=15, solid=False, strictly_convex=False)
 400         sage: actual = lyapunov_rank(K)
 401         sage: (phi1, _) = span_iso(K)
 402         sage: K_S = phi1(K)
 403         sage: (phi2, _) = span_iso(K_S.dual())
 404         sage: J_T = phi2(K_S.dual()).dual()
 405         sage: l = K.linear_subspace().dimension()
 406         sage: codim = K.lattice_dim() - K.dim()
 407         sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
 408         sage: actual == expected
 409         True
 410
 411     """
 412     return len(LL(K))