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 drop_dependent(vs):
  12     r"""
  13     Return the largest linearly-independent subset of ``vs``.
  14     """
  15     if len(vs) == 0:
  16         # ...for lazy enough definitions of linearly-independent
  17         return vs
  18
  19     result = []
  20     old_V = VectorSpace(vs[0].parent().base_field(), 0)
  21
  22     for v in vs:
  23         new_V = span(result + [v])
  24         if new_V.dimension() > old_V.dimension():
  25             result.append(v)
  26             old_V = new_V
  27
  28     return result
  29
  30
  31 def basically_the_same(K1,K2):
  32     r"""
  33     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  34     otherwise.
  35     """
  36     if K1.lattice_dim() != K2.lattice_dim():
  37         return False
  38
  39     if K1.nrays() != K2.nrays():
  40         return False
  41
  42     if K1.dim() != K2.dim():
  43         return False
  44
  45     if lineality(K1) != lineality(K2):
  46         return False
  47
  48     if K1.is_solid() != K2.is_solid():
  49         return False
  50
  51     if K1.is_strictly_convex() != K2.is_strictly_convex():
  52         return False
  53
  54     if len(LL(K1)) != len(LL(K2)):
  55         return False
  56
  57     C_of_K1 = discrete_complementarity_set(K1)
  58     C_of_K2 = discrete_complementarity_set(K2)
  59     if len(C_of_K1) != len(C_of_K2):
  60         return False
  61
  62     if len(K1.facets()) != len(K2.facets()):
  63         return False
  64
  65     return True
  66
  67
  68
  69 def iso_space(K):
  70     r"""
  71     Construct the space `W \times W^{\perp}` isomorphic to the ambient space
  72     of ``K`` where `W` is equal to the span of ``K``.
  73     """
  74     V = K.lattice().vector_space()
  75
  76     # Create the space W \times W^{\perp} isomorphic to V.
  77     # First we get an orthogonal (but not normal) basis...
  78     M = matrix(V.base_field(), K.rays())
  79     W_basis = drop_dependent(K.rays())
  80
  81     W = V.subspace_with_basis(W_basis)
  82     W_perp = W.complement()
  83
  84     return W.cartesian_product(W_perp)
  85
  86
  87 def ips_iso(K):
  88     r"""
  89     Construct the IPS isomorphism and its inverse from our paper.
  90
  91     Given a cone ``K``, the returned isomorphism will split its ambient
  92     vector space `V` into a cartesian product `W \times W^{\perp}` where
  93     `W` equals the span of ``K``.
  94     """
  95     V = K.lattice().vector_space()
  96     V_iso = iso_space(K)
  97     (W, W_perp) = V_iso.cartesian_factors()
  98
  99     # A space equivalent to V, but using our basis.
 100     V_user = V.subspace_with_basis( W.basis() + W_perp.basis() )
 101
 102     def phi(v):
 103         # Write v in terms of our custom basis, where the first dim(W)
 104         # coordinates are for the W-part of the basis.
 105         cs = V_user.coordinates(v)
 106
 107         w1 = sum([ V_user.basis()[idx]*cs[idx]
 108                     for idx in range(0, W.dimension()) ])
 109         w2 = sum([ V_user.basis()[idx]*cs[idx]
 110                     for idx in range(W.dimension(), V.dimension()) ])
 111
 112         return V_iso( (w1, w2) )
 113
 114
 115     def phi_inv( pair ):
 116         # Crash if the arguments are in the wrong spaces.
 117         V_iso(pair)
 118
 119         #w = sum([ sub_w[idx]*W.basis()[idx] for idx in range(0,m) ])
 120         #w_prime = sum([ sub_w_prime[idx]*W_perp.basis()[idx]
 121         #             for idx in range(0,n-m) ])
 122
 123         return sum( pair.cartesian_factors() )
 124
 125
 126     return (phi,phi_inv)
 127
 128
 129 def rho(K, K2=None):
 130     r"""
 131     Restrict ``K`` into its own span, or the span of another cone.
 132
 133     INPUT:
 134
 135     - ``K2`` -- another cone whose lattice has the same rank as this cone.
 136
 137     OUTPUT:
 138
 139     A new cone in a sublattice.
 140
 141     EXAMPLES::
 142
 143         sage: K = Cone([(1,)])
 144         sage: restrict_span(K) == K
 145         True
 146
 147         sage: K2 = Cone([(1,0)])
 148         sage: restrict_span(K2).rays()
 149         N(1)
 150         in 1-d lattice N
 151         sage: K3 = Cone([(1,0,0)])
 152         sage: restrict_span(K3).rays()
 153         N(1)
 154         in 1-d lattice N
 155         sage: restrict_span(K2) == restrict_span(K3)
 156         True
 157
 158     TESTS:
 159
 160     The projected cone should always be solid::
 161
 162         sage: set_random_seed()
 163         sage: K = random_cone(max_dim = 8)
 164         sage: K_S = restrict_span(K)
 165         sage: K_S.is_solid()
 166         True
 167
 168     And the resulting cone should live in a space having the same
 169     dimension as the space we restricted it to::
 170
 171         sage: set_random_seed()
 172         sage: K = random_cone(max_dim = 8)
 173         sage: K_S = restrict_span(K, K.dual() )
 174         sage: K_S.lattice_dim() == K.dual().dim()
 175         True
 176
 177     This function has ``unrestrict_span()`` as its inverse::
 178
 179         sage: set_random_seed()
 180         sage: K = random_cone(max_dim = 8, solid=True)
 181         sage: J = restrict_span(K)
 182         sage: K == unrestrict_span(J,K)
 183         True
 184
 185     This function should not affect the dimension of a cone::
 186
 187         sage: set_random_seed()
 188         sage: K = random_cone(max_dim = 8)
 189         sage: K.dim() == restrict_span(K).dim()
 190         True
 191
 192     Nor should it affect the lineality of a cone::
 193
 194         sage: set_random_seed()
 195         sage: K = random_cone(max_dim = 8)
 196         sage: lineality(K) == lineality(restrict_span(K))
 197         True
 198
 199     No matter which space we restrict to, the lineality should not
 200     increase::
 201
 202         sage: set_random_seed()
 203         sage: K = random_cone(max_dim = 8)
 204         sage: lineality(K) >= lineality(restrict_span(K))
 205         True
 206         sage: lineality(K) >= lineality(restrict_span(K, K.dual()))
 207         True
 208
 209     If we do this according to our paper, then the result is proper::
 210
 211         sage: set_random_seed()
 212         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
 213         sage: K_S = restrict_span(K)
 214         sage: P = restrict_span(K_S.dual()).dual()
 215         sage: P.is_proper()
 216         True
 217         sage: P = restrict_span(K_S, K_S.dual())
 218         sage: P.is_proper()
 219         True
 220
 221     ::
 222
 223         sage: set_random_seed()
 224         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
 225         sage: K_S = restrict_span(K)
 226         sage: P = restrict_span(K_S.dual()).dual()
 227         sage: P.is_proper()
 228         True
 229         sage: P = restrict_span(K_S, K_S.dual())
 230         sage: P.is_proper()
 231         True
 232
 233     ::
 234
 235         sage: set_random_seed()
 236         sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
 237         sage: K_S = restrict_span(K)
 238         sage: P = restrict_span(K_S.dual()).dual()
 239         sage: P.is_proper()
 240         True
 241         sage: P = restrict_span(K_S, K_S.dual())
 242         sage: P.is_proper()
 243         True
 244
 245     ::
 246
 247         sage: set_random_seed()
 248         sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
 249         sage: K_S = restrict_span(K)
 250         sage: P = restrict_span(K_S.dual()).dual()
 251         sage: P.is_proper()
 252         True
 253         sage: P = restrict_span(K_S, K_S.dual())
 254         sage: P.is_proper()
 255         True
 256
 257     Test the proposition in our paper concerning the duals, where the
 258     subspace `W` is the span of `K^{*}`::
 259
 260         sage: set_random_seed()
 261         sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=False)
 262         sage: K_W = restrict_span(K, K.dual())
 263         sage: K_star_W_star = restrict_span(K.dual()).dual()
 264         sage: basically_the_same(K_W, K_star_W_star)
 265         True
 266
 267     ::
 268
 269         sage: set_random_seed()
 270         sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=False)
 271         sage: K_W = restrict_span(K, K.dual())
 272         sage: K_star_W_star = restrict_span(K.dual()).dual()
 273         sage: basically_the_same(K_W, K_star_W_star)
 274         True
 275
 276     ::
 277
 278         sage: set_random_seed()
 279         sage: K = random_cone(max_dim = 8, solid=False, strictly_convex=True)
 280         sage: K_W = restrict_span(K, K.dual())
 281         sage: K_star_W_star = restrict_span(K.dual()).dual()
 282         sage: basically_the_same(K_W, K_star_W_star)
 283         True
 284
 285     ::
 286
 287         sage: set_random_seed()
 288         sage: K = random_cone(max_dim = 8, solid=True, strictly_convex=True)
 289         sage: K_W = restrict_span(K, K.dual())
 290         sage: K_star_W_star = restrict_span(K.dual()).dual()
 291         sage: basically_the_same(K_W, K_star_W_star)
 292         True
 293
 294     """
 295     if K2 is None:
 296         K2 = K
 297
 298     phi,_ = ips_iso(K2)
 299     (W, W_perp) = iso_space(K2).cartesian_factors()
 300
 301     ray_pairs = [ phi(r) for r in K.rays() ]
 302
 303     # Shouldn't matter?
 304     #
 305     #if any([ w2 != W_perp.zero() for (_, w2) in ray_pairs ]):
 306     #    msg = 'Cone has nonzero components in W-perp!'
 307     #    raise ValueError(msg)
 308
 309     # Represent the cone in terms of a basis for W, i.e. with smaller
 310     # vectors.
 311     ws = [ W.coordinate_vector(w1) for (w1, _) in ray_pairs ]
 312
 313     L = ToricLattice(W.dimension())
 314
 315     return Cone(ws, lattice=L)
 316
 317
 318
 319 def lineality(K):
 320     r"""
 321     Compute the lineality of this cone.
 322
 323     The lineality of a cone is the dimension of the largest linear
 324     subspace contained in that cone.
 325
 326     OUTPUT:
 327
 328     A nonnegative integer; the dimension of the largest subspace
 329     contained within this cone.
 330
 331     REFERENCES:
 332
 333     .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
 334        University Press, Princeton, 1970.
 335
 336     EXAMPLES:
 337
 338     The lineality of the nonnegative orthant is zero, since it clearly
 339     contains no lines::
 340
 341         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
 342         sage: lineality(K)
 343         0
 344
 345     However, if we add another ray so that the entire `x`-axis belongs
 346     to the cone, then the resulting cone will have lineality one::
 347
 348         sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
 349         sage: lineality(K)
 350         1
 351
 352     If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
 353     to the dimension of the ambient space (i.e. two)::
 354
 355         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 356         sage: lineality(K)
 357         2
 358
 359     Per the definition, the lineality of the trivial cone in a trivial
 360     space is zero::
 361
 362         sage: K = Cone([], lattice=ToricLattice(0))
 363         sage: lineality(K)
 364         0
 365
 366     TESTS:
 367
 368     The lineality of a cone should be an integer between zero and the
 369     dimension of the ambient space, inclusive::
 370
 371         sage: set_random_seed()
 372         sage: K = random_cone(max_dim = 8)
 373         sage: l = lineality(K)
 374         sage: l in ZZ
 375         True
 376         sage: (0 <= l) and (l <= K.lattice_dim())
 377         True
 378
 379     A strictly convex cone should have lineality zero::
 380
 381         sage: set_random_seed()
 382         sage: K = random_cone(max_dim = 8, strictly_convex = True)
 383         sage: lineality(K)
 384         0
 385
 386     """
 387     return K.linear_subspace().dimension()
 388
 389
 390 def discrete_complementarity_set(K):
 391     r"""
 392     Compute the discrete complementarity set of this cone.
 393
 394     The complementarity set of this cone is the set of all orthogonal
 395     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 396     dual. The discrete complementarity set restricts `x` and `s` to be
 397     generators of their respective cones.
 398
 399     OUTPUT:
 400
 401     A list of pairs `(x,s)` such that,
 402
 403       * `x` is in this cone.
 404       * `x` is a generator of this cone.
 405       * `s` is in this cone's dual.
 406       * `s` is a generator of this cone's dual.
 407       * `x` and `s` are orthogonal.
 408
 409     EXAMPLES:
 410
 411     The discrete complementarity set of the nonnegative orthant consists
 412     of pairs of standard basis vectors::
 413
 414         sage: K = Cone([(1,0),(0,1)])
 415         sage: discrete_complementarity_set(K)
 416         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 417
 418     If the cone consists of a single ray, the second components of the
 419     discrete complementarity set should generate the orthogonal
 420     complement of that ray::
 421
 422         sage: K = Cone([(1,0)])
 423         sage: discrete_complementarity_set(K)
 424         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 425         sage: K = Cone([(1,0,0)])
 426         sage: discrete_complementarity_set(K)
 427         [((1, 0, 0), (0, 1, 0)),
 428           ((1, 0, 0), (0, -1, 0)),
 429           ((1, 0, 0), (0, 0, 1)),
 430           ((1, 0, 0), (0, 0, -1))]
 431
 432     When the cone is the entire space, its dual is the trivial cone, so
 433     the discrete complementarity set is empty::
 434
 435         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 436         sage: discrete_complementarity_set(K)
 437         []
 438
 439     TESTS:
 440
 441     The complementarity set of the dual can be obtained by switching the
 442     components of the complementarity set of the original cone::
 443
 444         sage: set_random_seed()
 445         sage: K1 = random_cone(max_dim=6)
 446         sage: K2 = K1.dual()
 447         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 448         sage: actual = discrete_complementarity_set(K1)
 449         sage: sorted(actual) == sorted(expected)
 450         True
 451
 452     """
 453     V = K.lattice().vector_space()
 454
 455     # Convert the rays to vectors so that we can compute inner
 456     # products.
 457     xs = [V(x) for x in K.rays()]
 458     ss = [V(s) for s in K.dual().rays()]
 459
 460     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 461
 462
 463 def LL(K):
 464     r"""
 465     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 466     on this cone.
 467
 468     OUTPUT:
 469
 470     A list of matrices forming a basis for the space of all
 471     Lyapunov-like transformations on the given cone.
 472
 473     EXAMPLES:
 474
 475     The trivial cone has no Lyapunov-like transformations::
 476
 477         sage: L = ToricLattice(0)
 478         sage: K = Cone([], lattice=L)
 479         sage: LL(K)
 480         []
 481
 482     The Lyapunov-like transformations on the nonnegative orthant are
 483     simply diagonal matrices::
 484
 485         sage: K = Cone([(1,)])
 486         sage: LL(K)
 487         [[1]]
 488
 489         sage: K = Cone([(1,0),(0,1)])
 490         sage: LL(K)
 491         [
 492         [1 0]  [0 0]
 493         [0 0], [0 1]
 494         ]
 495
 496         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 497         sage: LL(K)
 498         [
 499         [1 0 0]  [0 0 0]  [0 0 0]
 500         [0 0 0]  [0 1 0]  [0 0 0]
 501         [0 0 0], [0 0 0], [0 0 1]
 502         ]
 503
 504     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 505     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 506
 507         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 508         sage: LL(L31)
 509         [
 510         [1 0 0]
 511         [0 1 0]
 512         [0 0 1]
 513         ]
 514
 515         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 516         sage: LL(L3infty)
 517         [
 518         [1 0 0]
 519         [0 1 0]
 520         [0 0 1]
 521         ]
 522
 523     If our cone is the entire space, then every transformation on it is
 524     Lyapunov-like::
 525
 526         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 527         sage: M = MatrixSpace(QQ,2)
 528         sage: M.basis() == LL(K)
 529         True
 530
 531     TESTS:
 532
 533     The inner product `\left< L\left(x\right), s \right>` is zero for
 534     every pair `\left( x,s \right)` in the discrete complementarity set
 535     of the cone::
 536
 537         sage: set_random_seed()
 538         sage: K = random_cone(max_dim=8)
 539         sage: C_of_K = discrete_complementarity_set(K)
 540         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 541         sage: sum(map(abs, l))
 542         0
 543
 544     The Lyapunov-like transformations on a cone and its dual are related
 545     by transposition, but we're not guaranteed to compute transposed
 546     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 547     \right)`
 548
 549         sage: set_random_seed()
 550         sage: K = random_cone(max_dim=8)
 551         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 552         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 553         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 554         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 555         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 556         True
 557
 558     """
 559     V = K.lattice().vector_space()
 560
 561     C_of_K = discrete_complementarity_set(K)
 562
 563     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 564
 565     # Sage doesn't think matrices are vectors, so we have to convert
 566     # our matrices to vectors explicitly before we can figure out how
 567     # many are linearly-indepenedent.
 568     #
 569     # The space W has the same base ring as V, but dimension
 570     # dim(V)^2. So it has the same dimension as the space of linear
 571     # transformations on V. In other words, it's just the right size
 572     # to create an isomorphism between it and our matrices.
 573     W = VectorSpace(V.base_ring(), V.dimension()**2)
 574
 575     # Turn our matrices into long vectors...
 576     vectors = [ W(m.list()) for m in tensor_products ]
 577
 578     # Vector space representation of Lyapunov-like matrices
 579     # (i.e. vec(L) where L is Luapunov-like).
 580     LL_vector = W.span(vectors).complement()
 581
 582     # Now construct an ambient MatrixSpace in which to stick our
 583     # transformations.
 584     M = MatrixSpace(V.base_ring(), V.dimension())
 585
 586     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 587
 588     return matrix_basis
 589
 590
 591
 592 def lyapunov_rank(K):
 593     r"""
 594     Compute the Lyapunov (or bilinearity) rank of this cone.
 595
 596     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 597
 598     1. The dimension of the Lie algebra of the automorphism group of the
 599        cone.
 600
 601     2. The dimension of the linear space of all Lyapunov-like
 602        transformations on the cone.
 603
 604     INPUT:
 605
 606     A closed, convex polyhedral cone.
 607
 608     OUTPUT:
 609
 610     An integer representing the Lyapunov rank of the cone. If the
 611     dimension of the ambient vector space is `n`, then the Lyapunov rank
 612     will be between `1` and `n` inclusive; however a rank of `n-1` is
 613     not possible (see the first reference).
 614
 615     .. note::
 616
 617         In the references, the cones are always assumed to be proper. We
 618         do not impose this restriction.
 619
 620     .. seealso::
 621
 622         :meth:`is_proper`
 623
 624     ALGORITHM:
 625
 626     The codimension formula from the second reference is used. We find
 627     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 628     and `s` are rays of our cone. It is known that these vectors are
 629     sufficient to apply the codimension formula. Once we have all such
 630     pairs, we "brute force" the codimension formula by finding all
 631     linearly-independent `xs^{T}`.
 632
 633     REFERENCES:
 634
 635     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 636        cone and Lyapunov-like transformations, Mathematical Programming, 147
 637        (2014) 155-170.
 638
 639     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 640        Improper Cone. Work in-progress.
 641
 642     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 643        optimality constraints for the cone of positive polynomials,
 644        Mathematical Programming, Series B, 129 (2011) 5-31.
 645
 646     EXAMPLES:
 647
 648     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 649     [Rudolf et al.]_::
 650
 651         sage: positives = Cone([(1,)])
 652         sage: lyapunov_rank(positives)
 653         1
 654         sage: quadrant = Cone([(1,0), (0,1)])
 655         sage: lyapunov_rank(quadrant)
 656         2
 657        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 658         sage: lyapunov_rank(octant)
 659         3
 660
 661     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 662     [Orlitzky/Gowda]_::
 663
 664         sage: R5 = VectorSpace(QQ, 5)
 665         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 666         sage: K = Cone(gs)
 667         sage: lyapunov_rank(K)
 668         25
 669
 670     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 671     [Rudolf et al.]_::
 672
 673         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 674         sage: lyapunov_rank(L31)
 675         1
 676
 677     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 678
 679         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 680         sage: lyapunov_rank(L3infty)
 681         1
 682
 683     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 684     + 1` [Orlitzky/Gowda]_::
 685
 686         sage: K = Cone([(1,0,0,0,0)])
 687         sage: lyapunov_rank(K)
 688         21
 689         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 690         21
 691
 692     A subspace (of dimension `m`) in `n` dimensions should have a
 693     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 694
 695         sage: e1 = (1,0,0,0,0)
 696         sage: neg_e1 = (-1,0,0,0,0)
 697         sage: e2 = (0,1,0,0,0)
 698         sage: neg_e2 = (0,-1,0,0,0)
 699         sage: z = (0,0,0,0,0)
 700         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 701         sage: lyapunov_rank(K)
 702         19
 703         sage: K.lattice_dim()**2 - K.dim()*codim(K)
 704         19
 705
 706     The Lyapunov rank should be additive on a product of proper cones
 707     [Rudolf et al.]_::
 708
 709         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 710         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 711         sage: K = L31.cartesian_product(octant)
 712         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 713         True
 714
 715     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 716     The cone ``K`` in the following example is isomorphic to the nonnegative
 717     octant in `\mathbb{R}^{3}`::
 718
 719         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 720         sage: lyapunov_rank(K)
 721         3
 722
 723     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 724     itself [Rudolf et al.]_::
 725
 726         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 727         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 728         True
 729
 730     TESTS:
 731
 732     The Lyapunov rank should be additive on a product of proper cones
 733     [Rudolf et al.]_::
 734
 735         sage: set_random_seed()
 736         sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 737         sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 738         sage: K = K1.cartesian_product(K2)
 739         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 740         True
 741
 742     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 743     itself [Rudolf et al.]_::
 744
 745         sage: set_random_seed()
 746         sage: K = random_cone(max_dim=8)
 747         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 748         True
 749
 750     Make sure we exercise the non-strictly-convex/non-solid case::
 751
 752         sage: set_random_seed()
 753         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 754         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 755         True
 756
 757     Let's check the other permutations as well, just to be sure::
 758
 759         sage: set_random_seed()
 760         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 761         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 762         True
 763
 764     ::
 765
 766         sage: set_random_seed()
 767         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 768         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 769         True
 770
 771     ::
 772
 773         sage: set_random_seed()
 774         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 775         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 776         True
 777
 778     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 779     be any number between `1` and `n` inclusive, excluding `n-1`
 780     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 781     trivial cone in a trivial space as well. However, in zero dimensions,
 782     the Lyapunov rank of the trivial cone will be zero::
 783
 784         sage: set_random_seed()
 785         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 786         sage: b = lyapunov_rank(K)
 787         sage: n = K.lattice_dim()
 788         sage: (n == 0 or 1 <= b) and b <= n
 789         True
 790         sage: b == n-1
 791         False
 792
 793     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 794     Lyapunov rank `n-1` in `n` dimensions::
 795
 796         sage: set_random_seed()
 797         sage: K = random_cone(max_dim=8)
 798         sage: b = lyapunov_rank(K)
 799         sage: n = K.lattice_dim()
 800         sage: b == n-1
 801         False
 802
 803     The calculation of the Lyapunov rank of an improper cone can be
 804     reduced to that of a proper cone [Orlitzky/Gowda]_::
 805
 806         sage: set_random_seed()
 807         sage: K = random_cone(max_dim=8)
 808         sage: actual = lyapunov_rank(K)
 809         sage: K_S = restrict_span(K)
 810         sage: P = restrict_span(K_S.dual()).dual()
 811         sage: l = lineality(K)
 812         sage: c = codim(K)
 813         sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
 814         sage: actual == expected
 815         True
 816
 817     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 818
 819         sage: set_random_seed()
 820         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 821         sage: lyapunov_rank(K) == len(LL(K))
 822         True
 823
 824     In fact the same can be said of any cone. These additional tests
 825     just increase our confidence that the reduction scheme works::
 826
 827         sage: set_random_seed()
 828         sage: K = random_cone(max_dim=8, strictly_convex=True, 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_dim=8, strictly_convex=False, solid=True)
 836         sage: lyapunov_rank(K) == len(LL(K))
 837         True
 838
 839     ::
 840
 841         sage: set_random_seed()
 842         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 843         sage: lyapunov_rank(K) == len(LL(K))
 844         True
 845
 846     """
 847     K_orig = K
 848     beta = 0
 849
 850     m = K.dim()
 851     n = K.lattice_dim()
 852     l = lineality(K)
 853
 854     if m < n:
 855         # K is not solid, project onto its span.
 856         K = restrict_span(K)
 857
 858         # Lemma 2
 859         beta += m*(n - m) + (n - m)**2
 860
 861     if l > 0:
 862         # K is not pointed, project its dual onto its span.
 863         # Uses a proposition from our paper, i.e. this is
 864         # equivalent to K = restrict_span(K.dual()).dual()
 865         #K = restrict_span(intersect_span(K,K.dual()), K.dual())
 866         K = restrict_span(K, K.dual())
 867
 868         # Lemma 3
 869         beta += m * l
 870
 871     beta += len(LL(K))
 872     return beta