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 codim(K):
 412     r"""
 413     Compute the codimension of this cone.
 414
 415     The codimension of a cone is the dimension of the space of all
 416     elements perpendicular to every element of the cone. In other words,
 417     the codimension is the difference between the dimension of the
 418     ambient space and the dimension of the cone itself.
 419
 420     OUTPUT:
 421
 422     A nonnegative integer representing the dimension of the space of all
 423     elements perpendicular to this cone.
 424
 425     .. seealso::
 426
 427         :meth:`dim`, :meth:`lattice_dim`
 428
 429     EXAMPLES:
 430
 431     The codimension of the nonnegative orthant is zero, since the span of
 432     its generators equals the entire ambient space::
 433
 434         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
 435         sage: codim(K)
 436         0
 437
 438     However, if we remove a ray so that the entire cone is contained
 439     within the `x-y`-plane, then the resulting cone will have
 440     codimension one, because the `z`-axis is perpendicular to every
 441     element of the cone::
 442
 443         sage: K = Cone([(1,0,0), (0,1,0)])
 444         sage: codim(K)
 445         1
 446
 447     If our cone is all of `\mathbb{R}^{2}`, then its codimension is zero::
 448
 449         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 450         sage: codim(K)
 451         0
 452
 453     And if the cone is trivial in any space, then its codimension is
 454     equal to the dimension of the ambient space::
 455
 456         sage: K = Cone([], lattice=ToricLattice(0))
 457         sage: K.lattice_dim()
 458         0
 459         sage: codim(K)
 460         0
 461
 462         sage: K = Cone([(0,)])
 463         sage: K.lattice_dim()
 464         1
 465         sage: codim(K)
 466         1
 467
 468         sage: K = Cone([(0,0)])
 469         sage: K.lattice_dim()
 470         2
 471         sage: codim(K)
 472         2
 473
 474     TESTS:
 475
 476     The codimension of a cone should be an integer between zero and
 477     the dimension of the ambient space, inclusive::
 478
 479         sage: set_random_seed()
 480         sage: K = random_cone(max_dim = 8)
 481         sage: c = codim(K)
 482         sage: c in ZZ
 483         True
 484         sage: (0 <= c) and (c <= K.lattice_dim())
 485         True
 486
 487     A solid cone should have codimension zero::
 488
 489         sage: set_random_seed()
 490         sage: K = random_cone(max_dim = 8, solid = True)
 491         sage: codim(K)
 492         0
 493
 494     The codimension of a cone is equal to the lineality of its dual::
 495
 496         sage: set_random_seed()
 497         sage: K = random_cone(max_dim = 8, solid = True)
 498         sage: codim(K) == lineality(K.dual())
 499         True
 500
 501     """
 502     return (K.lattice_dim() - K.dim())
 503
 504
 505 def discrete_complementarity_set(K):
 506     r"""
 507     Compute the discrete complementarity set of this cone.
 508
 509     The complementarity set of this cone is the set of all orthogonal
 510     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 511     dual. The discrete complementarity set restricts `x` and `s` to be
 512     generators of their respective cones.
 513
 514     OUTPUT:
 515
 516     A list of pairs `(x,s)` such that,
 517
 518       * `x` is in this cone.
 519       * `x` is a generator of this cone.
 520       * `s` is in this cone's dual.
 521       * `s` is a generator of this cone's dual.
 522       * `x` and `s` are orthogonal.
 523
 524     EXAMPLES:
 525
 526     The discrete complementarity set of the nonnegative orthant consists
 527     of pairs of standard basis vectors::
 528
 529         sage: K = Cone([(1,0),(0,1)])
 530         sage: discrete_complementarity_set(K)
 531         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 532
 533     If the cone consists of a single ray, the second components of the
 534     discrete complementarity set should generate the orthogonal
 535     complement of that ray::
 536
 537         sage: K = Cone([(1,0)])
 538         sage: discrete_complementarity_set(K)
 539         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 540         sage: K = Cone([(1,0,0)])
 541         sage: discrete_complementarity_set(K)
 542         [((1, 0, 0), (0, 1, 0)),
 543           ((1, 0, 0), (0, -1, 0)),
 544           ((1, 0, 0), (0, 0, 1)),
 545           ((1, 0, 0), (0, 0, -1))]
 546
 547     When the cone is the entire space, its dual is the trivial cone, so
 548     the discrete complementarity set is empty::
 549
 550         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 551         sage: discrete_complementarity_set(K)
 552         []
 553
 554     TESTS:
 555
 556     The complementarity set of the dual can be obtained by switching the
 557     components of the complementarity set of the original cone::
 558
 559         sage: set_random_seed()
 560         sage: K1 = random_cone(max_dim=6)
 561         sage: K2 = K1.dual()
 562         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 563         sage: actual = discrete_complementarity_set(K1)
 564         sage: sorted(actual) == sorted(expected)
 565         True
 566
 567     """
 568     V = K.lattice().vector_space()
 569
 570     # Convert the rays to vectors so that we can compute inner
 571     # products.
 572     xs = [V(x) for x in K.rays()]
 573     ss = [V(s) for s in K.dual().rays()]
 574
 575     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 576
 577
 578 def LL(K):
 579     r"""
 580     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 581     on this cone.
 582
 583     OUTPUT:
 584
 585     A list of matrices forming a basis for the space of all
 586     Lyapunov-like transformations on the given cone.
 587
 588     EXAMPLES:
 589
 590     The trivial cone has no Lyapunov-like transformations::
 591
 592         sage: L = ToricLattice(0)
 593         sage: K = Cone([], lattice=L)
 594         sage: LL(K)
 595         []
 596
 597     The Lyapunov-like transformations on the nonnegative orthant are
 598     simply diagonal matrices::
 599
 600         sage: K = Cone([(1,)])
 601         sage: LL(K)
 602         [[1]]
 603
 604         sage: K = Cone([(1,0),(0,1)])
 605         sage: LL(K)
 606         [
 607         [1 0]  [0 0]
 608         [0 0], [0 1]
 609         ]
 610
 611         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 612         sage: LL(K)
 613         [
 614         [1 0 0]  [0 0 0]  [0 0 0]
 615         [0 0 0]  [0 1 0]  [0 0 0]
 616         [0 0 0], [0 0 0], [0 0 1]
 617         ]
 618
 619     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 620     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 621
 622         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 623         sage: LL(L31)
 624         [
 625         [1 0 0]
 626         [0 1 0]
 627         [0 0 1]
 628         ]
 629
 630         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 631         sage: LL(L3infty)
 632         [
 633         [1 0 0]
 634         [0 1 0]
 635         [0 0 1]
 636         ]
 637
 638     If our cone is the entire space, then every transformation on it is
 639     Lyapunov-like::
 640
 641         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 642         sage: M = MatrixSpace(QQ,2)
 643         sage: M.basis() == LL(K)
 644         True
 645
 646     TESTS:
 647
 648     The inner product `\left< L\left(x\right), s \right>` is zero for
 649     every pair `\left( x,s \right)` in the discrete complementarity set
 650     of the cone::
 651
 652         sage: set_random_seed()
 653         sage: K = random_cone(max_dim=8)
 654         sage: C_of_K = discrete_complementarity_set(K)
 655         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 656         sage: sum(map(abs, l))
 657         0
 658
 659     The Lyapunov-like transformations on a cone and its dual are related
 660     by transposition, but we're not guaranteed to compute transposed
 661     elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
 662     \right)`
 663
 664         sage: set_random_seed()
 665         sage: K = random_cone(max_dim=8)
 666         sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
 667         sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
 668         sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
 669         sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
 670         sage: V.span(LL1_vecs) == V.span(LL2_vecs)
 671         True
 672
 673     """
 674     V = K.lattice().vector_space()
 675
 676     C_of_K = discrete_complementarity_set(K)
 677
 678     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 679
 680     # Sage doesn't think matrices are vectors, so we have to convert
 681     # our matrices to vectors explicitly before we can figure out how
 682     # many are linearly-indepenedent.
 683     #
 684     # The space W has the same base ring as V, but dimension
 685     # dim(V)^2. So it has the same dimension as the space of linear
 686     # transformations on V. In other words, it's just the right size
 687     # to create an isomorphism between it and our matrices.
 688     W = VectorSpace(V.base_ring(), V.dimension()**2)
 689
 690     # Turn our matrices into long vectors...
 691     vectors = [ W(m.list()) for m in tensor_products ]
 692
 693     # Vector space representation of Lyapunov-like matrices
 694     # (i.e. vec(L) where L is Luapunov-like).
 695     LL_vector = W.span(vectors).complement()
 696
 697     # Now construct an ambient MatrixSpace in which to stick our
 698     # transformations.
 699     M = MatrixSpace(V.base_ring(), V.dimension())
 700
 701     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 702
 703     return matrix_basis
 704
 705
 706
 707 def lyapunov_rank(K):
 708     r"""
 709     Compute the Lyapunov (or bilinearity) rank of this cone.
 710
 711     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 712
 713     1. The dimension of the Lie algebra of the automorphism group of the
 714        cone.
 715
 716     2. The dimension of the linear space of all Lyapunov-like
 717        transformations on the cone.
 718
 719     INPUT:
 720
 721     A closed, convex polyhedral cone.
 722
 723     OUTPUT:
 724
 725     An integer representing the Lyapunov rank of the cone. If the
 726     dimension of the ambient vector space is `n`, then the Lyapunov rank
 727     will be between `1` and `n` inclusive; however a rank of `n-1` is
 728     not possible (see the first reference).
 729
 730     .. note::
 731
 732         In the references, the cones are always assumed to be proper. We
 733         do not impose this restriction.
 734
 735     .. seealso::
 736
 737         :meth:`is_proper`
 738
 739     ALGORITHM:
 740
 741     The codimension formula from the second reference is used. We find
 742     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 743     and `s` are rays of our cone. It is known that these vectors are
 744     sufficient to apply the codimension formula. Once we have all such
 745     pairs, we "brute force" the codimension formula by finding all
 746     linearly-independent `xs^{T}`.
 747
 748     REFERENCES:
 749
 750     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 751        cone and Lyapunov-like transformations, Mathematical Programming, 147
 752        (2014) 155-170.
 753
 754     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 755        Improper Cone. Work in-progress.
 756
 757     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 758        optimality constraints for the cone of positive polynomials,
 759        Mathematical Programming, Series B, 129 (2011) 5-31.
 760
 761     EXAMPLES:
 762
 763     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 764     [Rudolf et al.]_::
 765
 766         sage: positives = Cone([(1,)])
 767         sage: lyapunov_rank(positives)
 768         1
 769         sage: quadrant = Cone([(1,0), (0,1)])
 770         sage: lyapunov_rank(quadrant)
 771         2
 772        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 773         sage: lyapunov_rank(octant)
 774         3
 775
 776     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 777     [Orlitzky/Gowda]_::
 778
 779         sage: R5 = VectorSpace(QQ, 5)
 780         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 781         sage: K = Cone(gs)
 782         sage: lyapunov_rank(K)
 783         25
 784
 785     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 786     [Rudolf et al.]_::
 787
 788         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 789         sage: lyapunov_rank(L31)
 790         1
 791
 792     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 793
 794         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 795         sage: lyapunov_rank(L3infty)
 796         1
 797
 798     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 799     + 1` [Orlitzky/Gowda]_::
 800
 801         sage: K = Cone([(1,0,0,0,0)])
 802         sage: lyapunov_rank(K)
 803         21
 804         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 805         21
 806
 807     A subspace (of dimension `m`) in `n` dimensions should have a
 808     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 809
 810         sage: e1 = (1,0,0,0,0)
 811         sage: neg_e1 = (-1,0,0,0,0)
 812         sage: e2 = (0,1,0,0,0)
 813         sage: neg_e2 = (0,-1,0,0,0)
 814         sage: z = (0,0,0,0,0)
 815         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 816         sage: lyapunov_rank(K)
 817         19
 818         sage: K.lattice_dim()**2 - K.dim()*codim(K)
 819         19
 820
 821     The Lyapunov rank should be additive on a product of proper cones
 822     [Rudolf et al.]_::
 823
 824         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 825         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 826         sage: K = L31.cartesian_product(octant)
 827         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 828         True
 829
 830     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 831     The cone ``K`` in the following example is isomorphic to the nonnegative
 832     octant in `\mathbb{R}^{3}`::
 833
 834         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 835         sage: lyapunov_rank(K)
 836         3
 837
 838     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 839     itself [Rudolf et al.]_::
 840
 841         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 842         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 843         True
 844
 845     TESTS:
 846
 847     The Lyapunov rank should be additive on a product of proper cones
 848     [Rudolf et al.]_::
 849
 850         sage: set_random_seed()
 851         sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 852         sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
 853         sage: K = K1.cartesian_product(K2)
 854         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 855         True
 856
 857     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 858     itself [Rudolf et al.]_::
 859
 860         sage: set_random_seed()
 861         sage: K = random_cone(max_dim=8)
 862         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 863         True
 864
 865     Make sure we exercise the non-strictly-convex/non-solid case::
 866
 867         sage: set_random_seed()
 868         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 869         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 870         True
 871
 872     Let's check the other permutations as well, just to be sure::
 873
 874         sage: set_random_seed()
 875         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 876         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 877         True
 878
 879     ::
 880
 881         sage: set_random_seed()
 882         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 883         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 884         True
 885
 886     ::
 887
 888         sage: set_random_seed()
 889         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 890         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 891         True
 892
 893     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 894     be any number between `1` and `n` inclusive, excluding `n-1`
 895     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 896     trivial cone in a trivial space as well. However, in zero dimensions,
 897     the Lyapunov rank of the trivial cone will be zero::
 898
 899         sage: set_random_seed()
 900         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 901         sage: b = lyapunov_rank(K)
 902         sage: n = K.lattice_dim()
 903         sage: (n == 0 or 1 <= b) and b <= n
 904         True
 905         sage: b == n-1
 906         False
 907
 908     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 909     Lyapunov rank `n-1` in `n` dimensions::
 910
 911         sage: set_random_seed()
 912         sage: K = random_cone(max_dim=8)
 913         sage: b = lyapunov_rank(K)
 914         sage: n = K.lattice_dim()
 915         sage: b == n-1
 916         False
 917
 918     The calculation of the Lyapunov rank of an improper cone can be
 919     reduced to that of a proper cone [Orlitzky/Gowda]_::
 920
 921         sage: set_random_seed()
 922         sage: K = random_cone(max_dim=8)
 923         sage: actual = lyapunov_rank(K)
 924         sage: K_S = restrict_span(K)
 925         sage: P = restrict_span(K_S.dual()).dual()
 926         sage: l = lineality(K)
 927         sage: c = codim(K)
 928         sage: expected = lyapunov_rank(P) + K.dim()*(l + c) + c**2
 929         sage: actual == expected
 930         True
 931
 932     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 933
 934         sage: set_random_seed()
 935         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
 936         sage: lyapunov_rank(K) == len(LL(K))
 937         True
 938
 939     In fact the same can be said of any cone. These additional tests
 940     just increase our confidence that the reduction scheme works::
 941
 942         sage: set_random_seed()
 943         sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
 944         sage: lyapunov_rank(K) == len(LL(K))
 945         True
 946
 947     ::
 948
 949         sage: set_random_seed()
 950         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
 951         sage: lyapunov_rank(K) == len(LL(K))
 952         True
 953
 954     ::
 955
 956         sage: set_random_seed()
 957         sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
 958         sage: lyapunov_rank(K) == len(LL(K))
 959         True
 960
 961     """
 962     K_orig = K
 963     beta = 0
 964
 965     m = K.dim()
 966     n = K.lattice_dim()
 967     l = lineality(K)
 968
 969     if m < n:
 970         # K is not solid, project onto its span.
 971         K = restrict_span(K)
 972
 973         # Lemma 2
 974         beta += m*(n - m) + (n - m)**2
 975
 976     if l > 0:
 977         # K is not pointed, project its dual onto its span.
 978         # Uses a proposition from our paper, i.e. this is
 979         # equivalent to K = restrict_span(K.dual()).dual()
 980         #K = restrict_span(intersect_span(K,K.dual()), K.dual())
 981         K = restrict_span(K, K.dual())
 982
 983         # Lemma 3
 984         beta += m * l
 985
 986     beta += len(LL(K))
 987     return beta