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