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.lyapunov_like_basis()) != len(K2.lyapunov_like_basis()):
  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
 418         ``K.lyapunov_like_basis()``::
 419
 420         sage: set_random_seed()
 421         sage: K = random_cone(max_ambient_dim=8)
 422         sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
 423         True
 424
 425     We can make an imperfect cone perfect by adding a slack variable
 426     (a Theorem in [Orlitzky/Gowda]_)::
 427
 428         sage: set_random_seed()
 429         sage: K = random_cone(max_ambient_dim=8,
 430         ....:                 strictly_convex=True,
 431         ....:                 solid=True)
 432         sage: L = ToricLattice(K.lattice_dim() + 1)
 433         sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
 434         sage: lyapunov_rank(K) >= K.lattice_dim()
 435         True
 436
 437     """
 438     beta = 0
 439
 440     m = K.dim()
 441     n = K.lattice_dim()
 442     l = K.lineality()
 443
 444     if m < n:
 445         # K is not solid, restrict to its span.
 446         K = _restrict_to_space(K, K.span())
 447
 448         # Non-solid reduction lemma.
 449         beta += (n - m)*n
 450
 451     if l > 0:
 452         # K is not pointed, restrict to the span of its dual. Uses a
 453         # proposition from our paper, i.e. this is equivalent to K =
 454         # _rho(K.dual()).dual().
 455         K = _restrict_to_space(K, K.dual().span())
 456
 457         # Non-pointed reduction lemma.
 458         beta += l * m
 459
 460     beta += len(K.lyapunov_like_basis())
 461     return beta
 462
 463
 464
 465 def is_lyapunov_like(L,K):
 466     r"""
 467     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 468
 469     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
 470     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
 471     `\left\langle x,s \right\rangle` in the complementarity set of
 472     ``K``. It is known [Orlitzky]_ that this property need only be
 473     checked for generators of ``K`` and its dual.
 474
 475     INPUT:
 476
 477     - ``L`` -- A linear transformation or matrix.
 478
 479     - ``K`` -- A polyhedral closed convex cone.
 480
 481     OUTPUT:
 482
 483     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
 484     and ``False`` otherwise.
 485
 486     .. WARNING::
 487
 488         If this function returns ``True``, then ``L`` is Lyapunov-like
 489         on ``K``. However, if ``False`` is returned, that could mean one
 490         of two things. The first is that ``L`` is definitely not
 491         Lyapunov-like on ``K``. The second is more of an "I don't know"
 492         answer, returned (for example) if we cannot prove that an inner
 493         product is zero.
 494
 495     REFERENCES:
 496
 497     .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
 498        improper cone (preprint).
 499
 500     EXAMPLES:
 501
 502     The identity is always Lyapunov-like in a nontrivial space::
 503
 504         sage: set_random_seed()
 505         sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
 506         sage: L = identity_matrix(K.lattice_dim())
 507         sage: is_lyapunov_like(L,K)
 508         True
 509
 510     As is the "zero" transformation::
 511
 512         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 513         sage: R = K.lattice().vector_space().base_ring()
 514         sage: L = zero_matrix(R, K.lattice_dim())
 515         sage: is_lyapunov_like(L,K)
 516         True
 517
 518         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 519         on ``K``::
 520
 521         sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
 522         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
 523         True
 524
 525     """
 526     return all([(L*x).inner_product(s) == 0
 527                 for (x,s) in K.discrete_complementarity_set()])
 528
 529
 530 def random_element(K):
 531     r"""
 532     Return a random element of ``K`` from its ambient vector space.
 533
 534     ALGORITHM:
 535
 536     The cone ``K`` is specified in terms of its generators, so that
 537     ``K`` is equal to the convex conic combination of those generators.
 538     To choose a random element of ``K``, we assign random nonnegative
 539     coefficients to each generator of ``K`` and construct a new vector
 540     from the scaled rays.
 541
 542     A vector, rather than a ray, is returned so that the element may
 543     have non-integer coordinates. Thus the element may have an
 544     arbitrarily small norm.
 545
 546     EXAMPLES:
 547
 548     A random element of the trivial cone is zero::
 549
 550         sage: set_random_seed()
 551         sage: K = Cone([], ToricLattice(0))
 552         sage: random_element(K)
 553         ()
 554         sage: K = Cone([(0,)])
 555         sage: random_element(K)
 556         (0)
 557         sage: K = Cone([(0,0)])
 558         sage: random_element(K)
 559         (0, 0)
 560         sage: K = Cone([(0,0,0)])
 561         sage: random_element(K)
 562         (0, 0, 0)
 563
 564     TESTS:
 565
 566     Any cone should contain an element of itself::
 567
 568         sage: set_random_seed()
 569         sage: K = random_cone(max_rays = 8)
 570         sage: K.contains(random_element(K))
 571         True
 572
 573     """
 574     V = K.lattice().vector_space()
 575     F = V.base_ring()
 576     coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
 577     vector_gens  = map(V, K.rays())
 578     scaled_gens  = [ coefficients[i]*vector_gens[i]
 579                          for i in range(len(vector_gens)) ]
 580
 581     # Make sure we return a vector. Without the coercion, we might
 582     # return ``0`` when ``K`` has no rays.
 583     v = V(sum(scaled_gens))
 584     return v
 585
 586
 587 def positive_operators(K):
 588     r"""
 589     Compute generators of the cone of positive operators on this cone.
 590
 591     OUTPUT:
 592
 593     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 594     Each matrix ``P`` in the list should have the property that ``P*x``
 595     is an element of ``K`` whenever ``x`` is an element of
 596     ``K``. Moreover, any nonnegative linear combination of these
 597     matrices shares the same property.
 598
 599     EXAMPLES:
 600
 601     The trivial cone in a trivial space has no positive operators::
 602
 603         sage: K = Cone([], ToricLattice(0))
 604         sage: positive_operators(K)
 605         []
 606
 607     Positive operators on the nonnegative orthant are nonnegative matrices::
 608
 609         sage: K = Cone([(1,)])
 610         sage: positive_operators(K)
 611         [[1]]
 612
 613         sage: K = Cone([(1,0),(0,1)])
 614         sage: positive_operators(K)
 615         [
 616         [1 0]  [0 1]  [0 0]  [0 0]
 617         [0 0], [0 0], [1 0], [0 1]
 618         ]
 619
 620     Every operator is positive on the ambient vector space::
 621
 622         sage: K = Cone([(1,),(-1,)])
 623         sage: K.is_full_space()
 624         True
 625         sage: positive_operators(K)
 626         [[1], [-1]]
 627
 628         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 629         sage: K.is_full_space()
 630         True
 631         sage: positive_operators(K)
 632         [
 633         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 634         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 635         ]
 636
 637     TESTS:
 638
 639     A positive operator on a cone should send its generators into the cone::
 640
 641         sage: K = random_cone(max_ambient_dim = 6)
 642         sage: pi_of_K = positive_operators(K)
 643         sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
 644         True
 645
 646     """
 647     # Sage doesn't think matrices are vectors, so we have to convert
 648     # our matrices to vectors explicitly before we can figure out how
 649     # many are linearly-indepenedent.
 650     #
 651     # The space W has the same base ring as V, but dimension
 652     # dim(V)^2. So it has the same dimension as the space of linear
 653     # transformations on V. In other words, it's just the right size
 654     # to create an isomorphism between it and our matrices.
 655     V = K.lattice().vector_space()
 656     W = VectorSpace(V.base_ring(), V.dimension()**2)
 657
 658     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 659
 660     # Turn our matrices into long vectors...
 661     vectors = [ W(m.list()) for m in tensor_products ]
 662
 663     # Create the *dual* cone of the positive operators, expressed as
 664     # long vectors..
 665     L = ToricLattice(W.dimension())
 666     pi_dual = Cone(vectors, lattice=L)
 667
 668     # Now compute the desired cone from its dual...
 669     pi_cone = pi_dual.dual()
 670
 671     # And finally convert its rays back to matrix representations.
 672     M = MatrixSpace(V.base_ring(), V.dimension())
 673
 674     return [ M(v.list()) for v in pi_cone.rays() ]
 675
 676
 677 def Z_transformations(K):
 678     r"""
 679     Compute generators of the cone of Z-transformations on this cone.
 680
 681     OUTPUT:
 682
 683     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 684     Each matrix ``L`` in the list should have the property that
 685     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 686     discrete complementarity set of ``K``. Moreover, any nonnegative
 687     linear combination of these matrices shares the same property.
 688
 689     EXAMPLES:
 690
 691     Z-transformations on the nonnegative orthant are just Z-matrices.
 692     That is, matrices whose off-diagonal elements are nonnegative::
 693
 694         sage: K = Cone([(1,0),(0,1)])
 695         sage: Z_transformations(K)
 696         [
 697         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 698         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 699         ]
 700         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 701         sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
 702         ....:                    for i in range(z.nrows())
 703         ....:                    for j in range(z.ncols())
 704         ....:                    if i != j ])
 705         True
 706
 707     The trivial cone in a trivial space has no Z-transformations::
 708
 709         sage: K = Cone([], ToricLattice(0))
 710         sage: Z_transformations(K)
 711         []
 712
 713     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 714
 715         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 716         sage: K.is_full_space()
 717         True
 718         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 719         sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
 720         sage: zs == lls
 721         True
 722
 723     TESTS:
 724
 725     The Z-property is possessed by every Z-transformation::
 726
 727         sage: set_random_seed()
 728         sage: K = random_cone(max_ambient_dim = 6)
 729         sage: Z_of_K = Z_transformations(K)
 730         sage: dcs = K.discrete_complementarity_set()
 731         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 732         ....:                                  for (x,s) in dcs])
 733         True
 734
 735     The lineality space of Z is LL::
 736
 737         sage: set_random_seed()
 738         sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
 739         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 740         sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
 741         sage: z_cone.linear_subspace() == lls
 742         True
 743
 744     """
 745     # Sage doesn't think matrices are vectors, so we have to convert
 746     # our matrices to vectors explicitly before we can figure out how
 747     # many are linearly-indepenedent.
 748     #
 749     # The space W has the same base ring as V, but dimension
 750     # dim(V)^2. So it has the same dimension as the space of linear
 751     # transformations on V. In other words, it's just the right size
 752     # to create an isomorphism between it and our matrices.
 753     V = K.lattice().vector_space()
 754     W = VectorSpace(V.base_ring(), V.dimension()**2)
 755
 756     C_of_K = K.discrete_complementarity_set()
 757     tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
 758
 759     # Turn our matrices into long vectors...
 760     vectors = [ W(m.list()) for m in tensor_products ]
 761
 762     # Create the *dual* cone of the cross-positive operators,
 763     # expressed as long vectors..
 764     L = ToricLattice(W.dimension())
 765     Sigma_dual = Cone(vectors, lattice=L)
 766
 767     # Now compute the desired cone from its dual...
 768     Sigma_cone = Sigma_dual.dual()
 769
 770     # And finally convert its rays back to matrix representations.
 771     # But first, make them negative, so we get Z-transformations and
 772     # not cross-positive ones.
 773     M = MatrixSpace(V.base_ring(), V.dimension())
 774
 775     return [ -M(v.list()) for v in Sigma_cone.rays() ]