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