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 discrete_complementarity_set(K):
 317     r"""
 318     Compute the discrete complementarity set of this cone.
 319
 320     The complementarity set of this cone is the set of all orthogonal
 321     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 322     dual. The discrete complementarity set restricts `x` and `s` to be
 323     generators of their respective cones.
 324
 325     OUTPUT:
 326
 327     A list of pairs `(x,s)` such that,
 328
 329       * `x` is in this cone.
 330       * `x` is a generator of this cone.
 331       * `s` is in this cone's dual.
 332       * `s` is a generator of this cone's dual.
 333       * `x` and `s` are orthogonal.
 334
 335     EXAMPLES:
 336
 337     The discrete complementarity set of the nonnegative orthant consists
 338     of pairs of standard basis vectors::
 339
 340         sage: K = Cone([(1,0),(0,1)])
 341         sage: discrete_complementarity_set(K)
 342         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 343
 344     If the cone consists of a single ray, the second components of the
 345     discrete complementarity set should generate the orthogonal
 346     complement of that ray::
 347
 348         sage: K = Cone([(1,0)])
 349         sage: discrete_complementarity_set(K)
 350         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 351         sage: K = Cone([(1,0,0)])
 352         sage: discrete_complementarity_set(K)
 353         [((1, 0, 0), (0, 1, 0)),
 354           ((1, 0, 0), (0, -1, 0)),
 355           ((1, 0, 0), (0, 0, 1)),
 356           ((1, 0, 0), (0, 0, -1))]
 357
 358     When the cone is the entire space, its dual is the trivial cone, so
 359     the discrete complementarity set is empty::
 360
 361         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 362         sage: discrete_complementarity_set(K)
 363         []
 364
 365     TESTS:
 366
 367     The complementarity set of the dual can be obtained by switching the
 368     components of the complementarity set of the original cone::
 369
 370         sage: set_random_seed()
 371         sage: K1 = random_cone(max_dim=6)
 372         sage: K2 = K1.dual()
 373         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 374         sage: actual = discrete_complementarity_set(K1)
 375         sage: sorted(actual) == sorted(expected)
 376         True
 377
 378     """
 379     V = K.lattice().vector_space()
 380
 381     # Convert the rays to vectors so that we can compute inner
 382     # products.
 383     xs = [V(x) for x in K.rays()]
 384     ss = [V(s) for s in K.dual().rays()]
 385
 386     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 387
 388
 389 def LL(K):
 390     r"""
 391     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 392     on this cone.
 393
 394     OUTPUT:
 395
 396     A list of matrices forming a basis for the space of all
 397     Lyapunov-like transformations on the given cone.
 398
 399     EXAMPLES:
 400
 401     The trivial cone has no Lyapunov-like transformations::
 402
 403         sage: L = ToricLattice(0)
 404         sage: K = Cone([], lattice=L)
 405         sage: LL(K)
 406         []
 407
 408     The Lyapunov-like transformations on the nonnegative orthant are
 409     simply diagonal matrices::
 410
 411         sage: K = Cone([(1,)])
 412         sage: LL(K)
 413         [[1]]
 414
 415         sage: K = Cone([(1,0),(0,1)])
 416         sage: LL(K)
 417         [
 418         [1 0]  [0 0]
 419         [0 0], [0 1]
 420         ]
 421
 422         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 423         sage: LL(K)
 424         [
 425         [1 0 0]  [0 0 0]  [0 0 0]
 426         [0 0 0]  [0 1 0]  [0 0 0]
 427         [0 0 0], [0 0 0], [0 0 1]
 428         ]
 429
 430     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 431     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 432
 433         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 434         sage: LL(L31)
 435         [
 436         [1 0 0]
 437         [0 1 0]
 438         [0 0 1]
 439         ]
 440
 441         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 442         sage: LL(L3infty)
 443         [
 444         [1 0 0]
 445         [0 1 0]
 446         [0 0 1]
 447         ]
 448
 449     If our cone is the entire space, then every transformation on it is
 450     Lyapunov-like::
 451
 452         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 453         sage: M = MatrixSpace(QQ,2)
 454         sage: M.basis() == LL(K)
 455         True
 456
 457     TESTS:
 458
 459     The inner product `\left< L\left(x\right), s \right>` is zero for
 460     every pair `\left( x,s \right)` in the discrete complementarity set
 461     of the cone::
 462
 463         sage: set_random_seed()
 464         sage: K = random_cone(max_dim=8)
 465         sage: C_of_K = discrete_complementarity_set(K)
 466         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 467         sage: sum(map(abs, l))
 468         0
 469
 470     The Lyapunov-like transformations on a cone and its dual are related
 471     by transposition, but we're not guaranteed to compute transposed
 472     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 473     \right)`
 474
 475         sage: set_random_seed()
 476         sage: K = random_cone(max_dim=8)
 477         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 478         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 479         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 480         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 481         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 482         True
 483
 484     """
 485     V = K.lattice().vector_space()
 486
 487     C_of_K = discrete_complementarity_set(K)
 488
 489     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 490
 491     # Sage doesn't think matrices are vectors, so we have to convert
 492     # our matrices to vectors explicitly before we can figure out how
 493     # many are linearly-indepenedent.
 494     #
 495     # The space W has the same base ring as V, but dimension
 496     # dim(V)^2. So it has the same dimension as the space of linear
 497     # transformations on V. In other words, it's just the right size
 498     # to create an isomorphism between it and our matrices.
 499     W = VectorSpace(V.base_ring(), V.dimension()**2)
 500
 501     # Turn our matrices into long vectors...
 502     vectors = [ W(m.list()) for m in tensor_products ]
 503
 504     # Vector space representation of Lyapunov-like matrices
 505     # (i.e. vec(L) where L is Luapunov-like).
 506     LL_vector = W.span(vectors).complement()
 507
 508     # Now construct an ambient MatrixSpace in which to stick our
 509     # transformations.
 510     M = MatrixSpace(V.base_ring(), V.dimension())
 511
 512     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 513
 514     return matrix_basis
 515
 516
 517
 518 def lyapunov_rank(K):
 519     r"""
 520     Compute the Lyapunov (or bilinearity) rank of this cone.
 521
 522     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 523
 524     1. The dimension of the Lie algebra of the automorphism group of the
 525        cone.
 526
 527     2. The dimension of the linear space of all Lyapunov-like
 528        transformations on the cone.
 529
 530     INPUT:
 531
 532     A closed, convex polyhedral cone.
 533
 534     OUTPUT:
 535
 536     An integer representing the Lyapunov rank of the cone. If the
 537     dimension of the ambient vector space is `n`, then the Lyapunov rank
 538     will be between `1` and `n` inclusive; however a rank of `n-1` is
 539     not possible (see the first reference).
 540
 541     .. note::
 542
 543         In the references, the cones are always assumed to be proper. We
 544         do not impose this restriction.
 545
 546     .. seealso::
 547
 548         :meth:`is_proper`
 549
 550     ALGORITHM:
 551
 552     The codimension formula from the second reference is used. We find
 553     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 554     and `s` are rays of our cone. It is known that these vectors are
 555     sufficient to apply the codimension formula. Once we have all such
 556     pairs, we "brute force" the codimension formula by finding all
 557     linearly-independent `xs^{T}`.
 558
 559     REFERENCES:
 560
 561     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 562        cone and Lyapunov-like transformations, Mathematical Programming, 147
 563        (2014) 155-170.
 564
 565     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 566        Improper Cone. Work in-progress.
 567
 568     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 569        optimality constraints for the cone of positive polynomials,
 570        Mathematical Programming, Series B, 129 (2011) 5-31.
 571
 572     EXAMPLES:
 573
 574     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 575     [Rudolf et al.]_::
 576
 577         sage: positives = Cone([(1,)])
 578         sage: lyapunov_rank(positives)
 579         1
 580         sage: quadrant = Cone([(1,0), (0,1)])
 581         sage: lyapunov_rank(quadrant)
 582         2
 583        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 584         sage: lyapunov_rank(octant)
 585         3
 586
 587     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 588     [Orlitzky/Gowda]_::
 589
 590         sage: R5 = VectorSpace(QQ, 5)
 591         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 592         sage: K = Cone(gs)
 593         sage: lyapunov_rank(K)
 594         25
 595
 596     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 597     [Rudolf et al.]_::
 598
 599         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 600         sage: lyapunov_rank(L31)
 601         1
 602
 603     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 604
 605         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 606         sage: lyapunov_rank(L3infty)
 607         1
 608
 609     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 610     + 1` [Orlitzky/Gowda]_::
 611
 612         sage: K = Cone([(1,0,0,0,0)])
 613         sage: lyapunov_rank(K)
 614         21
 615         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 616         21
 617
 618     A subspace (of dimension `m`) in `n` dimensions should have a
 619     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 620
 621         sage: e1 = (1,0,0,0,0)
 622         sage: neg_e1 = (-1,0,0,0,0)
 623         sage: e2 = (0,1,0,0,0)
 624         sage: neg_e2 = (0,-1,0,0,0)
 625         sage: z = (0,0,0,0,0)
 626         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 627         sage: lyapunov_rank(K)
 628         19
 629         sage: K.lattice_dim()**2 - K.dim()*codim(K)
 630         19
 631
 632     The Lyapunov rank should be additive on a product of proper cones
 633     [Rudolf et al.]_::
 634
 635         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 636         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 637         sage: K = L31.cartesian_product(octant)
 638         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 639         True
 640
 641     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 642     The cone ``K`` in the following example is isomorphic to the nonnegative
 643     octant in `\mathbb{R}^{3}`::
 644
 645         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 646         sage: lyapunov_rank(K)
 647         3
 648
 649     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 650     itself [Rudolf et al.]_::
 651
 652         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 653         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 654         True
 655
 656     TESTS:
 657
 658     The Lyapunov rank should be additive on a product of proper cones
 659     [Rudolf et al.]_::
 660
 661         sage: set_random_seed()
 662         sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 663         sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 664         sage: K = K1.cartesian_product(K2)
 665         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 666         True
 667
 668     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 669     itself [Rudolf et al.]_::
 670
 671         sage: set_random_seed()
 672         sage: K = random_cone(max_dim=8)
 673         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 674         True
 675
 676     Make sure we exercise the non-strictly-convex/non-solid case::
 677
 678         sage: set_random_seed()
 679         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 680         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 681         True
 682
 683     Let's check the other permutations as well, just to be sure::
 684
 685         sage: set_random_seed()
 686         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 687         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 688         True
 689
 690     ::
 691
 692         sage: set_random_seed()
 693         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 694         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 695         True
 696
 697     ::
 698
 699         sage: set_random_seed()
 700         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 701         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 702         True
 703
 704     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 705     be any number between `1` and `n` inclusive, excluding `n-1`
 706     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 707     trivial cone in a trivial space as well. However, in zero dimensions,
 708     the Lyapunov rank of the trivial cone will be zero::
 709
 710         sage: set_random_seed()
 711         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 712         sage: b = lyapunov_rank(K)
 713         sage: n = K.lattice_dim()
 714         sage: (n == 0 or 1 <= b) and b <= n
 715         True
 716         sage: b == n-1
 717         False
 718
 719     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 720     Lyapunov rank `n-1` in `n` dimensions::
 721
 722         sage: set_random_seed()
 723         sage: K = random_cone(max_dim=8)
 724         sage: b = lyapunov_rank(K)
 725         sage: n = K.lattice_dim()
 726         sage: b == n-1
 727         False
 728
 729     The calculation of the Lyapunov rank of an improper cone can be
 730     reduced to that of a proper cone [Orlitzky/Gowda]_::
 731
 732         sage: set_random_seed()
 733         sage: K = random_cone(max_dim=8)
 734         sage: actual = lyapunov_rank(K)
 735         sage: K_S = rho(K)
 736         sage: P = rho(K_S.dual()).dual()
 737         sage: l = lineality(K)
 738         sage: c = codim(K)
 739         sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
 740         sage: actual == expected
 741         True
 742
 743     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 744
 745         sage: set_random_seed()
 746         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 747         sage: lyapunov_rank(K) == len(LL(K))
 748         True
 749
 750     In fact the same can be said of any cone. These additional tests
 751     just increase our confidence that the reduction scheme works::
 752
 753         sage: set_random_seed()
 754         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 755         sage: lyapunov_rank(K) == len(LL(K))
 756         True
 757
 758     ::
 759
 760         sage: set_random_seed()
 761         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 762         sage: lyapunov_rank(K) == len(LL(K))
 763         True
 764
 765     ::
 766
 767         sage: set_random_seed()
 768         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 769         sage: lyapunov_rank(K) == len(LL(K))
 770         True
 771
 772     """
 773     K_orig = K
 774     beta = 0
 775
 776     m = K.dim()
 777     n = K.lattice_dim()
 778     l = lineality(K)
 779
 780     if m < n:
 781         # K is not solid, project onto its span.
 782         K = rho(K)
 783
 784         # Lemma 2
 785         beta += m*(n - m) + (n - m)**2
 786
 787     if l > 0:
 788         # K is not pointed, project its dual onto its span.
 789         # Uses a proposition from our paper, i.e. this is
 790         # equivalent to K = rho(K.dual()).dual()
 791         K = rho(K, K.dual())
 792
 793         # Lemma 3
 794         beta += m * l
 795
 796     beta += len(LL(K))
 797     return beta