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     Test whether or not ``K1`` and ``K2`` are "basically the same."
  14
  15     This is a hack to get around the fact that it's difficult to tell
  16     when two cones are linearly isomorphic. We have a proposition that
  17     equates two cones, but represented over `\mathbb{Q}`, they are
  18     merely linearly isomorphic (not equal). So rather than test for
  19     equality, we test a list of properties that should be preserved
  20     under an invertible linear transformation.
  21
  22     OUTPUT:
  23
  24     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  25     otherwise.
  26
  27     EXAMPLES:
  28
  29     Any proper cone with three generators in `\mathbb{R}^{3}` is
  30     basically the same as the nonnegative orthant::
  31
  32         sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
  33         sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
  34         sage: _basically_the_same(K1, K2)
  35         True
  36
  37     Negating a cone gives you another cone that is basically the same::
  38
  39         sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
  40         sage: _basically_the_same(K, -K)
  41         True
  42
  43     TESTS:
  44
  45     Any cone is basically the same as itself::
  46
  47         sage: K = random_cone(max_ambient_dim = 8)
  48         sage: _basically_the_same(K, K)
  49         True
  50
  51     After applying an invertible matrix to the rows of a cone, the
  52     result should be basically the same as the cone we started with::
  53
  54         sage: K1 = random_cone(max_ambient_dim = 8)
  55         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
  56         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
  57         sage: _basically_the_same(K1, K2)
  58         True
  59
  60     """
  61     if K1.lattice_dim() != K2.lattice_dim():
  62         return False
  63
  64     if K1.nrays() != K2.nrays():
  65         return False
  66
  67     if K1.dim() != K2.dim():
  68         return False
  69
  70     if K1.lineality() != K2.lineality():
  71         return False
  72
  73     if K1.is_solid() != K2.is_solid():
  74         return False
  75
  76     if K1.is_strictly_convex() != K2.is_strictly_convex():
  77         return False
  78
  79     if len(LL(K1)) != len(LL(K2)):
  80         return False
  81
  82     C_of_K1 = discrete_complementarity_set(K1)
  83     C_of_K2 = discrete_complementarity_set(K2)
  84     if len(C_of_K1) != len(C_of_K2):
  85         return False
  86
  87     if len(K1.facets()) != len(K2.facets()):
  88         return False
  89
  90     return True
  91
  92
  93
  94 def _restrict_to_space(K, W):
  95     r"""
  96     Restrict this cone a subspace of its ambient space.
  97
  98     INPUT:
  99
 100     - ``W`` -- The subspace into which this cone will be restricted.
 101
 102     OUTPUT:
 103
 104     A new cone in a sublattice corresponding to ``W``.
 105
 106     EXAMPLES:
 107
 108     When this cone is solid, restricting it into its own span should do
 109     nothing::
 110
 111         sage: K = Cone([(1,)])
 112         sage: _restrict_to_space(K, K.span()) == K
 113         True
 114
 115     A single ray restricted into its own span gives the same output
 116     regardless of the ambient space::
 117
 118         sage: K2 = Cone([(1,0)])
 119         sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
 120         sage: K2_S
 121         N(1)
 122         in 1-d lattice N
 123         sage: K3 = Cone([(1,0,0)])
 124         sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
 125         sage: K3_S
 126         N(1)
 127         in 1-d lattice N
 128         sage: K2_S == K3_S
 129         True
 130
 131     TESTS:
 132
 133     The projected cone should always be solid::
 134
 135         sage: set_random_seed()
 136         sage: K = random_cone(max_ambient_dim = 8)
 137         sage: _restrict_to_space(K, K.span()).is_solid()
 138         True
 139
 140     And the resulting cone should live in a space having the same
 141     dimension as the space we restricted it to::
 142
 143         sage: set_random_seed()
 144         sage: K = random_cone(max_ambient_dim = 8)
 145         sage: K_P = _restrict_to_space(K, K.dual().span())
 146         sage: K_P.lattice_dim() == K.dual().dim()
 147         True
 148
 149     This function should not affect the dimension of a cone::
 150
 151         sage: set_random_seed()
 152         sage: K = random_cone(max_ambient_dim = 8)
 153         sage: K.dim() == _restrict_to_space(K,K.span()).dim()
 154         True
 155
 156     Nor should it affect the lineality of a cone::
 157
 158         sage: set_random_seed()
 159         sage: K = random_cone(max_ambient_dim = 8)
 160         sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
 161         True
 162
 163     No matter which space we restrict to, the lineality should not
 164     increase::
 165
 166         sage: set_random_seed()
 167         sage: K = random_cone(max_ambient_dim = 8)
 168         sage: S = K.span(); P = K.dual().span()
 169         sage: K.lineality() >= _restrict_to_space(K,S).lineality()
 170         True
 171         sage: K.lineality() >= _restrict_to_space(K,P).lineality()
 172         True
 173
 174     If we do this according to our paper, then the result is proper::
 175
 176         sage: set_random_seed()
 177         sage: K = random_cone(max_ambient_dim = 8)
 178         sage: K_S = _restrict_to_space(K, K.span())
 179         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 180         sage: K_SP.is_proper()
 181         True
 182         sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
 183         sage: K_SP.is_proper()
 184         True
 185
 186     Test the proposition in our paper concerning the duals and
 187     restrictions. Generate a random cone, then create a subcone of
 188     it. The operation of dual-taking should then commute with
 189     _restrict_to_space::
 190
 191         sage: set_random_seed()
 192         sage: J = random_cone(max_ambient_dim = 8)
 193         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 194         sage: K_W_star = _restrict_to_space(K, J.span()).dual()
 195         sage: K_star_W = _restrict_to_space(K.dual(), J.span())
 196         sage: _basically_the_same(K_W_star, K_star_W)
 197         True
 198
 199     """
 200     # First we want to intersect ``K`` with ``W``. The easiest way to
 201     # do this is via cone intersection, so we turn the subspace ``W``
 202     # into a cone.
 203     W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice())
 204     K = K.intersection(W_cone)
 205
 206     # We've already intersected K with the span of K2, so every
 207     # generator of K should belong to W now.
 208     K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 209
 210     L = ToricLattice(W.dimension())
 211     return Cone(K_W_rays, lattice=L)
 212
 213
 214
 215 def discrete_complementarity_set(K):
 216     r"""
 217     Compute the discrete complementarity set of this cone.
 218
 219     The complementarity set of a cone is the set of all orthogonal pairs
 220     `(x,s)` such that `x` is in the cone, and `s` is in its dual. The
 221     discrete complementarity set is a subset of the complementarity set
 222     where `x` and `s` are required to be generators of their respective
 223     cones.
 224
 225     For polyhedral cones, the discrete complementarity set is always
 226     finite.
 227
 228     OUTPUT:
 229
 230     A list of pairs `(x,s)` such that,
 231
 232       * Both `x` and `s` are vectors (not rays).
 233       * `x` is a generator of this cone.
 234       * `s` is a generator of this cone's dual.
 235       * `x` and `s` are orthogonal.
 236
 237     REFERENCES:
 238
 239     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 240        Improper Cone. Work in-progress.
 241
 242     EXAMPLES:
 243
 244     The discrete complementarity set of the nonnegative orthant consists
 245     of pairs of standard basis vectors::
 246
 247         sage: K = Cone([(1,0),(0,1)])
 248         sage: discrete_complementarity_set(K)
 249         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 250
 251     If the cone consists of a single ray, the second components of the
 252     discrete complementarity set should generate the orthogonal
 253     complement of that ray::
 254
 255         sage: K = Cone([(1,0)])
 256         sage: discrete_complementarity_set(K)
 257         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 258         sage: K = Cone([(1,0,0)])
 259         sage: discrete_complementarity_set(K)
 260         [((1, 0, 0), (0, 1, 0)),
 261           ((1, 0, 0), (0, -1, 0)),
 262           ((1, 0, 0), (0, 0, 1)),
 263           ((1, 0, 0), (0, 0, -1))]
 264
 265     When the cone is the entire space, its dual is the trivial cone, so
 266     the discrete complementarity set is empty::
 267
 268         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 269         sage: discrete_complementarity_set(K)
 270         []
 271
 272     Likewise when this cone is trivial (its dual is the entire space)::
 273
 274         sage: L = ToricLattice(0)
 275         sage: K = Cone([], ToricLattice(0))
 276         sage: discrete_complementarity_set(K)
 277         []
 278
 279     TESTS:
 280
 281     The complementarity set of the dual can be obtained by switching the
 282     components of the complementarity set of the original cone::
 283
 284         sage: set_random_seed()
 285         sage: K1 = random_cone(max_ambient_dim=6)
 286         sage: K2 = K1.dual()
 287         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 288         sage: actual = discrete_complementarity_set(K1)
 289         sage: sorted(actual) == sorted(expected)
 290         True
 291
 292     The pairs in the discrete complementarity set are in fact
 293     complementary::
 294
 295         sage: set_random_seed()
 296         sage: K = random_cone(max_ambient_dim=6)
 297         sage: dcs = discrete_complementarity_set(K)
 298         sage: sum([x.inner_product(s).abs() for (x,s) in dcs])
 299         0
 300
 301     """
 302     V = K.lattice().vector_space()
 303
 304     # Convert rays to vectors so that we can compute inner products.
 305     xs = [V(x) for x in K.rays()]
 306
 307     # We also convert the generators of the dual cone so that we
 308     # return pairs of vectors and not (vector, ray) pairs.
 309     ss = [V(s) for s in K.dual().rays()]
 310
 311     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 312
 313
 314 def LL(K):
 315     r"""
 316     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 317     on this cone.
 318
 319     OUTPUT:
 320
 321     A list of matrices forming a basis for the space of all
 322     Lyapunov-like transformations on the given cone.
 323
 324     EXAMPLES:
 325
 326     The trivial cone has no Lyapunov-like transformations::
 327
 328         sage: L = ToricLattice(0)
 329         sage: K = Cone([], lattice=L)
 330         sage: LL(K)
 331         []
 332
 333     The Lyapunov-like transformations on the nonnegative orthant are
 334     simply diagonal matrices::
 335
 336         sage: K = Cone([(1,)])
 337         sage: LL(K)
 338         [[1]]
 339
 340         sage: K = Cone([(1,0),(0,1)])
 341         sage: LL(K)
 342         [
 343         [1 0]  [0 0]
 344         [0 0], [0 1]
 345         ]
 346
 347         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 348         sage: LL(K)
 349         [
 350         [1 0 0]  [0 0 0]  [0 0 0]
 351         [0 0 0]  [0 1 0]  [0 0 0]
 352         [0 0 0], [0 0 0], [0 0 1]
 353         ]
 354
 355     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 356     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 357
 358         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 359         sage: LL(L31)
 360         [
 361         [1 0 0]
 362         [0 1 0]
 363         [0 0 1]
 364         ]
 365
 366         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 367         sage: LL(L3infty)
 368         [
 369         [1 0 0]
 370         [0 1 0]
 371         [0 0 1]
 372         ]
 373
 374     If our cone is the entire space, then every transformation on it is
 375     Lyapunov-like::
 376
 377         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 378         sage: M = MatrixSpace(QQ,2)
 379         sage: M.basis() == LL(K)
 380         True
 381
 382     TESTS:
 383
 384     The inner product `\left< L\left(x\right), s \right>` is zero for
 385     every pair `\left( x,s \right)` in the discrete complementarity set
 386     of the cone::
 387
 388         sage: set_random_seed()
 389         sage: K = random_cone(max_ambient_dim=8)
 390         sage: C_of_K = discrete_complementarity_set(K)
 391         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 392         sage: sum(map(abs, l))
 393         0
 394
 395     The Lyapunov-like transformations on a cone and its dual are related
 396     by transposition, but we're not guaranteed to compute transposed
 397     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 398     \right)`
 399
 400         sage: set_random_seed()
 401         sage: K = random_cone(max_ambient_dim=8)
 402         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 403         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 404         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 405         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 406         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 407         True
 408
 409     """
 410     V = K.lattice().vector_space()
 411
 412     C_of_K = discrete_complementarity_set(K)
 413
 414     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 415
 416     # Sage doesn't think matrices are vectors, so we have to convert
 417     # our matrices to vectors explicitly before we can figure out how
 418     # many are linearly-indepenedent.
 419     #
 420     # The space W has the same base ring as V, but dimension
 421     # dim(V)^2. So it has the same dimension as the space of linear
 422     # transformations on V. In other words, it's just the right size
 423     # to create an isomorphism between it and our matrices.
 424     W = VectorSpace(V.base_ring(), V.dimension()**2)
 425
 426     # Turn our matrices into long vectors...
 427     vectors = [ W(m.list()) for m in tensor_products ]
 428
 429     # Vector space representation of Lyapunov-like matrices
 430     # (i.e. vec(L) where L is Luapunov-like).
 431     LL_vector = W.span(vectors).complement()
 432
 433     # Now construct an ambient MatrixSpace in which to stick our
 434     # transformations.
 435     M = MatrixSpace(V.base_ring(), V.dimension())
 436
 437     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 438
 439     return matrix_basis
 440
 441
 442
 443 def lyapunov_rank(K):
 444     r"""
 445     Compute the Lyapunov rank (or bilinearity rank) of this cone.
 446
 447     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 448
 449     1. The dimension of the Lie algebra of the automorphism group of the
 450        cone.
 451
 452     2. The dimension of the linear space of all Lyapunov-like
 453        transformations on the cone.
 454
 455     INPUT:
 456
 457     A closed, convex polyhedral cone.
 458
 459     OUTPUT:
 460
 461     An integer representing the Lyapunov rank of the cone. If the
 462     dimension of the ambient vector space is `n`, then the Lyapunov rank
 463     will be between `1` and `n` inclusive; however a rank of `n-1` is
 464     not possible (see [Orlitzky/Gowda]_).
 465
 466     ALGORITHM:
 467
 468     The codimension formula from the second reference is used. We find
 469     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 470     and `s` are rays of our cone. It is known that these vectors are
 471     sufficient to apply the codimension formula. Once we have all such
 472     pairs, we "brute force" the codimension formula by finding all
 473     linearly-independent `xs^{T}`.
 474
 475     REFERENCES:
 476
 477     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 478        cone and Lyapunov-like transformations, Mathematical Programming, 147
 479        (2014) 155-170.
 480
 481     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 482        Improper Cone. Work in-progress.
 483
 484     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 485        optimality constraints for the cone of positive polynomials,
 486        Mathematical Programming, Series B, 129 (2011) 5-31.
 487
 488     EXAMPLES:
 489
 490     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 491     [Rudolf et al.]_::
 492
 493         sage: positives = Cone([(1,)])
 494         sage: lyapunov_rank(positives)
 495         1
 496         sage: quadrant = Cone([(1,0), (0,1)])
 497         sage: lyapunov_rank(quadrant)
 498         2
 499        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 500         sage: lyapunov_rank(octant)
 501         3
 502
 503     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 504     [Orlitzky/Gowda]_::
 505
 506         sage: R5 = VectorSpace(QQ, 5)
 507         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 508         sage: K = Cone(gs)
 509         sage: lyapunov_rank(K)
 510         25
 511
 512     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 513     [Rudolf et al.]_::
 514
 515         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 516         sage: lyapunov_rank(L31)
 517         1
 518
 519     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 520
 521         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 522         sage: lyapunov_rank(L3infty)
 523         1
 524
 525     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 526     + 1` [Orlitzky/Gowda]_::
 527
 528         sage: K = Cone([(1,0,0,0,0)])
 529         sage: lyapunov_rank(K)
 530         21
 531         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 532         21
 533
 534     A subspace (of dimension `m`) in `n` dimensions should have a
 535     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 536
 537         sage: e1 = (1,0,0,0,0)
 538         sage: neg_e1 = (-1,0,0,0,0)
 539         sage: e2 = (0,1,0,0,0)
 540         sage: neg_e2 = (0,-1,0,0,0)
 541         sage: z = (0,0,0,0,0)
 542         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 543         sage: lyapunov_rank(K)
 544         19
 545         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 546         19
 547
 548     The Lyapunov rank should be additive on a product of proper cones
 549     [Rudolf et al.]_::
 550
 551         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 552         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 553         sage: K = L31.cartesian_product(octant)
 554         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 555         True
 556
 557     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 558     The cone ``K`` in the following example is isomorphic to the nonnegative
 559     octant in `\mathbb{R}^{3}`::
 560
 561         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 562         sage: lyapunov_rank(K)
 563         3
 564
 565     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 566     itself [Rudolf et al.]_::
 567
 568         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 569         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 570         True
 571
 572     TESTS:
 573
 574     The Lyapunov rank should be additive on a product of proper cones
 575     [Rudolf et al.]_::
 576
 577         sage: set_random_seed()
 578         sage: K1 = random_cone(max_ambient_dim=8,
 579         ....:                  strictly_convex=True,
 580         ....:                  solid=True)
 581         sage: K2 = random_cone(max_ambient_dim=8,
 582         ....:                  strictly_convex=True,
 583         ....:                  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 Lyapunov rank is invariant under a linear isomorphism
 589     [Orlitzky/Gowda]_::
 590
 591         sage: K1 = random_cone(max_ambient_dim = 8)
 592         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 593         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 594         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 595         True
 596
 597     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 598     itself [Rudolf et al.]_::
 599
 600         sage: set_random_seed()
 601         sage: K = random_cone(max_ambient_dim=8)
 602         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 603         True
 604
 605     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 606     be any number between `1` and `n` inclusive, excluding `n-1`
 607     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 608     trivial cone in a trivial space as well. However, in zero dimensions,
 609     the Lyapunov rank of the trivial cone will be zero::
 610
 611         sage: set_random_seed()
 612         sage: K = random_cone(max_ambient_dim=8,
 613         ....:                 strictly_convex=True,
 614         ....:                 solid=True)
 615         sage: b = lyapunov_rank(K)
 616         sage: n = K.lattice_dim()
 617         sage: (n == 0 or 1 <= b) and b <= n
 618         True
 619         sage: b == n-1
 620         False
 621
 622     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 623     Lyapunov rank `n-1` in `n` dimensions::
 624
 625         sage: set_random_seed()
 626         sage: K = random_cone(max_ambient_dim=8)
 627         sage: b = lyapunov_rank(K)
 628         sage: n = K.lattice_dim()
 629         sage: b == n-1
 630         False
 631
 632     The calculation of the Lyapunov rank of an improper cone can be
 633     reduced to that of a proper cone [Orlitzky/Gowda]_::
 634
 635         sage: set_random_seed()
 636         sage: K = random_cone(max_ambient_dim=8)
 637         sage: actual = lyapunov_rank(K)
 638         sage: K_S = _restrict_to_space(K, K.span())
 639         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 640         sage: l = K.lineality()
 641         sage: c = K.codim()
 642         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 643         sage: actual == expected
 644         True
 645
 646     The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
 647
 648         sage: set_random_seed()
 649         sage: K = random_cone(max_ambient_dim=8)
 650         sage: lyapunov_rank(K) == len(LL(K))
 651         True
 652
 653     We can make an imperfect cone perfect by adding a slack variable
 654     (a Theorem in [Orlitzky/Gowda]_)::
 655
 656         sage: set_random_seed()
 657         sage: K = random_cone(max_ambient_dim=8,
 658         ....:                 strictly_convex=True,
 659         ....:                 solid=True)
 660         sage: L = ToricLattice(K.lattice_dim() + 1)
 661         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 662         sage: lyapunov_rank(K) >= K.lattice_dim()
 663         True
 664
 665     """
 666     beta = 0
 667
 668     m = K.dim()
 669     n = K.lattice_dim()
 670     l = K.lineality()
 671
 672     if m < n:
 673         # K is not solid, restrict to its span.
 674         K = _restrict_to_space(K, K.span())
 675
 676         # Non-solid reduction lemma.
 677         beta += (n - m)*n
 678
 679     if l > 0:
 680         # K is not pointed, restrict to the span of its dual. Uses a
 681         # proposition from our paper, i.e. this is equivalent to K =
 682         # _rho(K.dual()).dual().
 683         K = _restrict_to_space(K, K.dual().span())
 684
 685         # Non-pointed reduction lemma.
 686         beta += l * m
 687
 688     beta += len(LL(K))
 689     return beta