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