mjo/cone/cone.py

   1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
   2 # have to explicitly mangle our sitedir here so that "mjo.cone"
   3 # resolves.
   4 from os.path import abspath
   5 from site import addsitedir
   6 addsitedir(abspath('../../'))
   7
   8 from sage.all import *
   9
  10
  11 def _basically_the_same(K1, K2):
  12     r"""
  13     Test whether or not ``K1`` and ``K2`` are "basically the same."
  14
  15     This is a hack to get around the fact that it's difficult to tell
  16     when two cones are linearly isomorphic. We have a proposition that
  17     equates two cones, but represented over `\mathbb{Q}`, they are
  18     merely linearly isomorphic (not equal). So rather than test for
  19     equality, we test a list of properties that should be preserved
  20     under an invertible linear transformation.
  21
  22     OUTPUT:
  23
  24     ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
  25     otherwise.
  26
  27     EXAMPLES:
  28
  29     Any proper cone with three generators in `\mathbb{R}^{3}` is
  30     basically the same as the nonnegative orthant::
  31
  32         sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
  33         sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
  34         sage: _basically_the_same(K1, K2)
  35         True
  36
  37     Negating a cone gives you another cone that is basically the same::
  38
  39         sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
  40         sage: _basically_the_same(K, -K)
  41         True
  42
  43     TESTS:
  44
  45     Any cone is basically the same as itself::
  46
  47         sage: K = random_cone(max_ambient_dim = 8)
  48         sage: _basically_the_same(K, K)
  49         True
  50
  51     After applying an invertible matrix to the rows of a cone, the
  52     result should be basically the same as the cone we started with::
  53
  54         sage: K1 = random_cone(max_ambient_dim = 8)
  55         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
  56         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
  57         sage: _basically_the_same(K1, K2)
  58         True
  59
  60     """
  61     if K1.lattice_dim() != K2.lattice_dim():
  62         return False
  63
  64     if K1.nrays() != K2.nrays():
  65         return False
  66
  67     if K1.dim() != K2.dim():
  68         return False
  69
  70     if K1.lineality() != K2.lineality():
  71         return False
  72
  73     if K1.is_solid() != K2.is_solid():
  74         return False
  75
  76     if K1.is_strictly_convex() != K2.is_strictly_convex():
  77         return False
  78
  79     if len(K1.LL()) != len(K2.LL()):
  80         return False
  81
  82     C_of_K1 = K1.discrete_complementarity_set()
  83     C_of_K2 = K2.discrete_complementarity_set()
  84     if len(C_of_K1) != len(C_of_K2):
  85         return False
  86
  87     if len(K1.facets()) != len(K2.facets()):
  88         return False
  89
  90     return True
  91
  92
  93
  94 def _restrict_to_space(K, W):
  95     r"""
  96     Restrict this cone a subspace of its ambient space.
  97
  98     INPUT:
  99
 100     - ``W`` -- The subspace into which this cone will be restricted.
 101
 102     OUTPUT:
 103
 104     A new cone in a sublattice corresponding to ``W``.
 105
 106     EXAMPLES:
 107
 108     When this cone is solid, restricting it into its own span should do
 109     nothing::
 110
 111         sage: K = Cone([(1,)])
 112         sage: _restrict_to_space(K, K.span()) == K
 113         True
 114
 115     A single ray restricted into its own span gives the same output
 116     regardless of the ambient space::
 117
 118         sage: K2 = Cone([(1,0)])
 119         sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
 120         sage: K2_S
 121         N(1)
 122         in 1-d lattice N
 123         sage: K3 = Cone([(1,0,0)])
 124         sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
 125         sage: K3_S
 126         N(1)
 127         in 1-d lattice N
 128         sage: K2_S == K3_S
 129         True
 130
 131     TESTS:
 132
 133     The projected cone should always be solid::
 134
 135         sage: set_random_seed()
 136         sage: K = random_cone(max_ambient_dim = 8)
 137         sage: _restrict_to_space(K, K.span()).is_solid()
 138         True
 139
 140     And the resulting cone should live in a space having the same
 141     dimension as the space we restricted it to::
 142
 143         sage: set_random_seed()
 144         sage: K = random_cone(max_ambient_dim = 8)
 145         sage: K_P = _restrict_to_space(K, K.dual().span())
 146         sage: K_P.lattice_dim() == K.dual().dim()
 147         True
 148
 149     This function should not affect the dimension of a cone::
 150
 151         sage: set_random_seed()
 152         sage: K = random_cone(max_ambient_dim = 8)
 153         sage: K.dim() == _restrict_to_space(K,K.span()).dim()
 154         True
 155
 156     Nor should it affect the lineality of a cone::
 157
 158         sage: set_random_seed()
 159         sage: K = random_cone(max_ambient_dim = 8)
 160         sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
 161         True
 162
 163     No matter which space we restrict to, the lineality should not
 164     increase::
 165
 166         sage: set_random_seed()
 167         sage: K = random_cone(max_ambient_dim = 8)
 168         sage: S = K.span(); P = K.dual().span()
 169         sage: K.lineality() >= _restrict_to_space(K,S).lineality()
 170         True
 171         sage: K.lineality() >= _restrict_to_space(K,P).lineality()
 172         True
 173
 174     If we do this according to our paper, then the result is proper::
 175
 176         sage: set_random_seed()
 177         sage: K = random_cone(max_ambient_dim = 8)
 178         sage: K_S = _restrict_to_space(K, K.span())
 179         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 180         sage: K_SP.is_proper()
 181         True
 182         sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
 183         sage: K_SP.is_proper()
 184         True
 185
 186     Test the proposition in our paper concerning the duals and
 187     restrictions. Generate a random cone, then create a subcone of
 188     it. The operation of dual-taking should then commute with
 189     _restrict_to_space::
 190
 191         sage: set_random_seed()
 192         sage: J = random_cone(max_ambient_dim = 8)
 193         sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
 194         sage: K_W_star = _restrict_to_space(K, J.span()).dual()
 195         sage: K_star_W = _restrict_to_space(K.dual(), J.span())
 196         sage: _basically_the_same(K_W_star, K_star_W)
 197         True
 198
 199     """
 200     # First we want to intersect ``K`` with ``W``. The easiest way to
 201     # do this is via cone intersection, so we turn the subspace ``W``
 202     # into a cone.
 203     W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice())
 204     K = K.intersection(W_cone)
 205
 206     # We've already intersected K with the span of K2, so every
 207     # generator of K should belong to W now.
 208     K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
 209
 210     L = ToricLattice(W.dimension())
 211     return Cone(K_W_rays, lattice=L)
 212
 213
 214 def lyapunov_rank(K):
 215     r"""
 216     Compute the Lyapunov rank (or bilinearity rank) of this cone.
 217
 218     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 219
 220     1. The dimension of the Lie algebra of the automorphism group of the
 221        cone.
 222
 223     2. The dimension of the linear space of all Lyapunov-like
 224        transformations on the cone.
 225
 226     INPUT:
 227
 228     A closed, convex polyhedral cone.
 229
 230     OUTPUT:
 231
 232     An integer representing the Lyapunov rank of the cone. If the
 233     dimension of the ambient vector space is `n`, then the Lyapunov rank
 234     will be between `1` and `n` inclusive; however a rank of `n-1` is
 235     not possible (see [Orlitzky/Gowda]_).
 236
 237     ALGORITHM:
 238
 239     The codimension formula from the second reference is used. We find
 240     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 241     and `s` are rays of our cone. It is known that these vectors are
 242     sufficient to apply the codimension formula. Once we have all such
 243     pairs, we "brute force" the codimension formula by finding all
 244     linearly-independent `xs^{T}`.
 245
 246     REFERENCES:
 247
 248     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 249        cone and Lyapunov-like transformations, Mathematical Programming, 147
 250        (2014) 155-170.
 251
 252     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 253        Improper Cone. Work in-progress.
 254
 255     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 256        optimality constraints for the cone of positive polynomials,
 257        Mathematical Programming, Series B, 129 (2011) 5-31.
 258
 259     EXAMPLES:
 260
 261     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 262     [Rudolf et al.]_::
 263
 264         sage: positives = Cone([(1,)])
 265         sage: lyapunov_rank(positives)
 266         1
 267         sage: quadrant = Cone([(1,0), (0,1)])
 268         sage: lyapunov_rank(quadrant)
 269         2
 270        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 271         sage: lyapunov_rank(octant)
 272         3
 273
 274     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 275     [Orlitzky/Gowda]_::
 276
 277         sage: R5 = VectorSpace(QQ, 5)
 278         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
 279         sage: K = Cone(gs)
 280         sage: lyapunov_rank(K)
 281         25
 282
 283     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 284     [Rudolf et al.]_::
 285
 286         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 287         sage: lyapunov_rank(L31)
 288         1
 289
 290     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 291
 292         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 293         sage: lyapunov_rank(L3infty)
 294         1
 295
 296     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 297     + 1` [Orlitzky/Gowda]_::
 298
 299         sage: K = Cone([(1,0,0,0,0)])
 300         sage: lyapunov_rank(K)
 301         21
 302         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 303         21
 304
 305     A subspace (of dimension `m`) in `n` dimensions should have a
 306     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 307
 308         sage: e1 = (1,0,0,0,0)
 309         sage: neg_e1 = (-1,0,0,0,0)
 310         sage: e2 = (0,1,0,0,0)
 311         sage: neg_e2 = (0,-1,0,0,0)
 312         sage: z = (0,0,0,0,0)
 313         sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
 314         sage: lyapunov_rank(K)
 315         19
 316         sage: K.lattice_dim()**2 - K.dim()*K.codim()
 317         19
 318
 319     The Lyapunov rank should be additive on a product of proper cones
 320     [Rudolf et al.]_::
 321
 322         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 323         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 324         sage: K = L31.cartesian_product(octant)
 325         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 326         True
 327
 328     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 329     The cone ``K`` in the following example is isomorphic to the nonnegative
 330     octant in `\mathbb{R}^{3}`::
 331
 332         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 333         sage: lyapunov_rank(K)
 334         3
 335
 336     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 337     itself [Rudolf et al.]_::
 338
 339         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 340         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 341         True
 342
 343     TESTS:
 344
 345     The Lyapunov rank should be additive on a product of proper cones
 346     [Rudolf et al.]_::
 347
 348         sage: set_random_seed()
 349         sage: K1 = random_cone(max_ambient_dim=8,
 350         ....:                  strictly_convex=True,
 351         ....:                  solid=True)
 352         sage: K2 = random_cone(max_ambient_dim=8,
 353         ....:                  strictly_convex=True,
 354         ....:                  solid=True)
 355         sage: K = K1.cartesian_product(K2)
 356         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 357         True
 358
 359     The Lyapunov rank is invariant under a linear isomorphism
 360     [Orlitzky/Gowda]_::
 361
 362         sage: K1 = random_cone(max_ambient_dim = 8)
 363         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
 364         sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
 365         sage: lyapunov_rank(K1) == lyapunov_rank(K2)
 366         True
 367
 368     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 369     itself [Rudolf et al.]_::
 370
 371         sage: set_random_seed()
 372         sage: K = random_cone(max_ambient_dim=8)
 373         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 374         True
 375
 376     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 377     be any number between `1` and `n` inclusive, excluding `n-1`
 378     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 379     trivial cone in a trivial space as well. However, in zero dimensions,
 380     the Lyapunov rank of the trivial cone will be zero::
 381
 382         sage: set_random_seed()
 383         sage: K = random_cone(max_ambient_dim=8,
 384         ....:                 strictly_convex=True,
 385         ....:                 solid=True)
 386         sage: b = lyapunov_rank(K)
 387         sage: n = K.lattice_dim()
 388         sage: (n == 0 or 1 <= b) and b <= n
 389         True
 390         sage: b == n-1
 391         False
 392
 393     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 394     Lyapunov rank `n-1` in `n` dimensions::
 395
 396         sage: set_random_seed()
 397         sage: K = random_cone(max_ambient_dim=8)
 398         sage: b = lyapunov_rank(K)
 399         sage: n = K.lattice_dim()
 400         sage: b == n-1
 401         False
 402
 403     The calculation of the Lyapunov rank of an improper cone can be
 404     reduced to that of a proper cone [Orlitzky/Gowda]_::
 405
 406         sage: set_random_seed()
 407         sage: K = random_cone(max_ambient_dim=8)
 408         sage: actual = lyapunov_rank(K)
 409         sage: K_S = _restrict_to_space(K, K.span())
 410         sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
 411         sage: l = K.lineality()
 412         sage: c = K.codim()
 413         sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
 414         sage: actual == expected
 415         True
 416
 417     The Lyapunov rank of any cone is just the dimension of ``K.LL()``::
 418
 419         sage: set_random_seed()
 420         sage: K = random_cone(max_ambient_dim=8)
 421         sage: lyapunov_rank(K) == len(K.LL())
 422         True
 423
 424     We can make an imperfect cone perfect by adding a slack variable
 425     (a Theorem in [Orlitzky/Gowda]_)::
 426
 427         sage: set_random_seed()
 428         sage: K = random_cone(max_ambient_dim=8,
 429         ....:                 strictly_convex=True,
 430         ....:                 solid=True)
 431         sage: L = ToricLattice(K.lattice_dim() + 1)
 432         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 433         sage: lyapunov_rank(K) >= K.lattice_dim()
 434         True
 435
 436     """
 437     beta = 0
 438
 439     m = K.dim()
 440     n = K.lattice_dim()
 441     l = K.lineality()
 442
 443     if m < n:
 444         # K is not solid, restrict to its span.
 445         K = _restrict_to_space(K, K.span())
 446
 447         # Non-solid reduction lemma.
 448         beta += (n - m)*n
 449
 450     if l > 0:
 451         # K is not pointed, restrict to the span of its dual. Uses a
 452         # proposition from our paper, i.e. this is equivalent to K =
 453         # _rho(K.dual()).dual().
 454         K = _restrict_to_space(K, K.dual().span())
 455
 456         # Non-pointed reduction lemma.
 457         beta += l * m
 458
 459     beta += len(K.LL())
 460     return beta
 461
 462
 463
 464 def is_lyapunov_like(L,K):
 465     r"""
 466     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 467
 468     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 469     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 470     `\left\langle x,s \right\rangle` in the complementarity set of
 471     ``K``. It is known [Orlitzky]_ that this property need only be
 472     checked for generators of ``K`` and its dual.
 473
 474     INPUT:
 475
 476     - ``L`` -- A linear transformation or matrix.
 477
 478     - ``K`` -- A polyhedral closed convex cone.
 479
 480     OUTPUT:
 481
 482     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 483     and ``False`` otherwise.
 484
 485     .. WARNING::
 486
 487         If this function returns ``True``, then ``L`` is Lyapunov-like
 488         on ``K``. However, if ``False`` is returned, that could mean one
 489         of two things. The first is that ``L`` is definitely not
 490         Lyapunov-like on ``K``. The second is more of an "I don't know"
 491         answer, returned (for example) if we cannot prove that an inner
 492         product is zero.
 493
 494     REFERENCES:
 495
 496     .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
 497        improper cone (preprint).
 498
 499     EXAMPLES:
 500
 501     The identity is always Lyapunov-like in a nontrivial space::
 502
 503         sage: set_random_seed()
 504         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
 505         sage: L = identity_matrix(K.lattice_dim())
 506         sage: is_lyapunov_like(L,K)
 507         True
 508
 509     As is the "zero" transformation::
 510
 511         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 512         sage: R = K.lattice().vector_space().base_ring()
 513         sage: L = zero_matrix(R, K.lattice_dim())
 514         sage: is_lyapunov_like(L,K)
 515         True
 516
 517     Everything in ``K.LL()`` should be Lyapunov-like on ``K``::
 518
 519         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 520         sage: all([is_lyapunov_like(L,K) for L in K.LL()])
 521         True
 522
 523     """
 524     return all([(L*x).inner_product(s) == 0
 525                 for (x,s) in K.discrete_complementarity_set()])
 526
 527
 528 def random_element(K):
 529     r"""
 530     Return a random element of ``K`` from its ambient vector space.
 531
 532     ALGORITHM:
 533
 534     The cone ``K`` is specified in terms of its generators, so that
 535     ``K`` is equal to the convex conic combination of those generators.
 536     To choose a random element of ``K``, we assign random nonnegative
 537     coefficients to each generator of ``K`` and construct a new vector
 538     from the scaled rays.
 539
 540     A vector, rather than a ray, is returned so that the element may
 541     have non-integer coordinates. Thus the element may have an
 542     arbitrarily small norm.
 543
 544     EXAMPLES:
 545
 546     A random element of the trivial cone is zero::
 547
 548         sage: set_random_seed()
 549         sage: K = Cone([], ToricLattice(0))
 550         sage: random_element(K)
 551         ()
 552         sage: K = Cone([(0,)])
 553         sage: random_element(K)
 554         (0)
 555         sage: K = Cone([(0,0)])
 556         sage: random_element(K)
 557         (0, 0)
 558         sage: K = Cone([(0,0,0)])
 559         sage: random_element(K)
 560         (0, 0, 0)
 561
 562     TESTS:
 563
 564     Any cone should contain an element of itself::
 565
 566         sage: set_random_seed()
 567         sage: K = random_cone(max_rays = 8)
 568         sage: K.contains(random_element(K))
 569         True
 570
 571     """
 572     V = K.lattice().vector_space()
 573     F = V.base_ring()
 574     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 575     vector_gens  = map(V, K.rays())
 576     scaled_gens  = [ coefficients[i]*vector_gens[i]
 577                          for i in range(len(vector_gens)) ]
 578
 579     # Make sure we return a vector. Without the coercion, we might
 580     # return ``0`` when ``K`` has no rays.
 581     v = V(sum(scaled_gens))
 582     return v
 583
 584
 585 def positive_operators(K):
 586     r"""
 587     Compute generators of the cone of positive operators on this cone.
 588
 589     OUTPUT:
 590
 591     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 592     Each matrix ``P`` in the list should have the property that ``P*x``
 593     is an element of ``K`` whenever ``x`` is an element of
 594     ``K``. Moreover, any nonnegative linear combination of these
 595     matrices shares the same property.
 596
 597     EXAMPLES:
 598
 599     The trivial cone in a trivial space has no positive operators::
 600
 601         sage: K = Cone([], ToricLattice(0))
 602         sage: positive_operators(K)
 603         []
 604
 605     Positive operators on the nonnegative orthant are nonnegative matrices::
 606
 607         sage: K = Cone([(1,)])
 608         sage: positive_operators(K)
 609         [[1]]
 610
 611         sage: K = Cone([(1,0),(0,1)])
 612         sage: positive_operators(K)
 613         [
 614         [1 0]  [0 1]  [0 0]  [0 0]
 615         [0 0], [0 0], [1 0], [0 1]
 616         ]
 617
 618     Every operator is positive on the ambient vector space::
 619
 620         sage: K = Cone([(1,),(-1,)])
 621         sage: K.is_full_space()
 622         True
 623         sage: positive_operators(K)
 624         [[1], [-1]]
 625
 626         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 627         sage: K.is_full_space()
 628         True
 629         sage: positive_operators(K)
 630         [
 631         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 632         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 633         ]
 634
 635     TESTS:
 636
 637     A positive operator on a cone should send its generators into the cone::
 638
 639         sage: K = random_cone(max_ambient_dim = 6)
 640         sage: pi_of_K = positive_operators(K)
 641         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 642         True
 643
 644     """
 645     # Sage doesn't think matrices are vectors, so we have to convert
 646     # our matrices to vectors explicitly before we can figure out how
 647     # many are linearly-indepenedent.
 648     #
 649     # The space W has the same base ring as V, but dimension
 650     # dim(V)^2. So it has the same dimension as the space of linear
 651     # transformations on V. In other words, it's just the right size
 652     # to create an isomorphism between it and our matrices.
 653     V = K.lattice().vector_space()
 654     W = VectorSpace(V.base_ring(), V.dimension()**2)
 655
 656     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 657
 658     # Turn our matrices into long vectors...
 659     vectors = [ W(m.list()) for m in tensor_products ]
 660
 661     # Create the *dual* cone of the positive operators, expressed as
 662     # long vectors..
 663     L = ToricLattice(W.dimension())
 664     pi_dual = Cone(vectors, lattice=L)
 665
 666     # Now compute the desired cone from its dual...
 667     pi_cone = pi_dual.dual()
 668
 669     # And finally convert its rays back to matrix representations.
 670     M = MatrixSpace(V.base_ring(), V.dimension())
 671
 672     return [ M(v.list()) for v in pi_cone.rays() ]