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