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 basically_the_same(K1,K2):
  12     r"""
  13     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  14     otherwise. This is intended as a lazy way to check whether or not
  15     ``K1`` and ``K2`` are linearly isomorphic (i.e. ``A(K1) == K2`` for
  16     some invertible linear transformation ``A``).
  17     """
  18     if K1.lattice_dim() != K2.lattice_dim():
  19         return False
  20
  21     if K1.nrays() != K2.nrays():
  22         return False
  23
  24     if K1.dim() != K2.dim():
  25         return False
  26
  27     if lineality(K1) != lineality(K2):
  28         return False
  29
  30     if K1.is_solid() != K2.is_solid():
  31         return False
  32
  33     if K1.is_strictly_convex() != K2.is_strictly_convex():
  34         return False
  35
  36     if len(LL(K1)) != len(LL(K2)):
  37         return False
  38
  39     C_of_K1 = discrete_complementarity_set(K1)
  40     C_of_K2 = discrete_complementarity_set(K2)
  41     if len(C_of_K1) != len(C_of_K2):
  42         return False
  43
  44     if len(K1.facets()) != len(K2.facets()):
  45         return False
  46
  47     return True
  48
  49
  50
  51 def rho(K, K2=None):
  52     r"""
  53     Restrict ``K`` into its own span, or the span of another cone.
  54
  55     INPUT:
  56
  57     - ``K2`` -- another cone whose lattice has the same rank as this
  58                 cone.
  59
  60     OUTPUT:
  61
  62     A new cone in a sublattice.
  63
  64     EXAMPLES::
  65
  66         sage: K = Cone([(1,)])
  67         sage: rho(K) == K
  68         True
  69
  70         sage: K2 = Cone([(1,0)])
  71         sage: rho(K2).rays()
  72         N(1)
  73         in 1-d lattice N
  74         sage: K3 = Cone([(1,0,0)])
  75         sage: rho(K3).rays()
  76         N(1)
  77         in 1-d lattice N
  78         sage: rho(K2) == rho(K3)
  79         True
  80
  81     TESTS:
  82
  83     The projected cone should always be solid::
  84
  85         sage: set_random_seed()
  86         sage: K = random_cone(max_dim = 8)
  87         sage: K_S = rho(K)
  88         sage: K_S.is_solid()
  89         True
  90
  91     And the resulting cone should live in a space having the same
  92     dimension as the space we restricted it to::
  93
  94         sage: set_random_seed()
  95         sage: K = random_cone(max_dim = 8)
  96         sage: K_S = rho(K, K.dual() )
  97         sage: K_S.lattice_dim() == K.dual().dim()
  98         True
  99
 100     This function should not affect the dimension of a cone::
 101
 102         sage: set_random_seed()
 103         sage: K = random_cone(max_dim = 8)
 104         sage: K.dim() == rho(K).dim()
 105         True
 106
 107     Nor should it affect the lineality of a cone::
 108
 109         sage: set_random_seed()
 110         sage: K = random_cone(max_dim = 8)
 111         sage: lineality(K) == lineality(rho(K))
 112         True
 113
 114     No matter which space we restrict to, the lineality should not
 115     increase::
 116
 117         sage: set_random_seed()
 118         sage: K = random_cone(max_dim = 8)
 119         sage: lineality(K) >= lineality(rho(K))
 120         True
 121         sage: lineality(K) >= lineality(rho(K, K.dual()))
 122         True
 123
 124     If we do this according to our paper, then the result is proper::
 125
 126         sage: set_random_seed()
 127         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
 128         sage: K_S = rho(K)
 129         sage: K_SP = rho(K_S.dual()).dual()
 130         sage: K_SP.is_proper()
 131         True
 132         sage: K_SP = rho(K_S, K_S.dual())
 133         sage: K_SP.is_proper()
 134         True
 135
 136     ::
 137
 138         sage: set_random_seed()
 139         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
 140         sage: K_S = rho(K)
 141         sage: K_SP = rho(K_S.dual()).dual()
 142         sage: K_SP.is_proper()
 143         True
 144         sage: K_SP = rho(K_S, K_S.dual())
 145         sage: K_SP.is_proper()
 146         True
 147
 148     ::
 149
 150         sage: set_random_seed()
 151         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
 152         sage: K_S = rho(K)
 153         sage: K_SP = rho(K_S.dual()).dual()
 154         sage: K_SP.is_proper()
 155         True
 156         sage: K_SP = rho(K_S, K_S.dual())
 157         sage: K_SP.is_proper()
 158         True
 159
 160     ::
 161
 162         sage: set_random_seed()
 163         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
 164         sage: K_S = rho(K)
 165         sage: K_SP = rho(K_S.dual()).dual()
 166         sage: K_SP.is_proper()
 167         True
 168         sage: K_SP = rho(K_S, K_S.dual())
 169         sage: K_SP.is_proper()
 170         True
 171
 172     Test the proposition in our paper concerning the duals and
 173     restrictions. Generate a random cone, then create a subcone of
 174     it. The operation of dual-taking should then commute with rho::
 175
 176         sage: set_random_seed()
 177         sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=False)
 178         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 179         sage: K_W = rho(K, J)
 180         sage: K_star_W_star = rho(K.dual(), J).dual()
 181         sage: basically_the_same(K_W, K_star_W_star)
 182         True
 183
 184     ::
 185
 186         sage: set_random_seed()
 187         sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=False)
 188         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 189         sage: K_W = rho(K, J)
 190         sage: K_star_W_star = rho(K.dual(), J).dual()
 191         sage: basically_the_same(K_W, K_star_W_star)
 192         True
 193
 194     ::
 195
 196         sage: set_random_seed()
 197         sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=True)
 198         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 199         sage: K_W = rho(K, J)
 200         sage: K_star_W_star = rho(K.dual(), J).dual()
 201         sage: basically_the_same(K_W, K_star_W_star)
 202         True
 203
 204     ::
 205
 206         sage: set_random_seed()
 207         sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=True)
 208         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 209         sage: K_W = rho(K, J)
 210         sage: K_star_W_star = rho(K.dual(), J).dual()
 211         sage: basically_the_same(K_W, K_star_W_star)
 212         True
 213
 214     """
 215     if K2 is None:
 216         K2 = K
 217
 218     # First we project K onto the span of K2. This will explode if the
 219     # rank of ``K2.lattice()`` doesn't match ours.
 220     span_K2 = Cone(K2.rays() + (-K2).rays(), lattice=K.lattice())
 221     K = K.intersection(span_K2)
 222
 223     # Cheat a little to get the subspace span(K2). The paper uses the
 224     # rays of K2 as a basis, but everything is invariant under linear
 225     # isomorphism (i.e. a change of basis), and this is a little
 226     # faster.
 227     W = span_K2.linear_subspace()
 228
 229     # We've already intersected K with the span of K2, so every
 230     # generator of K should belong to W now.
 231     W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 232
 233     L = ToricLattice(K2.dim())
 234     return Cone(W_rays, lattice=L)
 235
 236
 237
 238 def lineality(K):
 239     r"""
 240     Compute the lineality of this cone.
 241
 242     The lineality of a cone is the dimension of the largest linear
 243     subspace contained in that cone.
 244
 245     OUTPUT:
 246
 247     A nonnegative integer; the dimension of the largest subspace
 248     contained within this cone.
 249
 250     REFERENCES:
 251
 252     .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
 253        University Press, Princeton, 1970.
 254
 255     EXAMPLES:
 256
 257     The lineality of the nonnegative orthant is zero, since it clearly
 258     contains no lines::
 259
 260         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
 261         sage: lineality(K)
 262         0
 263
 264     However, if we add another ray so that the entire `x`-axis belongs
 265     to the cone, then the resulting cone will have lineality one::
 266
 267         sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
 268         sage: lineality(K)
 269         1
 270
 271     If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
 272     to the dimension of the ambient space (i.e. two)::
 273
 274         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 275         sage: lineality(K)
 276         2
 277
 278     Per the definition, the lineality of the trivial cone in a trivial
 279     space is zero::
 280
 281         sage: K = Cone([], lattice=ToricLattice(0))
 282         sage: lineality(K)
 283         0
 284
 285     TESTS:
 286
 287     The lineality of a cone should be an integer between zero and the
 288     dimension of the ambient space, inclusive::
 289
 290         sage: set_random_seed()
 291         sage: K = random_cone(max_dim = 8)
 292         sage: l = lineality(K)
 293         sage: l in ZZ
 294         True
 295         sage: (0 <= l) and (l <= K.lattice_dim())
 296         True
 297
 298     A strictly convex cone should have lineality zero::
 299
 300         sage: set_random_seed()
 301         sage: K = random_cone(max_dim = 8, strictly_convex = True)
 302         sage: lineality(K)
 303         0
 304
 305     """
 306     return K.linear_subspace().dimension()
 307
 308
 309 def codim(K):
 310     r"""
 311     Compute the codimension of this cone.
 312
 313     The codimension of a cone is the dimension of the space of all
 314     elements perpendicular to every element of the cone. In other words,
 315     the codimension is the difference between the dimension of the
 316     ambient space and the dimension of the cone itself.
 317
 318     OUTPUT:
 319
 320     A nonnegative integer representing the dimension of the space of all
 321     elements perpendicular to this cone.
 322
 323     .. seealso::
 324
 325         :meth:`dim`, :meth:`lattice_dim`
 326
 327     EXAMPLES:
 328
 329     The codimension of the nonnegative orthant is zero, since the span of
 330     its generators equals the entire ambient space::
 331
 332         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
 333         sage: codim(K)
 334         0
 335
 336     However, if we remove a ray so that the entire cone is contained
 337     within the `x-y`-plane, then the resulting cone will have
 338     codimension one, because the `z`-axis is perpendicular to every
 339     element of the cone::
 340
 341         sage: K = Cone([(1,0,0), (0,1,0)])
 342         sage: codim(K)
 343         1
 344
 345     If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
 346
 347         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 348         sage: codim(K)
 349         0
 350
 351     And if the cone is trivial in any space, then its codimension is
 352     equal to the dimension of the ambient space::
 353
 354         sage: K = Cone([], lattice=ToricLattice(0))
 355         sage: K.lattice_dim()
 356         0
 357         sage: codim(K)
 358         0
 359
 360         sage: K = Cone([(0,)])
 361         sage: K.lattice_dim()
 362         1
 363         sage: codim(K)
 364         1
 365
 366         sage: K = Cone([(0,0)])
 367         sage: K.lattice_dim()
 368         2
 369         sage: codim(K)
 370         2
 371
 372     TESTS:
 373
 374     The codimension of a cone should be an integer between zero and
 375     the dimension of the ambient space, inclusive::
 376
 377         sage: set_random_seed()
 378         sage: K = random_cone(max_dim = 8)
 379         sage: c = codim(K)
 380         sage: c in ZZ
 381         True
 382         sage: (0 <= c) and (c <= K.lattice_dim())
 383         True
 384
 385     A solid cone should have codimension zero::
 386
 387         sage: set_random_seed()
 388         sage: K = random_cone(max_dim = 8, solid = True)
 389         sage: codim(K)
 390         0
 391
 392     The codimension of a cone is equal to the lineality of its dual::
 393
 394         sage: set_random_seed()
 395         sage: K = random_cone(max_dim = 8, solid = True)
 396         sage: codim(K) == lineality(K.dual())
 397         True
 398
 399     """
 400     return (K.lattice_dim() - K.dim())
 401
 402
 403 def discrete_complementarity_set(K):
 404     r"""
 405     Compute the discrete complementarity set of this cone.
 406
 407     The complementarity set of this cone is the set of all orthogonal
 408     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 409     dual. The discrete complementarity set restricts `x` and `s` to be
 410     generators of their respective cones.
 411
 412     OUTPUT:
 413
 414     A list of pairs `(x,s)` such that,
 415
 416       * `x` is a generator of this cone.
 417       * `s` is a generator of this cone's dual.
 418       * `x` and `s` are orthogonal.
 419
 420     EXAMPLES:
 421
 422     The discrete complementarity set of the nonnegative orthant consists
 423     of pairs of standard basis vectors::
 424
 425         sage: K = Cone([(1,0),(0,1)])
 426         sage: discrete_complementarity_set(K)
 427         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 428
 429     If the cone consists of a single ray, the second components of the
 430     discrete complementarity set should generate the orthogonal
 431     complement of that ray::
 432
 433         sage: K = Cone([(1,0)])
 434         sage: discrete_complementarity_set(K)
 435         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 436         sage: K = Cone([(1,0,0)])
 437         sage: discrete_complementarity_set(K)
 438         [((1, 0, 0), (0, 1, 0)),
 439           ((1, 0, 0), (0, -1, 0)),
 440           ((1, 0, 0), (0, 0, 1)),
 441           ((1, 0, 0), (0, 0, -1))]
 442
 443     When the cone is the entire space, its dual is the trivial cone, so
 444     the discrete complementarity set is empty::
 445
 446         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 447         sage: discrete_complementarity_set(K)
 448         []
 449
 450     TESTS:
 451
 452     The complementarity set of the dual can be obtained by switching the
 453     components of the complementarity set of the original cone::
 454
 455         sage: set_random_seed()
 456         sage: K1 = random_cone(max_dim=6)
 457         sage: K2 = K1.dual()
 458         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 459         sage: actual = discrete_complementarity_set(K1)
 460         sage: sorted(actual) == sorted(expected)
 461         True
 462
 463     """
 464     V = K.lattice().vector_space()
 465
 466     # Convert the rays to vectors so that we can compute inner
 467     # products.
 468     xs = [V(x) for x in K.rays()]
 469     ss = [V(s) for s in K.dual().rays()]
 470
 471     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 472
 473
 474 def LL(K):
 475     r"""
 476     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 477     on this cone.
 478
 479     OUTPUT:
 480
 481     A list of matrices forming a basis for the space of all
 482     Lyapunov-like transformations on the given cone.
 483
 484     EXAMPLES:
 485
 486     The trivial cone has no Lyapunov-like transformations::
 487
 488         sage: L = ToricLattice(0)
 489         sage: K = Cone([], lattice=L)
 490         sage: LL(K)
 491         []
 492
 493     The Lyapunov-like transformations on the nonnegative orthant are
 494     simply diagonal matrices::
 495
 496         sage: K = Cone([(1,)])
 497         sage: LL(K)
 498         [[1]]
 499
 500         sage: K = Cone([(1,0),(0,1)])
 501         sage: LL(K)
 502         [
 503         [1 0]  [0 0]
 504         [0 0], [0 1]
 505         ]
 506
 507         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 508         sage: LL(K)
 509         [
 510         [1 0 0]  [0 0 0]  [0 0 0]
 511         [0 0 0]  [0 1 0]  [0 0 0]
 512         [0 0 0], [0 0 0], [0 0 1]
 513         ]
 514
 515     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 516     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 517
 518         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 519         sage: LL(L31)
 520         [
 521         [1 0 0]
 522         [0 1 0]
 523         [0 0 1]
 524         ]
 525
 526         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 527         sage: LL(L3infty)
 528         [
 529         [1 0 0]
 530         [0 1 0]
 531         [0 0 1]
 532         ]
 533
 534     If our cone is the entire space, then every transformation on it is
 535     Lyapunov-like::
 536
 537         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 538         sage: M = MatrixSpace(QQ,2)
 539         sage: M.basis() == LL(K)
 540         True
 541
 542     TESTS:
 543
 544     The inner product `\left< L\left(x\right), s \right>` is zero for
 545     every pair `\left( x,s \right)` in the discrete complementarity set
 546     of the cone::
 547
 548         sage: set_random_seed()
 549         sage: K = random_cone(max_dim=8)
 550         sage: C_of_K = discrete_complementarity_set(K)
 551         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 552         sage: sum(map(abs, l))
 553         0
 554
 555     The Lyapunov-like transformations on a cone and its dual are related
 556     by transposition, but we're not guaranteed to compute transposed
 557     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 558     \right)`
 559
 560         sage: set_random_seed()
 561         sage: K = random_cone(max_dim=8)
 562         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 563         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 564         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 565         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 566         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 567         True
 568
 569     """
 570     V = K.lattice().vector_space()
 571
 572     C_of_K = discrete_complementarity_set(K)
 573
 574     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 575
 576     # Sage doesn't think matrices are vectors, so we have to convert
 577     # our matrices to vectors explicitly before we can figure out how
 578     # many are linearly-indepenedent.
 579     #
 580     # The space W has the same base ring as V, but dimension
 581     # dim(V)^2. So it has the same dimension as the space of linear
 582     # transformations on V. In other words, it's just the right size
 583     # to create an isomorphism between it and our matrices.
 584     W = VectorSpace(V.base_ring(), V.dimension()**2)
 585
 586     # Turn our matrices into long vectors...
 587     vectors = [ W(m.list()) for m in tensor_products ]
 588
 589     # Vector space representation of Lyapunov-like matrices
 590     # (i.e. vec(L) where L is Luapunov-like).
 591     LL_vector = W.span(vectors).complement()
 592
 593     # Now construct an ambient MatrixSpace in which to stick our
 594     # transformations.
 595     M = MatrixSpace(V.base_ring(), V.dimension())
 596
 597     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 598
 599     return matrix_basis
 600
 601
 602
 603 def lyapunov_rank(K):
 604     r"""
 605     Compute the Lyapunov (or bilinearity) rank of this cone.
 606
 607     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 608
 609     1. The dimension of the Lie algebra of the automorphism group of the
 610        cone.
 611
 612     2. The dimension of the linear space of all Lyapunov-like
 613        transformations on the cone.
 614
 615     INPUT:
 616
 617     A closed, convex polyhedral cone.
 618
 619     OUTPUT:
 620
 621     An integer representing the Lyapunov rank of the cone. If the
 622     dimension of the ambient vector space is `n`, then the Lyapunov rank
 623     will be between `1` and `n` inclusive; however a rank of `n-1` is
 624     not possible (see the first reference).
 625
 626     .. note::
 627
 628         In the references, the cones are always assumed to be proper. We
 629         do not impose this restriction.
 630
 631     .. seealso::
 632
 633         :meth:`is_proper`
 634
 635     ALGORITHM:
 636
 637     The codimension formula from the second reference is used. We find
 638     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 639     and `s` are rays of our cone. It is known that these vectors are
 640     sufficient to apply the codimension formula. Once we have all such
 641     pairs, we "brute force" the codimension formula by finding all
 642     linearly-independent `xs^{T}`.
 643
 644     REFERENCES:
 645
 646     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 647        cone and Lyapunov-like transformations, Mathematical Programming, 147
 648        (2014) 155-170.
 649
 650     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 651        Improper Cone. Work in-progress.
 652
 653     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 654        optimality constraints for the cone of positive polynomials,
 655        Mathematical Programming, Series B, 129 (2011) 5-31.
 656
 657     EXAMPLES:
 658
 659     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 660     [Rudolf et al.]_::
 661
 662         sage: positives = Cone([(1,)])
 663         sage: lyapunov_rank(positives)
 664         1
 665         sage: quadrant = Cone([(1,0), (0,1)])
 666         sage: lyapunov_rank(quadrant)
 667         2
 668        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 669         sage: lyapunov_rank(octant)
 670         3
 671
 672     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 673     [Orlitzky/Gowda]_::
 674
 675         sage: R5 = VectorSpace(QQ, 5)
 676         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 677         sage: K = Cone(gs)
 678         sage: lyapunov_rank(K)
 679         25
 680
 681     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 682     [Rudolf et al.]_::
 683
 684         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 685         sage: lyapunov_rank(L31)
 686         1
 687
 688     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 689
 690         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 691         sage: lyapunov_rank(L3infty)
 692         1
 693
 694     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 695     + 1` [Orlitzky/Gowda]_::
 696
 697         sage: K = Cone([(1,0,0,0,0)])
 698         sage: lyapunov_rank(K)
 699         21
 700         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 701         21
 702
 703     A subspace (of dimension `m`) in `n` dimensions should have a
 704     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 705
 706         sage: e1 = (1,0,0,0,0)
 707         sage: neg_e1 = (-1,0,0,0,0)
 708         sage: e2 = (0,1,0,0,0)
 709         sage: neg_e2 = (0,-1,0,0,0)
 710         sage: z = (0,0,0,0,0)
 711         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 712         sage: lyapunov_rank(K)
 713         19
 714         sage: K.lattice_dim()**2 - K.dim()*codim(K)
 715         19
 716
 717     The Lyapunov rank should be additive on a product of proper cones
 718     [Rudolf et al.]_::
 719
 720         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 721         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 722         sage: K = L31.cartesian_product(octant)
 723         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 724         True
 725
 726     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 727     The cone ``K`` in the following example is isomorphic to the nonnegative
 728     octant in `\mathbb{R}^{3}`::
 729
 730         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 731         sage: lyapunov_rank(K)
 732         3
 733
 734     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 735     itself [Rudolf et al.]_::
 736
 737         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 738         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 739         True
 740
 741     TESTS:
 742
 743     The Lyapunov rank should be additive on a product of proper cones
 744     [Rudolf et al.]_::
 745
 746         sage: set_random_seed()
 747         sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 748         sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 749         sage: K = K1.cartesian_product(K2)
 750         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 751         True
 752
 753     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 754     itself [Rudolf et al.]_::
 755
 756         sage: set_random_seed()
 757         sage: K = random_cone(max_dim=8)
 758         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 759         True
 760
 761     Make sure we exercise the non-strictly-convex/non-solid case::
 762
 763         sage: set_random_seed()
 764         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 765         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 766         True
 767
 768     Let's check the other permutations as well, just to be sure::
 769
 770         sage: set_random_seed()
 771         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 772         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 773         True
 774
 775     ::
 776
 777         sage: set_random_seed()
 778         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 779         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 780         True
 781
 782     ::
 783
 784         sage: set_random_seed()
 785         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 786         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 787         True
 788
 789     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 790     be any number between `1` and `n` inclusive, excluding `n-1`
 791     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 792     trivial cone in a trivial space as well. However, in zero dimensions,
 793     the Lyapunov rank of the trivial cone will be zero::
 794
 795         sage: set_random_seed()
 796         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 797         sage: b = lyapunov_rank(K)
 798         sage: n = K.lattice_dim()
 799         sage: (n == 0 or 1 <= b) and b <= n
 800         True
 801         sage: b == n-1
 802         False
 803
 804     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 805     Lyapunov rank `n-1` in `n` dimensions::
 806
 807         sage: set_random_seed()
 808         sage: K = random_cone(max_dim=8)
 809         sage: b = lyapunov_rank(K)
 810         sage: n = K.lattice_dim()
 811         sage: b == n-1
 812         False
 813
 814     The calculation of the Lyapunov rank of an improper cone can be
 815     reduced to that of a proper cone [Orlitzky/Gowda]_::
 816
 817         sage: set_random_seed()
 818         sage: K = random_cone(max_dim=8)
 819         sage: actual = lyapunov_rank(K)
 820         sage: K_S = rho(K)
 821         sage: K_SP = rho(K_S.dual()).dual()
 822         sage: l = lineality(K)
 823         sage: c = codim(K)
 824         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 825         sage: actual == expected
 826         True
 827
 828     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 829
 830         sage: set_random_seed()
 831         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 832         sage: lyapunov_rank(K) == len(LL(K))
 833         True
 834
 835     In fact the same can be said of any cone. These additional tests
 836     just increase our confidence that the reduction scheme works::
 837
 838         sage: set_random_seed()
 839         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 840         sage: lyapunov_rank(K) == len(LL(K))
 841         True
 842
 843     ::
 844
 845         sage: set_random_seed()
 846         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 847         sage: lyapunov_rank(K) == len(LL(K))
 848         True
 849
 850     ::
 851
 852         sage: set_random_seed()
 853         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 854         sage: lyapunov_rank(K) == len(LL(K))
 855         True
 856
 857     """
 858     beta = 0
 859
 860     m = K.dim()
 861     n = K.lattice_dim()
 862     l = lineality(K)
 863
 864     if m < n:
 865         # K is not solid, restrict to its span.
 866         K = rho(K)
 867
 868         # Lemma 2
 869         beta += m*(n - m) + (n - m)**2
 870
 871     if l > 0:
 872         # K is not pointed, restrict to the span of its dual. Uses a
 873         # proposition from our paper, i.e. this is equivalent to K =
 874         # rho(K.dual()).dual().
 875         K = rho(K, K.dual())
 876
 877         # Lemma 3
 878         beta += m * l
 879
 880     beta += len(LL(K))
 881     return beta