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