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(K1.lyapunov_like_basis()) != len(K2.lyapunov_like_basis()):
  80         return False
  81
  82     C_of_K1 = K1.discrete_complementarity_set()
  83     C_of_K2 = K2.discrete_complementarity_set()
  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 def lyapunov_rank(K):
 215     r"""
 216     Compute the Lyapunov rank of this cone.
 217
 218     The Lyapunov rank of a cone is the dimension of the space of its
 219     Lyapunov-like transformations -- that is, the length of a
 220     :meth:`lyapunov_like_basis`. Equivalently, the Lyapunov rank is the
 221     dimension of the Lie algebra of the automorphism group of the cone.
 222
 223     OUTPUT:
 224
 225     A nonnegative integer representing the Lyapunov rank of this cone.
 226
 227     If the ambient space is trivial, the Lyapunov rank will be zero.
 228     Otherwise, if the dimension of the ambient vector space is `n`, then
 229     the resulting Lyapunov rank will be between `1` and `n` inclusive. A
 230     Lyapunov rank of `n-1` is not possible [Orlitzky]_.
 231
 232     ALGORITHM:
 233
 234     The codimension formula from the second reference is used. We find
 235     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 236     and `s` are rays of our cone. It is known that these vectors are
 237     sufficient to apply the codimension formula. Once we have all such
 238     pairs, we "brute force" the codimension formula by finding all
 239     linearly-independent `xs^{T}`.
 240
 241     REFERENCES:
 242
 243     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of
 244        a proper cone and Lyapunov-like transformations. Mathematical
 245        Programming, 147 (2014) 155-170.
 246
 247     M. Orlitzky. The Lyapunov rank of an improper cone.
 248     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 249
 250     G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 251     optimality constraints for the cone of positive polynomials,
 252     Mathematical Programming, Series B, 129 (2011) 5-31.
 253
 254     EXAMPLES:
 255
 256     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 257     [Rudolf]_::
 258
 259         sage: positives = Cone([(1,)])
 260         sage: lyapunov_rank(positives)
 261         1
 262         sage: quadrant = Cone([(1,0), (0,1)])
 263         sage: lyapunov_rank(quadrant)
 264         2
 265        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 266         sage: lyapunov_rank(octant)
 267         3
 268
 269     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 270     [Orlitzky]_::
 271
 272         sage: R5 = VectorSpace(QQ, 5)
 273         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 274         sage: K = Cone(gs)
 275         sage: lyapunov_rank(K)
 276         25
 277
 278     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 279     [Rudolf]_::
 280
 281         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 282         sage: lyapunov_rank(L31)
 283         1
 284
 285     Likewise for the `L^{3}_{\infty}` cone [Rudolf]_::
 286
 287         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 288         sage: lyapunov_rank(L3infty)
 289         1
 290
 291     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 292     + 1` [Orlitzky]_::
 293
 294         sage: K = Cone([(1,0,0,0,0)])
 295         sage: lyapunov_rank(K)
 296         21
 297         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 298         21
 299
 300     A subspace (of dimension `m`) in `n` dimensions should have a
 301     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky]_::
 302
 303         sage: e1 = (1,0,0,0,0)
 304         sage: neg_e1 = (-1,0,0,0,0)
 305         sage: e2 = (0,1,0,0,0)
 306         sage: neg_e2 = (0,-1,0,0,0)
 307         sage: z = (0,0,0,0,0)
 308         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 309         sage: lyapunov_rank(K)
 310         19
 311         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 312         19
 313
 314     The Lyapunov rank should be additive on a product of proper cones
 315     [Rudolf]_::
 316
 317         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 318         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 319         sage: K = L31.cartesian_product(octant)
 320         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 321         True
 322
 323     Two isomorphic cones should have the same Lyapunov rank [Rudolf]_.
 324     The cone ``K`` in the following example is isomorphic to the nonnegative
 325     octant in `\mathbb{R}^{3}`::
 326
 327         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 328         sage: lyapunov_rank(K)
 329         3
 330
 331     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 332     itself [Rudolf]_::
 333
 334         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 335         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 336         True
 337
 338     TESTS:
 339
 340     The Lyapunov rank should be additive on a product of proper cones
 341     [Rudolf]_::
 342
 343         sage: set_random_seed()
 344         sage: K1 = random_cone(max_ambient_dim=8,
 345         ....:                  strictly_convex=True,
 346         ....:                  solid=True)
 347         sage: K2 = random_cone(max_ambient_dim=8,
 348         ....:                  strictly_convex=True,
 349         ....:                  solid=True)
 350         sage: K = K1.cartesian_product(K2)
 351         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 352         True
 353
 354     The Lyapunov rank is invariant under a linear isomorphism
 355     [Orlitzky]_::
 356
 357         sage: K1 = random_cone(max_ambient_dim = 8)
 358         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 359         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 360         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 361         True
 362
 363     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 364     itself [Rudolf]_::
 365
 366         sage: set_random_seed()
 367         sage: K = random_cone(max_ambient_dim=8)
 368         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 369         True
 370
 371     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 372     be any number between `1` and `n` inclusive, excluding `n-1`
 373     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 374     trivial cone in a trivial space as well. However, in zero dimensions,
 375     the Lyapunov rank of the trivial cone will be zero::
 376
 377         sage: set_random_seed()
 378         sage: K = random_cone(max_ambient_dim=8,
 379         ....:                 strictly_convex=True,
 380         ....:                 solid=True)
 381         sage: b = lyapunov_rank(K)
 382         sage: n = K.lattice_dim()
 383         sage: (n == 0 or 1 <= b) and b <= n
 384         True
 385         sage: b == n-1
 386         False
 387
 388     In fact [Orlitzky]_, no closed convex polyhedral cone can have
 389     Lyapunov rank `n-1` in `n` dimensions::
 390
 391         sage: set_random_seed()
 392         sage: K = random_cone(max_ambient_dim=8)
 393         sage: b = lyapunov_rank(K)
 394         sage: n = K.lattice_dim()
 395         sage: b == n-1
 396         False
 397
 398     The calculation of the Lyapunov rank of an improper cone can be
 399     reduced to that of a proper cone [Orlitzky]_::
 400
 401         sage: set_random_seed()
 402         sage: K = random_cone(max_ambient_dim=8)
 403         sage: actual = lyapunov_rank(K)
 404         sage: K_S = _restrict_to_space(K, K.span())
 405         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 406         sage: l = K.lineality()
 407         sage: c = K.codim()
 408         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 409         sage: actual == expected
 410         True
 411
 412     The Lyapunov rank of a cone is the size of a :meth:`lyapunov_like_basis`::
 413
 414         sage: set_random_seed()
 415         sage: K = random_cone(max_ambient_dim=8)
 416         sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
 417         True
 418
 419     We can make an imperfect cone perfect by adding a slack variable
 420     (a Theorem in [Orlitzky]_)::
 421
 422         sage: set_random_seed()
 423         sage: K = random_cone(max_ambient_dim=8,
 424         ....:                 strictly_convex=True,
 425         ....:                 solid=True)
 426         sage: L = ToricLattice(K.lattice_dim() + 1)
 427         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 428         sage: lyapunov_rank(K) >= K.lattice_dim()
 429         True
 430
 431     """
 432     beta = 0 # running tally of the Lyapunov rank
 433
 434     m = K.dim()
 435     n = K.lattice_dim()
 436     l = K.lineality()
 437
 438     if m < n:
 439         # K is not solid, restrict to its span.
 440         K = _restrict_to_space(K, K.span())
 441
 442         # Non-solid reduction lemma.
 443         beta += (n - m)*n
 444
 445     if l > 0:
 446         # K is not pointed, restrict to the span of its dual. Uses a
 447         # proposition from our paper, i.e. this is equivalent to K =
 448         # _rho(K.dual()).dual().
 449         K = _restrict_to_space(K, K.dual().span())
 450
 451         # Non-pointed reduction lemma.
 452         beta += l * m
 453
 454     beta += len(K.lyapunov_like_basis())
 455     return beta
 456
 457
 458
 459 def is_lyapunov_like(L,K):
 460     r"""
 461     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 462
 463     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 464     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 465     `\left\langle x,s \right\rangle` in the complementarity set of
 466     ``K``. It is known [Orlitzky]_ that this property need only be
 467     checked for generators of ``K`` and its dual.
 468
 469     INPUT:
 470
 471     - ``L`` -- A linear transformation or matrix.
 472
 473     - ``K`` -- A polyhedral closed convex cone.
 474
 475     OUTPUT:
 476
 477     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 478     and ``False`` otherwise.
 479
 480     .. WARNING::
 481
 482         If this function returns ``True``, then ``L`` is Lyapunov-like
 483         on ``K``. However, if ``False`` is returned, that could mean one
 484         of two things. The first is that ``L`` is definitely not
 485         Lyapunov-like on ``K``. The second is more of an "I don't know"
 486         answer, returned (for example) if we cannot prove that an inner
 487         product is zero.
 488
 489     REFERENCES:
 490
 491     M. Orlitzky. The Lyapunov rank of an improper cone.
 492     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 493
 494     EXAMPLES:
 495
 496     The identity is always Lyapunov-like in a nontrivial space::
 497
 498         sage: set_random_seed()
 499         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
 500         sage: L = identity_matrix(K.lattice_dim())
 501         sage: is_lyapunov_like(L,K)
 502         True
 503
 504     As is the "zero" transformation::
 505
 506         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 507         sage: R = K.lattice().vector_space().base_ring()
 508         sage: L = zero_matrix(R, K.lattice_dim())
 509         sage: is_lyapunov_like(L,K)
 510         True
 511
 512         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 513         on ``K``::
 514
 515         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 516         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
 517         True
 518
 519     """
 520     return all([(L*x).inner_product(s) == 0
 521                 for (x,s) in K.discrete_complementarity_set()])
 522
 523
 524 def random_element(K):
 525     r"""
 526     Return a random element of ``K`` from its ambient vector space.
 527
 528     ALGORITHM:
 529
 530     The cone ``K`` is specified in terms of its generators, so that
 531     ``K`` is equal to the convex conic combination of those generators.
 532     To choose a random element of ``K``, we assign random nonnegative
 533     coefficients to each generator of ``K`` and construct a new vector
 534     from the scaled rays.
 535
 536     A vector, rather than a ray, is returned so that the element may
 537     have non-integer coordinates. Thus the element may have an
 538     arbitrarily small norm.
 539
 540     EXAMPLES:
 541
 542     A random element of the trivial cone is zero::
 543
 544         sage: set_random_seed()
 545         sage: K = Cone([], ToricLattice(0))
 546         sage: random_element(K)
 547         ()
 548         sage: K = Cone([(0,)])
 549         sage: random_element(K)
 550         (0)
 551         sage: K = Cone([(0,0)])
 552         sage: random_element(K)
 553         (0, 0)
 554         sage: K = Cone([(0,0,0)])
 555         sage: random_element(K)
 556         (0, 0, 0)
 557
 558     TESTS:
 559
 560     Any cone should contain an element of itself::
 561
 562         sage: set_random_seed()
 563         sage: K = random_cone(max_rays = 8)
 564         sage: K.contains(random_element(K))
 565         True
 566
 567     """
 568     V = K.lattice().vector_space()
 569     F = V.base_ring()
 570     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 571     vector_gens  = map(V, K.rays())
 572     scaled_gens  = [ coefficients[i]*vector_gens[i]
 573                          for i in range(len(vector_gens)) ]
 574
 575     # Make sure we return a vector. Without the coercion, we might
 576     # return ``0`` when ``K`` has no rays.
 577     v = V(sum(scaled_gens))
 578     return v
 579
 580
 581 def positive_operators(K):
 582     r"""
 583     Compute generators of the cone of positive operators on this cone.
 584
 585     OUTPUT:
 586
 587     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 588     Each matrix ``P`` in the list should have the property that ``P*x``
 589     is an element of ``K`` whenever ``x`` is an element of
 590     ``K``. Moreover, any nonnegative linear combination of these
 591     matrices shares the same property.
 592
 593     EXAMPLES:
 594
 595     The trivial cone in a trivial space has no positive operators::
 596
 597         sage: K = Cone([], ToricLattice(0))
 598         sage: positive_operators(K)
 599         []
 600
 601     Positive operators on the nonnegative orthant are nonnegative matrices::
 602
 603         sage: K = Cone([(1,)])
 604         sage: positive_operators(K)
 605         [[1]]
 606
 607         sage: K = Cone([(1,0),(0,1)])
 608         sage: positive_operators(K)
 609         [
 610         [1 0]  [0 1]  [0 0]  [0 0]
 611         [0 0], [0 0], [1 0], [0 1]
 612         ]
 613
 614     Every operator is positive on the ambient vector space::
 615
 616         sage: K = Cone([(1,),(-1,)])
 617         sage: K.is_full_space()
 618         True
 619         sage: positive_operators(K)
 620         [[1], [-1]]
 621
 622         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 623         sage: K.is_full_space()
 624         True
 625         sage: positive_operators(K)
 626         [
 627         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 628         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 629         ]
 630
 631     TESTS:
 632
 633     A positive operator on a cone should send its generators into the cone::
 634
 635         sage: K = random_cone(max_ambient_dim = 6)
 636         sage: pi_of_K = positive_operators(K)
 637         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 638         True
 639
 640     """
 641     # Sage doesn't think matrices are vectors, so we have to convert
 642     # our matrices to vectors explicitly before we can figure out how
 643     # many are linearly-indepenedent.
 644     #
 645     # The space W has the same base ring as V, but dimension
 646     # dim(V)^2. So it has the same dimension as the space of linear
 647     # transformations on V. In other words, it's just the right size
 648     # to create an isomorphism between it and our matrices.
 649     V = K.lattice().vector_space()
 650     W = VectorSpace(V.base_ring(), V.dimension()**2)
 651
 652     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 653
 654     # Turn our matrices into long vectors...
 655     vectors = [ W(m.list()) for m in tensor_products ]
 656
 657     # Create the *dual* cone of the positive operators, expressed as
 658     # long vectors..
 659     L = ToricLattice(W.dimension())
 660     pi_dual = Cone(vectors, lattice=L)
 661
 662     # Now compute the desired cone from its dual...
 663     pi_cone = pi_dual.dual()
 664
 665     # And finally convert its rays back to matrix representations.
 666     M = MatrixSpace(V.base_ring(), V.dimension())
 667
 668     return [ M(v.list()) for v in pi_cone.rays() ]
 669
 670
 671 def Z_transformations(K):
 672     r"""
 673     Compute generators of the cone of Z-transformations on this cone.
 674
 675     OUTPUT:
 676
 677     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 678     Each matrix ``L`` in the list should have the property that
 679     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 680     discrete complementarity set of ``K``. Moreover, any nonnegative
 681     linear combination of these matrices shares the same property.
 682
 683     EXAMPLES:
 684
 685     Z-transformations on the nonnegative orthant are just Z-matrices.
 686     That is, matrices whose off-diagonal elements are nonnegative::
 687
 688         sage: K = Cone([(1,0),(0,1)])
 689         sage: Z_transformations(K)
 690         [
 691         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 692         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 693         ]
 694         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 695         sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
 696         ....:                    for i in range(z.nrows())
 697         ....:                    for j in range(z.ncols())
 698         ....:                    if i != j ])
 699         True
 700
 701     The trivial cone in a trivial space has no Z-transformations::
 702
 703         sage: K = Cone([], ToricLattice(0))
 704         sage: Z_transformations(K)
 705         []
 706
 707     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 708
 709         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 710         sage: K.is_full_space()
 711         True
 712         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 713         sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
 714         sage: zs == lls
 715         True
 716
 717     TESTS:
 718
 719     The Z-property is possessed by every Z-transformation::
 720
 721         sage: set_random_seed()
 722         sage: K = random_cone(max_ambient_dim = 6)
 723         sage: Z_of_K = Z_transformations(K)
 724         sage: dcs = K.discrete_complementarity_set()
 725         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 726         ....:                                  for (x,s) in dcs])
 727         True
 728
 729     The lineality space of Z is LL::
 730
 731         sage: set_random_seed()
 732         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 733         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 734         sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
 735         sage: z_cone.linear_subspace() == lls
 736         True
 737
 738     """
 739     # Sage doesn't think matrices are vectors, so we have to convert
 740     # our matrices to vectors explicitly before we can figure out how
 741     # many are linearly-indepenedent.
 742     #
 743     # The space W has the same base ring as V, but dimension
 744     # dim(V)^2. So it has the same dimension as the space of linear
 745     # transformations on V. In other words, it's just the right size
 746     # to create an isomorphism between it and our matrices.
 747     V = K.lattice().vector_space()
 748     W = VectorSpace(V.base_ring(), V.dimension()**2)
 749
 750     C_of_K = K.discrete_complementarity_set()
 751     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 752
 753     # Turn our matrices into long vectors...
 754     vectors = [ W(m.list()) for m in tensor_products ]
 755
 756     # Create the *dual* cone of the cross-positive operators,
 757     # expressed as long vectors..
 758     L = ToricLattice(W.dimension())
 759     Sigma_dual = Cone(vectors, lattice=L)
 760
 761     # Now compute the desired cone from its dual...
 762     Sigma_cone = Sigma_dual.dual()
 763
 764     # And finally convert its rays back to matrix representations.
 765     # But first, make them negative, so we get Z-transformations and
 766     # not cross-positive ones.
 767     M = MatrixSpace(V.base_ring(), V.dimension())
 768
 769     return [ -M(v.list()) for v in Sigma_cone.rays() ]