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