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