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