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     REFERENCES:
  96
  97     .. [Orlitzky-Pi-Z]
  98        M. Orlitzky.
  99        Positive and Z-operators on closed convex cones.
 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     """
 437     if K2 is None:
 438         K2 = K1
 439
 440     # Matrices are not vectors in Sage, so we have to convert them
 441     # to vectors explicitly before we can find a basis. We need these
 442     # two values to construct the appropriate "long vector" space.
 443     F = K1.lattice().base_field()
 444     n = K1.lattice_dim()
 445     m = K2.lattice_dim()
 446
 447     tensor_products = [ s.tensor_product(x) for x in K1 for s in K2.dual() ]
 448
 449     # Convert those tensor products to long vectors.
 450     W = VectorSpace(F, n*m)
 451     vectors = [ W(tp.list()) for tp in tensor_products ]
 452
 453     check = True
 454     if K1.is_proper() and K2.is_proper():
 455         # All of the generators involved are extreme vectors and
 456         # therefore minimal [Tam]_. If this cone is neither solid nor
 457         # strictly convex, then the tensor product of ``s`` and ``x``
 458         # is the same as that of ``-s`` and ``-x``. However, as a
 459         # /set/, ``tensor_products`` may still be minimal.
 460         check = False
 461
 462     # Create the dual cone of the positive operators, expressed as
 463     # long vectors.
 464     pi_dual = Cone(vectors, ToricLattice(W.dimension()), check=check)
 465
 466     # Now compute the desired cone from its dual...
 467     pi_cone = pi_dual.dual()
 468
 469     # And finally convert its rays back to matrix representations.
 470     M = MatrixSpace(F, m, n)
 471     return [ M(v.list()) for v in pi_cone ]
 472
 473
 474 def Z_operator_gens(K):
 475     r"""
 476     Compute generators of the cone of Z-operators on this cone.
 477
 478     OUTPUT:
 479
 480     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 481     Each matrix ``L`` in the list should have the property that
 482     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element of
 483     this cone's :meth:`discrete_complementarity_set`. Moreover, any
 484     conic (nonnegative linear) combination of these matrices shares the
 485     same property.
 486
 487     REFERENCES:
 488
 489     M. Orlitzky.
 490     Positive and Z-operators on closed convex cones.
 491
 492     EXAMPLES:
 493
 494     Z-operators on the nonnegative orthant are just Z-matrices.
 495     That is, matrices whose off-diagonal elements are nonnegative::
 496
 497         sage: K = Cone([(1,0),(0,1)])
 498         sage: Z_operator_gens(K)
 499         [
 500         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 501         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 502         ]
 503         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 504         sage: all([ z[i][j] <= 0 for z in Z_operator_gens(K)
 505         ....:                    for i in range(z.nrows())
 506         ....:                    for j in range(z.ncols())
 507         ....:                    if i != j ])
 508         True
 509
 510     The trivial cone in a trivial space has no Z-operators::
 511
 512         sage: K = Cone([], ToricLattice(0))
 513         sage: Z_operator_gens(K)
 514         []
 515
 516     Every operator is a Z-operator on the ambient vector space::
 517
 518         sage: K = Cone([(1,),(-1,)])
 519         sage: K.is_full_space()
 520         True
 521         sage: Z_operator_gens(K)
 522         [[-1], [1]]
 523
 524         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 525         sage: K.is_full_space()
 526         True
 527         sage: Z_operator_gens(K)
 528         [
 529         [-1  0]  [1 0]  [ 0 -1]  [0 1]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
 530         [ 0  0], [0 0], [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
 531         ]
 532
 533     A non-obvious application is to find the Z-operators on the
 534     right half-plane::
 535
 536         sage: K = Cone([(1,0),(0,1),(0,-1)])
 537         sage: Z_operator_gens(K)
 538         [
 539         [-1  0]  [1 0]  [ 0  0]  [0 0]  [ 0  0]  [0 0]
 540         [ 0  0], [0 0], [-1  0], [1 0], [ 0 -1], [0 1]
 541         ]
 542
 543     Z-operators on a subspace are Lyapunov-like and vice-versa::
 544
 545         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 546         sage: K.is_full_space()
 547         True
 548         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 549         sage: zs  = span([ vector(z.list()) for z in Z_operator_gens(K) ])
 550         sage: zs == lls
 551         True
 552
 553     TESTS:
 554
 555     The Z-property is possessed by every Z-operator::
 556
 557         sage: set_random_seed()
 558         sage: K = random_cone(max_ambient_dim=4)
 559         sage: Z_of_K = Z_operator_gens(K)
 560         sage: dcs = K.discrete_complementarity_set()
 561         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 562         ....:                                  for (x,s) in dcs])
 563         True
 564
 565     The lineality space of the cone of Z-operators is the space of
 566     Lyapunov-like operators::
 567
 568         sage: set_random_seed()
 569         sage: K = random_cone(max_ambient_dim=4)
 570         sage: L = ToricLattice(K.lattice_dim()**2)
 571         sage: Z_cone = Cone([ z.list() for z in Z_operator_gens(K) ],
 572         ....:               lattice=L,
 573         ....:               check=False)
 574         sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ]
 575         sage: lls = L.vector_space().span(ll_basis)
 576         sage: Z_cone.linear_subspace() == lls
 577         True
 578
 579     The lineality of the Z-operators on a cone is the Lyapunov
 580     rank of that cone::
 581
 582         sage: set_random_seed()
 583         sage: K = random_cone(max_ambient_dim=4)
 584         sage: Z_of_K = Z_operator_gens(K)
 585         sage: L = ToricLattice(K.lattice_dim()**2)
 586         sage: Z_cone  = Cone([ z.list() for z in Z_of_K ],
 587         ....:                lattice=L,
 588         ....:                check=False)
 589         sage: Z_cone.lineality() == K.lyapunov_rank()
 590         True
 591
 592     The lineality spaces of the duals of the positive and Z-operator
 593     cones are equal. From this it follows that the dimensions of the
 594     Z-operator cone and positive operator cone are equal::
 595
 596         sage: set_random_seed()
 597         sage: K = random_cone(max_ambient_dim=4)
 598         sage: pi_of_K = positive_operator_gens(K)
 599         sage: Z_of_K = Z_operator_gens(K)
 600         sage: L = ToricLattice(K.lattice_dim()**2)
 601         sage: pi_cone = Cone([p.list() for p in pi_of_K],
 602         ....:                lattice=L,
 603         ....:                check=False)
 604         sage: Z_cone = Cone([ z.list() for z in Z_of_K],
 605         ....:               lattice=L,
 606         ....:               check=False)
 607         sage: pi_cone.dim() == Z_cone.dim()
 608         True
 609         sage: pi_star = pi_cone.dual()
 610         sage: z_star = Z_cone.dual()
 611         sage: pi_star.linear_subspace() == z_star.linear_subspace()
 612         True
 613
 614     The trivial cone, full space, and half-plane all give rise to the
 615     expected dimensions::
 616
 617         sage: n = ZZ.random_element().abs()
 618         sage: K = Cone([[0] * n], ToricLattice(n))
 619         sage: K.is_trivial()
 620         True
 621         sage: L = ToricLattice(n^2)
 622         sage: Z_of_K = Z_operator_gens(K)
 623         sage: Z_cone = Cone([z.list() for z in Z_of_K],
 624         ....:               lattice=L,
 625         ....:               check=False)
 626         sage: actual = Z_cone.dim()
 627         sage: actual == n^2
 628         True
 629         sage: K = K.dual()
 630         sage: K.is_full_space()
 631         True
 632         sage: Z_of_K = Z_operator_gens(K)
 633         sage: Z_cone = Cone([z.list() for z in Z_of_K],
 634         ....:                lattice=L,
 635         ....:                check=False)
 636         sage: actual = Z_cone.dim()
 637         sage: actual == n^2
 638         True
 639         sage: K = Cone([(1,0),(0,1),(0,-1)])
 640         sage: Z_of_K = Z_operator_gens(K)
 641         sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False)
 642         sage: Z_cone.dim() == 3
 643         True
 644
 645     The Z-operators of a permuted cone can be obtained by conjugation::
 646
 647         sage: set_random_seed()
 648         sage: K = random_cone(max_ambient_dim=4)
 649         sage: L = ToricLattice(K.lattice_dim()**2)
 650         sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix()
 651         sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False)
 652         sage: Z_of_pK = Z_operator_gens(pK)
 653         sage: actual = Cone([t.list() for t in Z_of_pK],
 654         ....:                lattice=L,
 655         ....:                check=False)
 656         sage: Z_of_K = Z_operator_gens(K)
 657         sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K],
 658         ....:                   lattice=L,
 659         ....:                   check=False)
 660         sage: actual.is_equivalent(expected)
 661         True
 662
 663     An operator is a Z-operator on a cone if and only if its
 664     adjoint is a Z-operator on the dual of that cone::
 665
 666         sage: set_random_seed()
 667         sage: K = random_cone(max_ambient_dim=4)
 668         sage: F = K.lattice().vector_space().base_field()
 669         sage: n = K.lattice_dim()
 670         sage: L = ToricLattice(n**2)
 671         sage: W = VectorSpace(F, n**2)
 672         sage: Z_of_K = Z_operator_gens(K)
 673         sage: Z_of_K_star = Z_operator_gens(K.dual())
 674         sage: Z_cone = Cone([p.list() for p in Z_of_K],
 675         ....:               lattice=L,
 676         ....:               check=False)
 677         sage: Z_star = Cone([p.list() for p in Z_of_K_star],
 678         ....:               lattice=L,
 679         ....:               check=False)
 680         sage: M = MatrixSpace(F, n)
 681         sage: L = M(Z_cone.random_element(ring=QQ).list())
 682         sage: Z_star.contains(W(L.transpose().list()))
 683         True
 684
 685         sage: L = W.random_element()
 686         sage: L_star = W(M(L.list()).transpose().list())
 687         sage: Z_cone.contains(L) ==  Z_star.contains(L_star)
 688         True
 689     """
 690     # Matrices are not vectors in Sage, so we have to convert them
 691     # to vectors explicitly before we can find a basis. We need these
 692     # two values to construct the appropriate "long vector" space.
 693     F = K.lattice().base_field()
 694     n = K.lattice_dim()
 695
 696     # These tensor products contain generators for the dual cone of
 697     # the cross-positive operators.
 698     tensor_products = [ s.tensor_product(x)
 699                         for (x,s) in K.discrete_complementarity_set() ]
 700
 701     # Turn our matrices into long vectors...
 702     W = VectorSpace(F, n**2)
 703     vectors = [ W(m.list()) for m in tensor_products ]
 704
 705     check = True
 706     if K.is_proper():
 707         # All of the generators involved are extreme vectors and
 708         # therefore minimal. If this cone is neither solid nor
 709         # strictly convex, then the tensor product of ``s`` and ``x``
 710         # is the same as that of ``-s`` and ``-x``. However, as a
 711         # /set/, ``tensor_products`` may still be minimal.
 712         check = False
 713
 714     # Create the dual cone of the cross-positive operators,
 715     # expressed as long vectors.
 716     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()), check=check)
 717
 718     # Now compute the desired cone from its dual...
 719     Sigma_cone = Sigma_dual.dual()
 720
 721     # And finally convert its rays back to matrix representations.
 722     # But first, make them negative, so we get Z-operators and
 723     # not cross-positive ones.
 724     M = MatrixSpace(F, n)
 725     return [ -M(v.list()) for v in Sigma_cone ]
 726
 727
 728 def LL_cone(K):
 729     gens = K.lyapunov_like_basis()
 730     L = ToricLattice(K.lattice_dim()**2)
 731     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 732
 733 def Z_cone(K):
 734     gens = Z_operator_gens(K)
 735     L = ToricLattice(K.lattice_dim()**2)
 736     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 737
 738 def pi_cone(K):
 739     gens = positive_operator_gens(K)
 740     L = ToricLattice(K.lattice_dim()**2)
 741     return Cone([ g.list() for g in gens ], lattice=L, check=False)