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