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