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