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