mjo/cone/cone.py

   1 from sage.all import *
   2
   3 def is_lyapunov_like(L,K):
   4     r"""
   5     Determine whether or not ``L`` is Lyapunov-like on ``K``.
   6
   7     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
   8     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
   9     `\left\langle x,s \right\rangle` in the complementarity set of
  10     ``K``. It is known [Orlitzky]_ that this property need only be
  11     checked for generators of ``K`` and its dual.
  12
  13     INPUT:
  14
  15     - ``L`` -- A linear transformation or matrix.
  16
  17     - ``K`` -- A polyhedral closed convex cone.
  18
  19     OUTPUT:
  20
  21     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
  22     and ``False`` otherwise.
  23
  24     .. WARNING::
  25
  26         If this function returns ``True``, then ``L`` is Lyapunov-like
  27         on ``K``. However, if ``False`` is returned, that could mean one
  28         of two things. The first is that ``L`` is definitely not
  29         Lyapunov-like on ``K``. The second is more of an "I don't know"
  30         answer, returned (for example) if we cannot prove that an inner
  31         product is zero.
  32
  33     REFERENCES:
  34
  35     M. Orlitzky. The Lyapunov rank of an improper cone.
  36     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
  37
  38     EXAMPLES:
  39
  40     The identity is always Lyapunov-like in a nontrivial space::
  41
  42         sage: set_random_seed()
  43         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  44         sage: L = identity_matrix(K.lattice_dim())
  45         sage: is_lyapunov_like(L,K)
  46         True
  47
  48     As is the "zero" transformation::
  49
  50         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  51         sage: R = K.lattice().vector_space().base_ring()
  52         sage: L = zero_matrix(R, K.lattice_dim())
  53         sage: is_lyapunov_like(L,K)
  54         True
  55
  56         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
  57         on ``K``::
  58
  59         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
  60         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
  61         True
  62
  63     """
  64     return all([(L*x).inner_product(s) == 0
  65                 for (x,s) in K.discrete_complementarity_set()])
  66
  67
  68 def motzkin_decomposition(K):
  69     r"""
  70     Return the pair of components in the Motzkin decomposition of this cone.
  71
  72     Every convex cone is the direct sum of a strictly convex cone and a
  73     linear subspace [Stoer-Witzgall]_. Return a pair ``(P,S)`` of cones
  74     such that ``P`` is strictly convex, ``S`` is a subspace, and ``K``
  75     is the direct sum of ``P`` and ``S``.
  76
  77     OUTPUT:
  78
  79     An ordered pair ``(P,S)`` of closed convex polyhedral cones where
  80     ``P`` is strictly convex, ``S`` is a subspace, and ``K`` is the
  81     direct sum of ``P`` and ``S``.
  82
  83     REFERENCES:
  84
  85     .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and
  86        Optimization in Finite Dimensions I. Springer-Verlag, New
  87        York, 1970.
  88
  89     EXAMPLES:
  90
  91     The nonnegative orthant is strictly convex, so it is its own
  92     strictly convex component and its subspace component is trivial::
  93
  94         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  95         sage: (P,S) = motzkin_decomposition(K)
  96         sage: K.is_equivalent(P)
  97         True
  98         sage: S.is_trivial()
  99         True
 100
 101     Likewise, full spaces are their own subspace components::
 102
 103         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 104         sage: K.is_full_space()
 105         True
 106         sage: (P,S) = motzkin_decomposition(K)
 107         sage: K.is_equivalent(S)
 108         True
 109         sage: P.is_trivial()
 110         True
 111
 112     TESTS:
 113
 114     A random point in the cone should belong to either the strictly
 115     convex component or the subspace component. If the point is nonzero,
 116     it cannot be in both::
 117
 118         sage: set_random_seed()
 119         sage: K = random_cone(max_ambient_dim=8)
 120         sage: (P,S) = motzkin_decomposition(K)
 121         sage: x = K.random_element(ring=QQ)
 122         sage: P.contains(x) or S.contains(x)
 123         True
 124         sage: x.is_zero() or (P.contains(x) != S.contains(x))
 125         True
 126
 127     The strictly convex component should always be strictly convex, and
 128     the subspace component should always be a subspace::
 129
 130         sage: set_random_seed()
 131         sage: K = random_cone(max_ambient_dim=8)
 132         sage: (P,S) = motzkin_decomposition(K)
 133         sage: P.is_strictly_convex()
 134         True
 135         sage: S.lineality() == S.dim()
 136         True
 137
 138     A strictly convex cone should be equal to its strictly convex component::
 139
 140         sage: set_random_seed()
 141         sage: K = random_cone(max_ambient_dim=8, strictly_convex=True)
 142         sage: (P,_) = motzkin_decomposition(K)
 143         sage: K.is_equivalent(P)
 144         True
 145
 146     The generators of the components are obtained from orthogonal
 147     projections of the original generators [Stoer-Witzgall]_::
 148
 149         sage: set_random_seed()
 150         sage: K = random_cone(max_ambient_dim=8)
 151         sage: (P,S) = motzkin_decomposition(K)
 152         sage: A = S.linear_subspace().complement().matrix()
 153         sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A
 154         sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice())
 155         sage: P.is_equivalent(expected_P)
 156         True
 157         sage: A = S.linear_subspace().matrix()
 158         sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A
 159         sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice())
 160         sage: S.is_equivalent(expected_S)
 161         True
 162     """
 163     # The lines() method only returns one generator per line. For a true
 164     # line, we also need a generator pointing in the opposite direction.
 165     S_gens = [ direction*gen for direction in [1,-1] for gen in K.lines() ]
 166     S = Cone(S_gens, K.lattice(), check=False)
 167
 168     # Since ``S`` is a subspace, the rays of its dual generate its
 169     # orthogonal complement.
 170     S_perp = Cone(S.dual(), K.lattice(), check=False)
 171     P = K.intersection(S_perp)
 172
 173     return (P,S)
 174
 175
 176 def positive_operator_gens(K):
 177     r"""
 178     Compute generators of the cone of positive operators on this cone.
 179
 180     OUTPUT:
 181
 182     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 183     Each matrix ``P`` in the list should have the property that ``P*x``
 184     is an element of ``K`` whenever ``x`` is an element of
 185     ``K``. Moreover, any nonnegative linear combination of these
 186     matrices shares the same property.
 187
 188     EXAMPLES:
 189
 190     Positive operators on the nonnegative orthant are nonnegative matrices::
 191
 192         sage: K = Cone([(1,)])
 193         sage: positive_operator_gens(K)
 194         [[1]]
 195
 196         sage: K = Cone([(1,0),(0,1)])
 197         sage: positive_operator_gens(K)
 198         [
 199         [1 0]  [0 1]  [0 0]  [0 0]
 200         [0 0], [0 0], [1 0], [0 1]
 201         ]
 202
 203     The trivial cone in a trivial space has no positive operators::
 204
 205         sage: K = Cone([], ToricLattice(0))
 206         sage: positive_operator_gens(K)
 207         []
 208
 209     Every operator is positive on the trivial cone::
 210
 211         sage: K = Cone([(0,)])
 212         sage: positive_operator_gens(K)
 213         [[1], [-1]]
 214
 215         sage: K = Cone([(0,0)])
 216         sage: K.is_trivial()
 217         True
 218         sage: positive_operator_gens(K)
 219         [
 220         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 221         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 222         ]
 223
 224     Every operator is positive on the ambient vector space::
 225
 226         sage: K = Cone([(1,),(-1,)])
 227         sage: K.is_full_space()
 228         True
 229         sage: positive_operator_gens(K)
 230         [[1], [-1]]
 231
 232         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 233         sage: K.is_full_space()
 234         True
 235         sage: positive_operator_gens(K)
 236         [
 237         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 238         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 239         ]
 240
 241     A non-obvious application is to find the positive operators on the
 242     right half-plane::
 243
 244         sage: K = Cone([(1,0),(0,1),(0,-1)])
 245         sage: positive_operator_gens(K)
 246         [
 247         [1 0]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 248         [0 0], [1 0], [-1  0], [0 1], [ 0 -1]
 249         ]
 250
 251     TESTS:
 252
 253     Each positive operator generator should send the generators of the
 254     cone into the cone::
 255
 256         sage: set_random_seed()
 257         sage: K = random_cone(max_ambient_dim=4)
 258         sage: pi_of_K = positive_operator_gens(K)
 259         sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
 260         True
 261
 262     Each positive operator generator should send a random element of the
 263     cone into the cone::
 264
 265         sage: set_random_seed()
 266         sage: K = random_cone(max_ambient_dim=4)
 267         sage: pi_of_K = positive_operator_gens(K)
 268         sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
 269         True
 270
 271     A random element of the positive operator cone should send the
 272     generators of the cone into the cone::
 273
 274         sage: set_random_seed()
 275         sage: K = random_cone(max_ambient_dim=4)
 276         sage: pi_of_K = positive_operator_gens(K)
 277         sage: L = ToricLattice(K.lattice_dim()**2)
 278         sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
 279         ....:                lattice=L,
 280         ....:                check=False)
 281         sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
 282         sage: all([ K.contains(P*x) for x in K ])
 283         True
 284
 285     A random element of the positive operator cone should send a random
 286     element of the cone into the cone::
 287
 288         sage: set_random_seed()
 289         sage: K = random_cone(max_ambient_dim=4)
 290         sage: pi_of_K = positive_operator_gens(K)
 291         sage: L = ToricLattice(K.lattice_dim()**2)
 292         sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
 293         ....:                lattice=L,
 294         ....:                check=False)
 295         sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
 296         sage: K.contains(P*K.random_element(ring=QQ))
 297         True
 298
 299     The lineality space of the dual of the cone of positive operators
 300     can be computed from the lineality spaces of the cone and its dual::
 301
 302         sage: set_random_seed()
 303         sage: K = random_cone(max_ambient_dim=4)
 304         sage: pi_of_K = positive_operator_gens(K)
 305         sage: L = ToricLattice(K.lattice_dim()**2)
 306         sage: pi_cone = Cone([ g.list() for g in pi_of_K ],
 307         ....:                lattice=L,
 308         ....:                check=False)
 309         sage: actual = pi_cone.dual().linear_subspace()
 310         sage: U1 = [ vector((s.tensor_product(x)).list())
 311         ....:        for x in K.lines()
 312         ....:        for s in K.dual() ]
 313         sage: U2 = [ vector((s.tensor_product(x)).list())
 314         ....:        for x in K
 315         ....:        for s in K.dual().lines() ]
 316         sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
 317         sage: actual == expected
 318         True
 319
 320     The lineality of the dual of the cone of positive operators
 321     is known from its lineality space::
 322
 323         sage: set_random_seed()
 324         sage: K = random_cone(max_ambient_dim=4)
 325         sage: n = K.lattice_dim()
 326         sage: m = K.dim()
 327         sage: l = K.lineality()
 328         sage: pi_of_K = positive_operator_gens(K)
 329         sage: L = ToricLattice(n**2)
 330         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 331         ....:                 lattice=L,
 332         ....:                 check=False)
 333         sage: actual = pi_cone.dual().lineality()
 334         sage: expected = l*(m - l) + m*(n - m)
 335         sage: actual == expected
 336         True
 337
 338     The dimension of the cone of positive operators is given by the
 339     corollary in my paper::
 340
 341         sage: set_random_seed()
 342         sage: K = random_cone(max_ambient_dim=4)
 343         sage: n = K.lattice_dim()
 344         sage: m = K.dim()
 345         sage: l = K.lineality()
 346         sage: pi_of_K = positive_operator_gens(K)
 347         sage: L = ToricLattice(n**2)
 348         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 349         ....:                lattice=L,
 350         ....:                check=False)
 351         sage: actual = pi_cone.dim()
 352         sage: expected = n**2 - l*(m - l) - (n - m)*m
 353         sage: actual == expected
 354         True
 355
 356     The trivial cone, full space, and half-plane all give rise to the
 357     expected dimensions::
 358
 359         sage: n = ZZ.random_element().abs()
 360         sage: K = Cone([[0] * n], ToricLattice(n))
 361         sage: K.is_trivial()
 362         True
 363         sage: L = ToricLattice(n^2)
 364         sage: pi_of_K = positive_operator_gens(K)
 365         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 366         ....:                lattice=L,
 367         ....:                check=False)
 368         sage: actual = pi_cone.dim()
 369         sage: actual == n^2
 370         True
 371         sage: K = K.dual()
 372         sage: K.is_full_space()
 373         True
 374         sage: pi_of_K = positive_operator_gens(K)
 375         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 376         ....:                lattice=L,
 377         ....:                check=False)
 378         sage: actual = pi_cone.dim()
 379         sage: actual == n^2
 380         True
 381         sage: K = Cone([(1,0),(0,1),(0,-1)])
 382         sage: pi_of_K = positive_operator_gens(K)
 383         sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim()
 384         sage: actual == 3
 385         True
 386
 387     The lineality of the cone of positive operators follows from the
 388     description of its generators::
 389
 390         sage: set_random_seed()
 391         sage: K = random_cone(max_ambient_dim=4)
 392         sage: n = K.lattice_dim()
 393         sage: pi_of_K = positive_operator_gens(K)
 394         sage: L = ToricLattice(n**2)
 395         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 396         ....:                lattice=L,
 397         ....:                check=False)
 398         sage: actual = pi_cone.lineality()
 399         sage: expected = n**2 - K.dim()*K.dual().dim()
 400         sage: actual == expected
 401         True
 402
 403     The trivial cone, full space, and half-plane all give rise to the
 404     expected linealities::
 405
 406         sage: n = ZZ.random_element().abs()
 407         sage: K = Cone([[0] * n], ToricLattice(n))
 408         sage: K.is_trivial()
 409         True
 410         sage: L = ToricLattice(n^2)
 411         sage: pi_of_K = positive_operator_gens(K)
 412         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 413         ....:                lattice=L,
 414         ....:                check=False)
 415         sage: actual = pi_cone.lineality()
 416         sage: actual == n^2
 417         True
 418         sage: K = K.dual()
 419         sage: K.is_full_space()
 420         True
 421         sage: pi_of_K = positive_operator_gens(K)
 422         sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
 423         sage: pi_cone.lineality() == n^2
 424         True
 425         sage: K = Cone([(1,0),(0,1),(0,-1)])
 426         sage: pi_of_K = positive_operator_gens(K)
 427         sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False)
 428         sage: actual = pi_cone.lineality()
 429         sage: actual == 2
 430         True
 431
 432     A cone is proper if and only if its cone of positive operators
 433     is proper::
 434
 435         sage: set_random_seed()
 436         sage: K = random_cone(max_ambient_dim=4)
 437         sage: pi_of_K = positive_operator_gens(K)
 438         sage: L = ToricLattice(K.lattice_dim()**2)
 439         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 440         ....:                lattice=L,
 441         ....:                check=False)
 442         sage: K.is_proper() == pi_cone.is_proper()
 443         True
 444
 445     The positive operators of a permuted cone can be obtained by
 446     conjugation::
 447
 448         sage: set_random_seed()
 449         sage: K = random_cone(max_ambient_dim=4)
 450         sage: L = ToricLattice(K.lattice_dim()**2)
 451         sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
 452         sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
 453         sage: pi_of_pK = positive_operator_gens(pK)
 454         sage: actual = Cone([t.list() for t in pi_of_pK],
 455         ....:                lattice=L,
 456         ....:                check=False)
 457         sage: pi_of_K = positive_operator_gens(K)
 458         sage: expected = Cone([(p*t*p.inverse()).list() for t in pi_of_K],
 459         ....:                   lattice=L,
 460         ....:                   check=False)
 461         sage: actual.is_equivalent(expected)
 462         True
 463
 464     A transformation is positive on a cone if and only if its adjoint is
 465     positive on the dual of that cone::
 466
 467         sage: set_random_seed()
 468         sage: K = random_cone(max_ambient_dim=4)
 469         sage: F = K.lattice().vector_space().base_field()
 470         sage: n = K.lattice_dim()
 471         sage: L = ToricLattice(n**2)
 472         sage: W = VectorSpace(F, n**2)
 473         sage: pi_of_K = positive_operator_gens(K)
 474         sage: pi_of_K_star = positive_operator_gens(K.dual())
 475         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 476         ....:                lattice=L,
 477         ....:                check=False)
 478         sage: pi_star = Cone([p.list() for p in pi_of_K_star],
 479         ....:                lattice=L,
 480         ....:                check=False)
 481         sage: M = MatrixSpace(F, n)
 482         sage: L = M(pi_cone.random_element(ring=QQ).list())
 483         sage: pi_star.contains(W(L.transpose().list()))
 484         True
 485
 486         sage: L = W.random_element()
 487         sage: L_star = W(M(L.list()).transpose().list())
 488         sage: pi_cone.contains(L) ==  pi_star.contains(L_star)
 489         True
 490     """
 491     # Matrices are not vectors in Sage, so we have to convert them
 492     # to vectors explicitly before we can find a basis. We need these
 493     # two values to construct the appropriate "long vector" space.
 494     F = K.lattice().base_field()
 495     n = K.lattice_dim()
 496
 497     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 498
 499     # Convert those tensor products to long vectors.
 500     W = VectorSpace(F, n**2)
 501     vectors = [ W(tp.list()) for tp in tensor_products ]
 502
 503     check = True
 504     if K.is_solid() or K.is_strictly_convex():
 505         # The lineality space of either ``K`` or ``K.dual()`` is
 506         # trivial and it's easy to show that our generating set is
 507         # minimal. I would love a proof that this works when ``K`` is
 508         # neither pointed nor solid.
 509         #
 510         # Note that in that case we can get *duplicates*, since the
 511         # tensor product of (x,s) is the same as that of (-x,-s).
 512         check = False
 513
 514     # Create the dual cone of the positive operators, expressed as
 515     # long vectors.
 516     pi_dual = Cone(vectors, ToricLattice(W.dimension()), check=check)
 517
 518     # Now compute the desired cone from its dual...
 519     pi_cone = pi_dual.dual()
 520
 521     # And finally convert its rays back to matrix representations.
 522     M = MatrixSpace(F, n)
 523     return [ M(v.list()) for v in pi_cone ]
 524
 525
 526 def Z_transformation_gens(K):
 527     r"""
 528     Compute generators of the cone of Z-transformations on this cone.
 529
 530     OUTPUT:
 531
 532     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 533     Each matrix ``L`` in the list should have the property that
 534     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 535     discrete complementarity set of ``K``. Moreover, any nonnegative
 536     linear combination of these matrices shares the same property.
 537
 538     EXAMPLES:
 539
 540     Z-transformations on the nonnegative orthant are just Z-matrices.
 541     That is, matrices whose off-diagonal elements are nonnegative::
 542
 543         sage: K = Cone([(1,0),(0,1)])
 544         sage: Z_transformation_gens(K)
 545         [
 546         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 547         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 548         ]
 549         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 550         sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
 551         ....:                    for i in range(z.nrows())
 552         ....:                    for j in range(z.ncols())
 553         ....:                    if i != j ])
 554         True
 555
 556     The trivial cone in a trivial space has no Z-transformations::
 557
 558         sage: K = Cone([], ToricLattice(0))
 559         sage: Z_transformation_gens(K)
 560         []
 561
 562     Every operator is a Z-transformation on the ambient vector space::
 563
 564         sage: K = Cone([(1,),(-1,)])
 565         sage: K.is_full_space()
 566         True
 567         sage: Z_transformation_gens(K)
 568         [[-1], [1]]
 569
 570         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 571         sage: K.is_full_space()
 572         True
 573         sage: Z_transformation_gens(K)
 574         [
 575         [-1  0]  [1 0]  [ 0 -1]  [0 1]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
 576         [ 0  0], [0 0], [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
 577         ]
 578
 579     A non-obvious application is to find the Z-transformations on the
 580     right half-plane::
 581
 582         sage: K = Cone([(1,0),(0,1),(0,-1)])
 583         sage: Z_transformation_gens(K)
 584         [
 585         [-1  0]  [1 0]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
 586         [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
 587         ]
 588
 589     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 590
 591         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 592         sage: K.is_full_space()
 593         True
 594         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 595         sage: zs  = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
 596         sage: zs == lls
 597         True
 598
 599     TESTS:
 600
 601     The Z-property is possessed by every Z-transformation::
 602
 603         sage: set_random_seed()
 604         sage: K = random_cone(max_ambient_dim=4)
 605         sage: Z_of_K = Z_transformation_gens(K)
 606         sage: dcs = K.discrete_complementarity_set()
 607         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 608         ....:                                  for (x,s) in dcs])
 609         True
 610
 611     The lineality space of the cone of Z-transformations is the space of
 612     Lyapunov-like transformations::
 613
 614         sage: set_random_seed()
 615         sage: K = random_cone(max_ambient_dim=4)
 616         sage: L = ToricLattice(K.lattice_dim()**2)
 617         sage: Z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ],
 618         ....:               lattice=L,
 619         ....:               check=False)
 620         sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
 621         sage: lls = L.vector_space().span(ll_basis)
 622         sage: Z_cone.linear_subspace() == lls
 623         True
 624
 625     The lineality of the Z-transformations on a cone is the Lyapunov
 626     rank of that cone::
 627
 628         sage: set_random_seed()
 629         sage: K = random_cone(max_ambient_dim=4)
 630         sage: Z_of_K = Z_transformation_gens(K)
 631         sage: L = ToricLattice(K.lattice_dim()**2)
 632         sage: Z_cone  = Cone([ z.list() for z in Z_of_K ],
 633         ....:                lattice=L,
 634         ....:                check=False)
 635         sage: Z_cone.lineality() == K.lyapunov_rank()
 636         True
 637
 638     The lineality spaces of the duals of the positive operator and
 639     Z-transformation cones are equal. From this it follows that the
 640     dimensions of the Z-transformation cone and positive operator cone
 641     are equal::
 642
 643         sage: set_random_seed()
 644         sage: K = random_cone(max_ambient_dim=4)
 645         sage: pi_of_K = positive_operator_gens(K)
 646         sage: Z_of_K = Z_transformation_gens(K)
 647         sage: L = ToricLattice(K.lattice_dim()**2)
 648         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 649         ....:                lattice=L,
 650         ....:                check=False)
 651         sage: Z_cone = Cone([ z.list() for z in Z_of_K],
 652         ....:               lattice=L,
 653         ....:               check=False)
 654         sage: pi_cone.dim() == Z_cone.dim()
 655         True
 656         sage: pi_star = pi_cone.dual()
 657         sage: z_star = Z_cone.dual()
 658         sage: pi_star.linear_subspace() == z_star.linear_subspace()
 659         True
 660
 661     The trivial cone, full space, and half-plane all give rise to the
 662     expected dimensions::
 663
 664         sage: n = ZZ.random_element().abs()
 665         sage: K = Cone([[0] * n], ToricLattice(n))
 666         sage: K.is_trivial()
 667         True
 668         sage: L = ToricLattice(n^2)
 669         sage: Z_of_K = Z_transformation_gens(K)
 670         sage: Z_cone = Cone([z.list() for z in Z_of_K],
 671         ....:               lattice=L,
 672         ....:               check=False)
 673         sage: actual = Z_cone.dim()
 674         sage: actual == n^2
 675         True
 676         sage: K = K.dual()
 677         sage: K.is_full_space()
 678         True
 679         sage: Z_of_K = Z_transformation_gens(K)
 680         sage: Z_cone = Cone([z.list() for z in Z_of_K],
 681         ....:                lattice=L,
 682         ....:                check=False)
 683         sage: actual = Z_cone.dim()
 684         sage: actual == n^2
 685         True
 686         sage: K = Cone([(1,0),(0,1),(0,-1)])
 687         sage: Z_of_K = Z_transformation_gens(K)
 688         sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
 689         sage: Z_cone.dim() == 3
 690         True
 691
 692     The Z-transformations of a permuted cone can be obtained by
 693     conjugation::
 694
 695         sage: set_random_seed()
 696         sage: K = random_cone(max_ambient_dim=4)
 697         sage: L = ToricLattice(K.lattice_dim()**2)
 698         sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
 699         sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
 700         sage: Z_of_pK = Z_transformation_gens(pK)
 701         sage: actual = Cone([t.list() for t in Z_of_pK],
 702         ....:                lattice=L,
 703         ....:                check=False)
 704         sage: Z_of_K = Z_transformation_gens(K)
 705         sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K],
 706         ....:                   lattice=L,
 707         ....:                   check=False)
 708         sage: actual.is_equivalent(expected)
 709         True
 710
 711     A transformation is a Z-transformation on a cone if and only if its
 712     adjoint is a Z-transformation on the dual of that cone::
 713
 714         sage: set_random_seed()
 715         sage: K = random_cone(max_ambient_dim=4)
 716         sage: F = K.lattice().vector_space().base_field()
 717         sage: n = K.lattice_dim()
 718         sage: L = ToricLattice(n**2)
 719         sage: W = VectorSpace(F, n**2)
 720         sage: Z_of_K = Z_transformation_gens(K)
 721         sage: Z_of_K_star = Z_transformation_gens(K.dual())
 722         sage: Z_cone = Cone([p.list() for p in Z_of_K],
 723         ....:               lattice=L,
 724         ....:               check=False)
 725         sage: Z_star = Cone([p.list() for p in Z_of_K_star],
 726         ....:               lattice=L,
 727         ....:               check=False)
 728         sage: M = MatrixSpace(F, n)
 729         sage: L = M(Z_cone.random_element(ring=QQ).list())
 730         sage: Z_star.contains(W(L.transpose().list()))
 731         True
 732
 733         sage: L = W.random_element()
 734         sage: L_star = W(M(L.list()).transpose().list())
 735         sage: Z_cone.contains(L) ==  Z_star.contains(L_star)
 736         True
 737     """
 738     # Matrices are not vectors in Sage, so we have to convert them
 739     # to vectors explicitly before we can find a basis. We need these
 740     # two values to construct the appropriate "long vector" space.
 741     F = K.lattice().base_field()
 742     n = K.lattice_dim()
 743
 744     # These tensor products contain generators for the dual cone of
 745     # the cross-positive transformations.
 746     tensor_products = [ s.tensor_product(x)
 747                         for (x,s) in K.discrete_complementarity_set() ]
 748
 749     # Turn our matrices into long vectors...
 750     W = VectorSpace(F, n**2)
 751     vectors = [ W(m.list()) for m in tensor_products ]
 752
 753     check = True
 754     if K.is_solid() or K.is_strictly_convex():
 755         # The lineality space of either ``K`` or ``K.dual()`` is
 756         # trivial and it's easy to show that our generating set is
 757         # minimal. I would love a proof that this works when ``K`` is
 758         # neither pointed nor solid.
 759         #
 760         # Note that in that case we can get *duplicates*, since the
 761         # tensor product of (x,s) is the same as that of (-x,-s).
 762         check = False
 763
 764     # Create the dual cone of the cross-positive operators,
 765     # expressed as long vectors.
 766     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()), check=check)
 767
 768     # Now compute the desired cone from its dual...
 769     Sigma_cone = Sigma_dual.dual()
 770
 771     # And finally convert its rays back to matrix representations.
 772     # But first, make them negative, so we get Z-transformations and
 773     # not cross-positive ones.
 774     M = MatrixSpace(F, n)
 775     return [ -M(v.list()) for v in Sigma_cone ]
 776
 777
 778 def Z_cone(K):
 779     gens = Z_transformation_gens(K)
 780     L = ToricLattice(K.lattice_dim()**2)
 781     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 782
 783 def pi_cone(K):
 784     gens = positive_operator_gens(K)
 785     L = ToricLattice(K.lattice_dim()**2)
 786     return Cone([ g.list() for g in gens ], lattice=L, check=False)