mjo/cone/cone.py

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