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