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