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 _rho(K, K2=None):
  95     r"""
  96     Restrict ``K`` into its own span, or the span of another cone.
  97
  98     INPUT:
  99
 100     - ``K2`` -- another cone whose lattice has the same rank as this
 101                 cone.
 102
 103     OUTPUT:
 104
 105     A new cone in a sublattice.
 106
 107     EXAMPLES::
 108
 109         sage: K = Cone([(1,)])
 110         sage: _rho(K) == K
 111         True
 112
 113         sage: K2 = Cone([(1,0)])
 114         sage: _rho(K2).rays()
 115         N(1)
 116         in 1-d lattice N
 117         sage: K3 = Cone([(1,0,0)])
 118         sage: _rho(K3).rays()
 119         N(1)
 120         in 1-d lattice N
 121         sage: _rho(K2) == _rho(K3)
 122         True
 123
 124     TESTS:
 125
 126     The projected cone should always be solid::
 127
 128         sage: set_random_seed()
 129         sage: K = random_cone(max_ambient_dim = 8)
 130         sage: K_S = _rho(K)
 131         sage: K_S.is_solid()
 132         True
 133
 134     And the resulting cone should live in a space having the same
 135     dimension as the space we restricted it to::
 136
 137         sage: set_random_seed()
 138         sage: K = random_cone(max_ambient_dim = 8)
 139         sage: K_S = _rho(K, K.dual() )
 140         sage: K_S.lattice_dim() == K.dual().dim()
 141         True
 142
 143     This function should not affect the dimension of a cone::
 144
 145         sage: set_random_seed()
 146         sage: K = random_cone(max_ambient_dim = 8)
 147         sage: K.dim() == _rho(K).dim()
 148         True
 149
 150     Nor should it affect the lineality of a cone::
 151
 152         sage: set_random_seed()
 153         sage: K = random_cone(max_ambient_dim = 8)
 154         sage: K.lineality() == _rho(K).lineality()
 155         True
 156
 157     No matter which space we restrict to, the lineality should not
 158     increase::
 159
 160         sage: set_random_seed()
 161         sage: K = random_cone(max_ambient_dim = 8)
 162         sage: K.lineality() >= _rho(K).lineality()
 163         True
 164         sage: K.lineality() >= _rho(K, K.dual()).lineality()
 165         True
 166
 167     If we do this according to our paper, then the result is proper::
 168
 169         sage: set_random_seed()
 170         sage: K = random_cone(max_ambient_dim = 8)
 171         sage: K_S = _rho(K)
 172         sage: K_SP = _rho(K_S.dual()).dual()
 173         sage: K_SP.is_proper()
 174         True
 175         sage: K_SP = _rho(K_S, K_S.dual())
 176         sage: K_SP.is_proper()
 177         True
 178
 179     Test the proposition in our paper concerning the duals and
 180     restrictions. Generate a random cone, then create a subcone of
 181     it. The operation of dual-taking should then commute with rho::
 182
 183         sage: set_random_seed()
 184         sage: J = random_cone(max_ambient_dim = 8)
 185         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 186         sage: K_W_star = _rho(K, J).dual()
 187         sage: K_star_W = _rho(K.dual(), J)
 188         sage: _basically_the_same(K_W_star, K_star_W)
 189         True
 190
 191     """
 192     if K2 is None:
 193         K2 = K
 194
 195     # First we project K onto the span of K2. This will explode if the
 196     # rank of ``K2.lattice()`` doesn't match ours.
 197     span_K2 = Cone(K2.rays() + (-K2).rays(), lattice=K.lattice())
 198     K = K.intersection(span_K2)
 199
 200     # Cheat a little to get the subspace span(K2). The paper uses the
 201     # rays of K2 as a basis, but everything is invariant under linear
 202     # isomorphism (i.e. a change of basis), and this is a little
 203     # faster.
 204     W = span_K2.linear_subspace()
 205
 206     # We've already intersected K with the span of K2, so every
 207     # generator of K should belong to W now.
 208     W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 209
 210     L = ToricLattice(K2.dim())
 211     return Cone(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 (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 the first reference).
 465
 466     .. note::
 467
 468         In the references, the cones are always assumed to be proper. We
 469         do not impose this restriction.
 470
 471     .. seealso::
 472
 473         :meth:`is_proper`
 474
 475     ALGORITHM:
 476
 477     The codimension formula from the second reference is used. We find
 478     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 479     and `s` are rays of our cone. It is known that these vectors are
 480     sufficient to apply the codimension formula. Once we have all such
 481     pairs, we "brute force" the codimension formula by finding all
 482     linearly-independent `xs^{T}`.
 483
 484     REFERENCES:
 485
 486     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 487        cone and Lyapunov-like transformations, Mathematical Programming, 147
 488        (2014) 155-170.
 489
 490     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 491        Improper Cone. Work in-progress.
 492
 493     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 494        optimality constraints for the cone of positive polynomials,
 495        Mathematical Programming, Series B, 129 (2011) 5-31.
 496
 497     EXAMPLES:
 498
 499     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 500     [Rudolf et al.]_::
 501
 502         sage: positives = Cone([(1,)])
 503         sage: lyapunov_rank(positives)
 504         1
 505         sage: quadrant = Cone([(1,0), (0,1)])
 506         sage: lyapunov_rank(quadrant)
 507         2
 508        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 509         sage: lyapunov_rank(octant)
 510         3
 511
 512     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 513     [Orlitzky/Gowda]_::
 514
 515         sage: R5 = VectorSpace(QQ, 5)
 516         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 517         sage: K = Cone(gs)
 518         sage: lyapunov_rank(K)
 519         25
 520
 521     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 522     [Rudolf et al.]_::
 523
 524         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 525         sage: lyapunov_rank(L31)
 526         1
 527
 528     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 529
 530         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 531         sage: lyapunov_rank(L3infty)
 532         1
 533
 534     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 535     + 1` [Orlitzky/Gowda]_::
 536
 537         sage: K = Cone([(1,0,0,0,0)])
 538         sage: lyapunov_rank(K)
 539         21
 540         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 541         21
 542
 543     A subspace (of dimension `m`) in `n` dimensions should have a
 544     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 545
 546         sage: e1 = (1,0,0,0,0)
 547         sage: neg_e1 = (-1,0,0,0,0)
 548         sage: e2 = (0,1,0,0,0)
 549         sage: neg_e2 = (0,-1,0,0,0)
 550         sage: z = (0,0,0,0,0)
 551         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 552         sage: lyapunov_rank(K)
 553         19
 554         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 555         19
 556
 557     The Lyapunov rank should be additive on a product of proper cones
 558     [Rudolf et al.]_::
 559
 560         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 561         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 562         sage: K = L31.cartesian_product(octant)
 563         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 564         True
 565
 566     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 567     The cone ``K`` in the following example is isomorphic to the nonnegative
 568     octant in `\mathbb{R}^{3}`::
 569
 570         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 571         sage: lyapunov_rank(K)
 572         3
 573
 574     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 575     itself [Rudolf et al.]_::
 576
 577         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 578         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 579         True
 580
 581     TESTS:
 582
 583     The Lyapunov rank should be additive on a product of proper cones
 584     [Rudolf et al.]_::
 585
 586         sage: set_random_seed()
 587         sage: K1 = random_cone(max_ambient_dim=8,
 588         ....:                  strictly_convex=True,
 589         ....:                  solid=True)
 590         sage: K2 = random_cone(max_ambient_dim=8,
 591         ....:                  strictly_convex=True,
 592         ....:                  solid=True)
 593         sage: K = K1.cartesian_product(K2)
 594         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 595         True
 596
 597     The Lyapunov rank is invariant under a linear isomorphism
 598     [Orlitzky/Gowda]_::
 599
 600         sage: K1 = random_cone(max_ambient_dim = 8)
 601         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 602         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 603         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 604         True
 605
 606     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 607     itself [Rudolf et al.]_::
 608
 609         sage: set_random_seed()
 610         sage: K = random_cone(max_ambient_dim=8)
 611         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 612         True
 613
 614     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 615     be any number between `1` and `n` inclusive, excluding `n-1`
 616     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 617     trivial cone in a trivial space as well. However, in zero dimensions,
 618     the Lyapunov rank of the trivial cone will be zero::
 619
 620         sage: set_random_seed()
 621         sage: K = random_cone(max_ambient_dim=8,
 622         ....:                 strictly_convex=True,
 623         ....:                 solid=True)
 624         sage: b = lyapunov_rank(K)
 625         sage: n = K.lattice_dim()
 626         sage: (n == 0 or 1 <= b) and b <= n
 627         True
 628         sage: b == n-1
 629         False
 630
 631     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 632     Lyapunov rank `n-1` in `n` dimensions::
 633
 634         sage: set_random_seed()
 635         sage: K = random_cone(max_ambient_dim=8)
 636         sage: b = lyapunov_rank(K)
 637         sage: n = K.lattice_dim()
 638         sage: b == n-1
 639         False
 640
 641     The calculation of the Lyapunov rank of an improper cone can be
 642     reduced to that of a proper cone [Orlitzky/Gowda]_::
 643
 644         sage: set_random_seed()
 645         sage: K = random_cone(max_ambient_dim=8)
 646         sage: actual = lyapunov_rank(K)
 647         sage: K_S = _rho(K)
 648         sage: K_SP = _rho(K_S.dual()).dual()
 649         sage: l = K.lineality()
 650         sage: c = K.codim()
 651         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 652         sage: actual == expected
 653         True
 654
 655     The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
 656
 657         sage: set_random_seed()
 658         sage: K = random_cone(max_ambient_dim=8)
 659         sage: lyapunov_rank(K) == len(LL(K))
 660         True
 661
 662     Test Theorem 3 in [Orlitzky/Gowda]_::
 663
 664         sage: set_random_seed()
 665         sage: K = random_cone(max_ambient_dim=8,
 666         ....:                 strictly_convex=True,
 667         ....:                 solid=True)
 668         sage: L = ToricLattice(K.lattice_dim() + 1)
 669         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 670         sage: lyapunov_rank(K) >= K.lattice_dim()
 671         True
 672
 673     """
 674     beta = 0
 675
 676     m = K.dim()
 677     n = K.lattice_dim()
 678     l = K.lineality()
 679
 680     if m < n:
 681         # K is not solid, restrict to its span.
 682         K = _rho(K)
 683
 684         # Lemma 2
 685         beta += m*(n - m) + (n - m)**2
 686
 687     if l > 0:
 688         # K is not pointed, restrict to the span of its dual. Uses a
 689         # proposition from our paper, i.e. this is equivalent to K =
 690         # _rho(K.dual()).dual().
 691         K = _rho(K, K.dual())
 692
 693         # Lemma 3
 694         beta += m * l
 695
 696     beta += len(LL(K))
 697     return beta