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 _restrict_to_space(K, W):
  12     r"""
  13     Restrict this cone a subspace of its ambient space.
  14
  15     INPUT:
  16
  17     - ``W`` -- The subspace into which this cone will be restricted.
  18
  19     OUTPUT:
  20
  21     A new cone in a sublattice corresponding to ``W``.
  22
  23     EXAMPLES:
  24
  25     When this cone is solid, restricting it into its own span should do
  26     nothing::
  27
  28         sage: K = Cone([(1,)])
  29         sage: _restrict_to_space(K, K.span()) == K
  30         True
  31
  32     A single ray restricted into its own span gives the same output
  33     regardless of the ambient space::
  34
  35         sage: K2 = Cone([(1,0)])
  36         sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
  37         sage: K2_S
  38         N(1)
  39         in 1-d lattice N
  40         sage: K3 = Cone([(1,0,0)])
  41         sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
  42         sage: K3_S
  43         N(1)
  44         in 1-d lattice N
  45         sage: K2_S == K3_S
  46         True
  47
  48     TESTS:
  49
  50     The projected cone should always be solid::
  51
  52         sage: set_random_seed()
  53         sage: K = random_cone(max_ambient_dim = 8)
  54         sage: _restrict_to_space(K, K.span()).is_solid()
  55         True
  56
  57     And the resulting cone should live in a space having the same
  58     dimension as the space we restricted it to::
  59
  60         sage: set_random_seed()
  61         sage: K = random_cone(max_ambient_dim = 8)
  62         sage: K_P = _restrict_to_space(K, K.dual().span())
  63         sage: K_P.lattice_dim() == K.dual().dim()
  64         True
  65
  66     This function should not affect the dimension of a cone::
  67
  68         sage: set_random_seed()
  69         sage: K = random_cone(max_ambient_dim = 8)
  70         sage: K.dim() == _restrict_to_space(K,K.span()).dim()
  71         True
  72
  73     Nor should it affect the lineality of a cone::
  74
  75         sage: set_random_seed()
  76         sage: K = random_cone(max_ambient_dim = 8)
  77         sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
  78         True
  79
  80     No matter which space we restrict to, the lineality should not
  81     increase::
  82
  83         sage: set_random_seed()
  84         sage: K = random_cone(max_ambient_dim = 8)
  85         sage: S = K.span(); P = K.dual().span()
  86         sage: K.lineality() >= _restrict_to_space(K,S).lineality()
  87         True
  88         sage: K.lineality() >= _restrict_to_space(K,P).lineality()
  89         True
  90
  91     If we do this according to our paper, then the result is proper::
  92
  93         sage: set_random_seed()
  94         sage: K = random_cone(max_ambient_dim = 8)
  95         sage: K_S = _restrict_to_space(K, K.span())
  96         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
  97         sage: K_SP.is_proper()
  98         True
  99         sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
 100         sage: K_SP.is_proper()
 101         True
 102
 103     Test the proposition in our paper concerning the duals and
 104     restrictions. Generate a random cone, then create a subcone of
 105     it. The operation of dual-taking should then commute with
 106     _restrict_to_space::
 107
 108         sage: set_random_seed()
 109         sage: J = random_cone(max_ambient_dim = 8)
 110         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 111         sage: K_W_star = _restrict_to_space(K, J.span()).dual()
 112         sage: K_star_W = _restrict_to_space(K.dual(), J.span())
 113         sage: _basically_the_same(K_W_star, K_star_W)
 114         True
 115
 116     """
 117     # First we want to intersect ``K`` with ``W``. The easiest way to
 118     # do this is via cone intersection, so we turn the subspace ``W``
 119     # into a cone.
 120     W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice())
 121     K = K.intersection(W_cone)
 122
 123     # We've already intersected K with the span of K2, so every
 124     # generator of K should belong to W now.
 125     K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 126
 127     L = ToricLattice(W.dimension())
 128     return Cone(K_W_rays, lattice=L)
 129
 130
 131 def lyapunov_rank(K):
 132     r"""
 133     Compute the Lyapunov rank of this cone.
 134
 135     The Lyapunov rank of a cone is the dimension of the space of its
 136     Lyapunov-like transformations -- that is, the length of a
 137     :meth:`lyapunov_like_basis`. Equivalently, the Lyapunov rank is the
 138     dimension of the Lie algebra of the automorphism group of the cone.
 139
 140     OUTPUT:
 141
 142     A nonnegative integer representing the Lyapunov rank of this cone.
 143
 144     If the ambient space is trivial, the Lyapunov rank will be zero.
 145     Otherwise, if the dimension of the ambient vector space is `n`, then
 146     the resulting Lyapunov rank will be between `1` and `n` inclusive. A
 147     Lyapunov rank of `n-1` is not possible [Orlitzky]_.
 148
 149     ALGORITHM:
 150
 151     The codimension formula from the second reference is used. We find
 152     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 153     and `s` are rays of our cone. It is known that these vectors are
 154     sufficient to apply the codimension formula. Once we have all such
 155     pairs, we "brute force" the codimension formula by finding all
 156     linearly-independent `xs^{T}`.
 157
 158     REFERENCES:
 159
 160     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of
 161        a proper cone and Lyapunov-like transformations. Mathematical
 162        Programming, 147 (2014) 155-170.
 163
 164     M. Orlitzky. The Lyapunov rank of an improper cone.
 165     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 166
 167     G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 168     optimality constraints for the cone of positive polynomials,
 169     Mathematical Programming, Series B, 129 (2011) 5-31.
 170
 171     EXAMPLES:
 172
 173     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 174     [Rudolf]_::
 175
 176         sage: positives = Cone([(1,)])
 177         sage: lyapunov_rank(positives)
 178         1
 179         sage: quadrant = Cone([(1,0), (0,1)])
 180         sage: lyapunov_rank(quadrant)
 181         2
 182        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 183         sage: lyapunov_rank(octant)
 184         3
 185
 186     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 187     [Orlitzky]_::
 188
 189         sage: R5 = VectorSpace(QQ, 5)
 190         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 191         sage: K = Cone(gs)
 192         sage: lyapunov_rank(K)
 193         25
 194
 195     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 196     [Rudolf]_::
 197
 198         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 199         sage: lyapunov_rank(L31)
 200         1
 201
 202     Likewise for the `L^{3}_{\infty}` cone [Rudolf]_::
 203
 204         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 205         sage: lyapunov_rank(L3infty)
 206         1
 207
 208     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 209     + 1` [Orlitzky]_::
 210
 211         sage: K = Cone([(1,0,0,0,0)])
 212         sage: lyapunov_rank(K)
 213         21
 214         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 215         21
 216
 217     A subspace (of dimension `m`) in `n` dimensions should have a
 218     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky]_::
 219
 220         sage: e1 = (1,0,0,0,0)
 221         sage: neg_e1 = (-1,0,0,0,0)
 222         sage: e2 = (0,1,0,0,0)
 223         sage: neg_e2 = (0,-1,0,0,0)
 224         sage: z = (0,0,0,0,0)
 225         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 226         sage: lyapunov_rank(K)
 227         19
 228         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 229         19
 230
 231     The Lyapunov rank should be additive on a product of proper cones
 232     [Rudolf]_::
 233
 234         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 235         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 236         sage: K = L31.cartesian_product(octant)
 237         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 238         True
 239
 240     Two isomorphic cones should have the same Lyapunov rank [Rudolf]_.
 241     The cone ``K`` in the following example is isomorphic to the nonnegative
 242     octant in `\mathbb{R}^{3}`::
 243
 244         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 245         sage: lyapunov_rank(K)
 246         3
 247
 248     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 249     itself [Rudolf]_::
 250
 251         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 252         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 253         True
 254
 255     TESTS:
 256
 257     The Lyapunov rank should be additive on a product of proper cones
 258     [Rudolf]_::
 259
 260         sage: set_random_seed()
 261         sage: K1 = random_cone(max_ambient_dim=8,
 262         ....:                  strictly_convex=True,
 263         ....:                  solid=True)
 264         sage: K2 = random_cone(max_ambient_dim=8,
 265         ....:                  strictly_convex=True,
 266         ....:                  solid=True)
 267         sage: K = K1.cartesian_product(K2)
 268         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 269         True
 270
 271     The Lyapunov rank is invariant under a linear isomorphism
 272     [Orlitzky]_::
 273
 274         sage: K1 = random_cone(max_ambient_dim = 8)
 275         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 276         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 277         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 278         True
 279
 280     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 281     itself [Rudolf]_::
 282
 283         sage: set_random_seed()
 284         sage: K = random_cone(max_ambient_dim=8)
 285         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 286         True
 287
 288     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 289     be any number between `1` and `n` inclusive, excluding `n-1`
 290     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 291     trivial cone in a trivial space as well. However, in zero dimensions,
 292     the Lyapunov rank of the trivial cone will be zero::
 293
 294         sage: set_random_seed()
 295         sage: K = random_cone(max_ambient_dim=8,
 296         ....:                 strictly_convex=True,
 297         ....:                 solid=True)
 298         sage: b = lyapunov_rank(K)
 299         sage: n = K.lattice_dim()
 300         sage: (n == 0 or 1 <= b) and b <= n
 301         True
 302         sage: b == n-1
 303         False
 304
 305     In fact [Orlitzky]_, no closed convex polyhedral cone can have
 306     Lyapunov rank `n-1` in `n` dimensions::
 307
 308         sage: set_random_seed()
 309         sage: K = random_cone(max_ambient_dim=8)
 310         sage: b = lyapunov_rank(K)
 311         sage: n = K.lattice_dim()
 312         sage: b == n-1
 313         False
 314
 315     The calculation of the Lyapunov rank of an improper cone can be
 316     reduced to that of a proper cone [Orlitzky]_::
 317
 318         sage: set_random_seed()
 319         sage: K = random_cone(max_ambient_dim=8)
 320         sage: actual = lyapunov_rank(K)
 321         sage: K_S = _restrict_to_space(K, K.span())
 322         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 323         sage: l = K.lineality()
 324         sage: c = K.codim()
 325         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 326         sage: actual == expected
 327         True
 328
 329     The Lyapunov rank of a cone is the size of a :meth:`lyapunov_like_basis`::
 330
 331         sage: set_random_seed()
 332         sage: K = random_cone(max_ambient_dim=8)
 333         sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
 334         True
 335
 336     We can make an imperfect cone perfect by adding a slack variable
 337     (a Theorem in [Orlitzky]_)::
 338
 339         sage: set_random_seed()
 340         sage: K = random_cone(max_ambient_dim=8,
 341         ....:                 strictly_convex=True,
 342         ....:                 solid=True)
 343         sage: L = ToricLattice(K.lattice_dim() + 1)
 344         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 345         sage: lyapunov_rank(K) >= K.lattice_dim()
 346         True
 347
 348     """
 349     beta = 0 # running tally of the Lyapunov rank
 350
 351     m = K.dim()
 352     n = K.lattice_dim()
 353     l = K.lineality()
 354
 355     if m < n:
 356         # K is not solid, restrict to its span.
 357         K = _restrict_to_space(K, K.span())
 358
 359         # Non-solid reduction lemma.
 360         beta += (n - m)*n
 361
 362     if l > 0:
 363         # K is not pointed, restrict to the span of its dual. Uses a
 364         # proposition from our paper, i.e. this is equivalent to K =
 365         # _rho(K.dual()).dual().
 366         K = _restrict_to_space(K, K.dual().span())
 367
 368         # Non-pointed reduction lemma.
 369         beta += l * m
 370
 371     beta += len(K.lyapunov_like_basis())
 372     return beta
 373
 374
 375
 376 def is_lyapunov_like(L,K):
 377     r"""
 378     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 379
 380     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 381     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 382     `\left\langle x,s \right\rangle` in the complementarity set of
 383     ``K``. It is known [Orlitzky]_ that this property need only be
 384     checked for generators of ``K`` and its dual.
 385
 386     INPUT:
 387
 388     - ``L`` -- A linear transformation or matrix.
 389
 390     - ``K`` -- A polyhedral closed convex cone.
 391
 392     OUTPUT:
 393
 394     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 395     and ``False`` otherwise.
 396
 397     .. WARNING::
 398
 399         If this function returns ``True``, then ``L`` is Lyapunov-like
 400         on ``K``. However, if ``False`` is returned, that could mean one
 401         of two things. The first is that ``L`` is definitely not
 402         Lyapunov-like on ``K``. The second is more of an "I don't know"
 403         answer, returned (for example) if we cannot prove that an inner
 404         product is zero.
 405
 406     REFERENCES:
 407
 408     M. Orlitzky. The Lyapunov rank of an improper cone.
 409     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 410
 411     EXAMPLES:
 412
 413     The identity is always Lyapunov-like in a nontrivial space::
 414
 415         sage: set_random_seed()
 416         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
 417         sage: L = identity_matrix(K.lattice_dim())
 418         sage: is_lyapunov_like(L,K)
 419         True
 420
 421     As is the "zero" transformation::
 422
 423         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 424         sage: R = K.lattice().vector_space().base_ring()
 425         sage: L = zero_matrix(R, K.lattice_dim())
 426         sage: is_lyapunov_like(L,K)
 427         True
 428
 429         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 430         on ``K``::
 431
 432         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 433         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
 434         True
 435
 436     """
 437     return all([(L*x).inner_product(s) == 0
 438                 for (x,s) in K.discrete_complementarity_set()])
 439
 440
 441 def random_element(K):
 442     r"""
 443     Return a random element of ``K`` from its ambient vector space.
 444
 445     ALGORITHM:
 446
 447     The cone ``K`` is specified in terms of its generators, so that
 448     ``K`` is equal to the convex conic combination of those generators.
 449     To choose a random element of ``K``, we assign random nonnegative
 450     coefficients to each generator of ``K`` and construct a new vector
 451     from the scaled rays.
 452
 453     A vector, rather than a ray, is returned so that the element may
 454     have non-integer coordinates. Thus the element may have an
 455     arbitrarily small norm.
 456
 457     EXAMPLES:
 458
 459     A random element of the trivial cone is zero::
 460
 461         sage: set_random_seed()
 462         sage: K = Cone([], ToricLattice(0))
 463         sage: random_element(K)
 464         ()
 465         sage: K = Cone([(0,)])
 466         sage: random_element(K)
 467         (0)
 468         sage: K = Cone([(0,0)])
 469         sage: random_element(K)
 470         (0, 0)
 471         sage: K = Cone([(0,0,0)])
 472         sage: random_element(K)
 473         (0, 0, 0)
 474
 475     TESTS:
 476
 477     Any cone should contain an element of itself::
 478
 479         sage: set_random_seed()
 480         sage: K = random_cone(max_rays = 8)
 481         sage: K.contains(random_element(K))
 482         True
 483
 484     """
 485     V = K.lattice().vector_space()
 486     F = V.base_ring()
 487     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 488     vector_gens  = map(V, K.rays())
 489     scaled_gens  = [ coefficients[i]*vector_gens[i]
 490                          for i in range(len(vector_gens)) ]
 491
 492     # Make sure we return a vector. Without the coercion, we might
 493     # return ``0`` when ``K`` has no rays.
 494     v = V(sum(scaled_gens))
 495     return v
 496
 497
 498 def positive_operators(K):
 499     r"""
 500     Compute generators of the cone of positive operators on this cone.
 501
 502     OUTPUT:
 503
 504     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 505     Each matrix ``P`` in the list should have the property that ``P*x``
 506     is an element of ``K`` whenever ``x`` is an element of
 507     ``K``. Moreover, any nonnegative linear combination of these
 508     matrices shares the same property.
 509
 510     EXAMPLES:
 511
 512     The trivial cone in a trivial space has no positive operators::
 513
 514         sage: K = Cone([], ToricLattice(0))
 515         sage: positive_operators(K)
 516         []
 517
 518     Positive operators on the nonnegative orthant are nonnegative matrices::
 519
 520         sage: K = Cone([(1,)])
 521         sage: positive_operators(K)
 522         [[1]]
 523
 524         sage: K = Cone([(1,0),(0,1)])
 525         sage: positive_operators(K)
 526         [
 527         [1 0]  [0 1]  [0 0]  [0 0]
 528         [0 0], [0 0], [1 0], [0 1]
 529         ]
 530
 531     Every operator is positive on the ambient vector space::
 532
 533         sage: K = Cone([(1,),(-1,)])
 534         sage: K.is_full_space()
 535         True
 536         sage: positive_operators(K)
 537         [[1], [-1]]
 538
 539         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 540         sage: K.is_full_space()
 541         True
 542         sage: positive_operators(K)
 543         [
 544         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 545         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 546         ]
 547
 548     TESTS:
 549
 550     A positive operator on a cone should send its generators into the cone::
 551
 552         sage: K = random_cone(max_ambient_dim = 6)
 553         sage: pi_of_K = positive_operators(K)
 554         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 555         True
 556
 557     """
 558     # Sage doesn't think matrices are vectors, so we have to convert
 559     # our matrices to vectors explicitly before we can figure out how
 560     # many are linearly-indepenedent.
 561     #
 562     # The space W has the same base ring as V, but dimension
 563     # dim(V)^2. So it has the same dimension as the space of linear
 564     # transformations on V. In other words, it's just the right size
 565     # to create an isomorphism between it and our matrices.
 566     V = K.lattice().vector_space()
 567     W = VectorSpace(V.base_ring(), V.dimension()**2)
 568
 569     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 570
 571     # Turn our matrices into long vectors...
 572     vectors = [ W(m.list()) for m in tensor_products ]
 573
 574     # Create the *dual* cone of the positive operators, expressed as
 575     # long vectors..
 576     L = ToricLattice(W.dimension())
 577     pi_dual = Cone(vectors, lattice=L)
 578
 579     # Now compute the desired cone from its dual...
 580     pi_cone = pi_dual.dual()
 581
 582     # And finally convert its rays back to matrix representations.
 583     M = MatrixSpace(V.base_ring(), V.dimension())
 584
 585     return [ M(v.list()) for v in pi_cone.rays() ]
 586
 587
 588 def Z_transformations(K):
 589     r"""
 590     Compute generators of the cone of Z-transformations on this cone.
 591
 592     OUTPUT:
 593
 594     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 595     Each matrix ``L`` in the list should have the property that
 596     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 597     discrete complementarity set of ``K``. Moreover, any nonnegative
 598     linear combination of these matrices shares the same property.
 599
 600     EXAMPLES:
 601
 602     Z-transformations on the nonnegative orthant are just Z-matrices.
 603     That is, matrices whose off-diagonal elements are nonnegative::
 604
 605         sage: K = Cone([(1,0),(0,1)])
 606         sage: Z_transformations(K)
 607         [
 608         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 609         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 610         ]
 611         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 612         sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
 613         ....:                    for i in range(z.nrows())
 614         ....:                    for j in range(z.ncols())
 615         ....:                    if i != j ])
 616         True
 617
 618     The trivial cone in a trivial space has no Z-transformations::
 619
 620         sage: K = Cone([], ToricLattice(0))
 621         sage: Z_transformations(K)
 622         []
 623
 624     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 625
 626         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 627         sage: K.is_full_space()
 628         True
 629         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 630         sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
 631         sage: zs == lls
 632         True
 633
 634     TESTS:
 635
 636     The Z-property is possessed by every Z-transformation::
 637
 638         sage: set_random_seed()
 639         sage: K = random_cone(max_ambient_dim = 6)
 640         sage: Z_of_K = Z_transformations(K)
 641         sage: dcs = K.discrete_complementarity_set()
 642         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 643         ....:                                  for (x,s) in dcs])
 644         True
 645
 646     The lineality space of Z is LL::
 647
 648         sage: set_random_seed()
 649         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 650         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 651         sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
 652         sage: z_cone.linear_subspace() == lls
 653         True
 654
 655     """
 656     # Sage doesn't think matrices are vectors, so we have to convert
 657     # our matrices to vectors explicitly before we can figure out how
 658     # many are linearly-indepenedent.
 659     #
 660     # The space W has the same base ring as V, but dimension
 661     # dim(V)^2. So it has the same dimension as the space of linear
 662     # transformations on V. In other words, it's just the right size
 663     # to create an isomorphism between it and our matrices.
 664     V = K.lattice().vector_space()
 665     W = VectorSpace(V.base_ring(), V.dimension()**2)
 666
 667     C_of_K = K.discrete_complementarity_set()
 668     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 669
 670     # Turn our matrices into long vectors...
 671     vectors = [ W(m.list()) for m in tensor_products ]
 672
 673     # Create the *dual* cone of the cross-positive operators,
 674     # expressed as long vectors..
 675     L = ToricLattice(W.dimension())
 676     Sigma_dual = Cone(vectors, lattice=L)
 677
 678     # Now compute the desired cone from its dual...
 679     Sigma_cone = Sigma_dual.dual()
 680
 681     # And finally convert its rays back to matrix representations.
 682     # But first, make them negative, so we get Z-transformations and
 683     # not cross-positive ones.
 684     M = MatrixSpace(V.base_ring(), V.dimension())
 685
 686     return [ -M(v.list()) for v in Sigma_cone.rays() ]