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