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