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