Fix the is_positive_on test and give better examples.

[sage.d.git] / mjo / cone / cone.py

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index c1b7ffdcbdc4e6a8e585b38e5493fe1d4060f3c8..01ef5f9d5cc56b087344237b9b5c915423ac0921 100644 (file)
--- a/mjo/cone/cone.py
+++ b/mjo/cone/cone.py
@@ -1,471 +1,172 @@
-# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
-# have to explicitly mangle our sitedir here so that "mjo.cone"
-# resolves.
-from os.path import abspath
-from site import addsitedir
-addsitedir(abspath('../../'))
-

 from sage.all import *
 
 from sage.all import *
 

-
-def _basically_the_same(K1, K2):
+def is_positive_on(L,K):

     r"""
     r"""

-    Test whether or not ``K1`` and ``K2`` are "basically the same."
-
-    This is a hack to get around the fact that it's difficult to tell
-    when two cones are linearly isomorphic. We have a proposition that
-    equates two cones, but represented over `\mathbb{Q}`, they are
-    merely linearly isomorphic (not equal). So rather than test for
-    equality, we test a list of properties that should be preserved
-    under an invertible linear transformation.
-
-    OUTPUT:
-
-    ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
-    otherwise.
-
-    EXAMPLES:
-
-    Any proper cone with three generators in `\mathbb{R}^{3}` is
-    basically the same as the nonnegative orthant::
-
-        sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
-        sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
-        sage: _basically_the_same(K1, K2)
-        True
-
-    Negating a cone gives you another cone that is basically the same::
-
-        sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
-        sage: _basically_the_same(K, -K)
-        True
-
-    TESTS:
-
-    Any cone is basically the same as itself::
-
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: _basically_the_same(K, K)
-        True
-
-    After applying an invertible matrix to the rows of a cone, the
-    result should be basically the same as the cone we started with::
-
-        sage: K1 = random_cone(max_ambient_dim = 8)
-        sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
-        sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
-        sage: _basically_the_same(K1, K2)
-        True
-
-    """
-    if K1.lattice_dim() != K2.lattice_dim():
-        return False
-
-    if K1.nrays() != K2.nrays():
-        return False
-
-    if K1.dim() != K2.dim():
-        return False
-
-    if K1.lineality() != K2.lineality():
-        return False
+    Determine whether or not ``L`` is positive on ``K``.

 
 

-    if K1.is_solid() != K2.is_solid():
-        return False
-
-    if K1.is_strictly_convex() != K2.is_strictly_convex():
-        return False
-
-    if len(K1.lyapunov_like_basis()) != len(K2.lyapunov_like_basis()):
-        return False
-
-    C_of_K1 = K1.discrete_complementarity_set()
-    C_of_K2 = K2.discrete_complementarity_set()
-    if len(C_of_K1) != len(C_of_K2):
-        return False
-
-    if len(K1.facets()) != len(K2.facets()):
-        return False
-
-    return True
+    We say that ``L`` is positive on ``K`` if `L\left\lparen x
+    \right\rparen` belongs to ``K`` for all `x` in ``K``. This
+    property need only be checked for generators of ``K``.

 
 

+    INPUT:

 
 

+    - ``L`` -- A linear transformation or matrix.

 
 

-def _restrict_to_space(K, W):
-    r"""
-    Restrict this cone a subspace of its ambient space.
+    - ``K`` -- A polyhedral closed convex cone.

 
 

-    INPUT:
+    OUTPUT:

 
 

-    - ``W`` -- The subspace into which this cone will be restricted.
+    ``True`` if it can be proven that ``L`` is positive on ``K``,
+    and ``False`` otherwise.

 
 

-    OUTPUT:
+    .. WARNING::

 
 

-    A new cone in a sublattice corresponding to ``W``.
+        If this function returns ``True``, then ``L`` is positive
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        positive on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is nonnegative.

 
     EXAMPLES:
 
 
     EXAMPLES:
 

-    When this cone is solid, restricting it into its own span should do
-    nothing::
+    Positive operators on the nonnegative orthant are nonnegative
+    matrices::

 
 

-        sage: K = Cone([(1,)])
-        sage: _restrict_to_space(K, K.span()) == K
-        True
-
-    A single ray restricted into its own span gives the same output
-    regardless of the ambient space::
-
-        sage: K2 = Cone([(1,0)])
-        sage: K2_S = _restrict_to_space(K2, K2.span()).rays()
-        sage: K2_S
-        N(1)
-        in 1-d lattice N
-        sage: K3 = Cone([(1,0,0)])
-        sage: K3_S = _restrict_to_space(K3, K3.span()).rays()
-        sage: K3_S
-        N(1)
-        in 1-d lattice N
-        sage: K2_S == K3_S
+        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
+        sage: L = random_matrix(QQ,3).apply_map(abs)
+        sage: is_positive_on(L,K)

         True
 
     TESTS:
 
         True
 
     TESTS:
 

-    The projected cone should always be solid::
+    The identity is always positive in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: _restrict_to_space(K, K.span()).is_solid()
-        True
-
-    And the resulting cone should live in a space having the same
-    dimension as the space we restricted it to::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: K_P = _restrict_to_space(K, K.dual().span())
-        sage: K_P.lattice_dim() == K.dual().dim()
-        True
-
-    This function should not affect the dimension of a cone::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: K.dim() == _restrict_to_space(K,K.span()).dim()
-        True
-
-    Nor should it affect the lineality of a cone::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: K.lineality() == _restrict_to_space(K, K.span()).lineality()
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_positive_on(L,K)

         True
 
         True
 

-    No matter which space we restrict to, the lineality should not
-    increase::
+    As is the "zero" transformation::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: S = K.span(); P = K.dual().span()
-        sage: K.lineality() >= _restrict_to_space(K,S).lineality()
-        True
-        sage: K.lineality() >= _restrict_to_space(K,P).lineality()
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_positive_on(L,K)

         True
 
         True
 

-    If we do this according to our paper, then the result is proper::
+    Everything in ``K.positive_operators_gens()`` should be
+    positive on ``K``::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 8)
-        sage: K_S = _restrict_to_space(K, K.span())
-        sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
-        sage: K_SP.is_proper()
-        True
-        sage: K_SP = _restrict_to_space(K_S, K_S.dual().span())
-        sage: K_SP.is_proper()
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_positive_on(L,K)
+        ....:       for L in K.positive_operators_gens() ])

         True
         True

-
-    Test the proposition in our paper concerning the duals and
-    restrictions. Generate a random cone, then create a subcone of
-    it. The operation of dual-taking should then commute with
-    _restrict_to_space::
-
-        sage: set_random_seed()
-        sage: J = random_cone(max_ambient_dim = 8)
-        sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
-        sage: K_W_star = _restrict_to_space(K, J.span()).dual()
-        sage: K_star_W = _restrict_to_space(K.dual(), J.span())
-        sage: _basically_the_same(K_W_star, K_star_W)
+        sage: all([ is_positive_on(L.change_ring(SR),K)
+        ....:       for L in K.positive_operators_gens() ])

         True
 
     """
         True
 
     """

-    # First we want to intersect ``K`` with ``W``. The easiest way to
-    # do this is via cone intersection, so we turn the subspace ``W``
-    # into a cone.
-    W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice())
-    K = K.intersection(W_cone)
-
-    # We've already intersected K with the span of K2, so every
-    # generator of K should belong to W now.
-    K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ]
-
-    L = ToricLattice(W.dimension())
-    return Cone(K_W_rays, lattice=L)
-
-
-def lyapunov_rank(K):
+    if L.base_ring().is_exact():
+        # This could potentially be extended to other types of ``K``...
+        return all([ L*x in K for x in K ])
+    elif L.base_ring() is SR:
+        # Fall back to inequality-checking when the entries of ``L``
+        # might be symbolic.
+        return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')
+
+
+def is_cross_positive_on(L,K):

     r"""
     r"""

-    Compute the Lyapunov rank of this cone.
+    Determine whether or not ``L`` is cross-positive on ``K``.

 
 

-    The Lyapunov rank of a cone is the dimension of the space of its
-    Lyapunov-like transformations -- that is, the length of a
-    :meth:`lyapunov_like_basis`. Equivalently, the Lyapunov rank is the
-    dimension of the Lie algebra of the automorphism group of the cone.
-
-    OUTPUT:
-
-    A nonnegative integer representing the Lyapunov rank of this cone.
+    We say that ``L`` is cross-positive on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle \ge 0` for all pairs
+    `\left\langle x,s \right\rangle` in the complementarity set of
+    ``K``. This property need only be checked for generators of
+    ``K`` and its dual.

 
 

-    If the ambient space is trivial, the Lyapunov rank will be zero.
-    Otherwise, if the dimension of the ambient vector space is `n`, then
-    the resulting Lyapunov rank will be between `1` and `n` inclusive. A
-    Lyapunov rank of `n-1` is not possible [Orlitzky]_.
+    INPUT:

 
 

-    ALGORITHM:
+    - ``L`` -- A linear transformation or matrix.

 
 

-    The codimension formula from the second reference is used. We find
-    all pairs `(x,s)` in the complementarity set of `K` such that `x`
-    and `s` are rays of our cone. It is known that these vectors are
-    sufficient to apply the codimension formula. Once we have all such
-    pairs, we "brute force" the codimension formula by finding all
-    linearly-independent `xs^{T}`.
+    - ``K`` -- A polyhedral closed convex cone.

 
 

-    REFERENCES:
+    OUTPUT:

 
 

-    .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of
-       a proper cone and Lyapunov-like transformations. Mathematical
-       Programming, 147 (2014) 155-170.
+    ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
+    and ``False`` otherwise.

 
 

-    M. Orlitzky. The Lyapunov rank of an improper cone.
-    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
+    .. WARNING::

 
 

-    G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
-    optimality constraints for the cone of positive polynomials,
-    Mathematical Programming, Series B, 129 (2011) 5-31.
+        If this function returns ``True``, then ``L`` is cross-positive
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        cross-positive on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is nonnegative.

 
     EXAMPLES:
 
 
     EXAMPLES:
 

-    The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
-    [Rudolf]_::
-
-        sage: positives = Cone([(1,)])
-        sage: lyapunov_rank(positives)
-        1
-        sage: quadrant = Cone([(1,0), (0,1)])
-        sage: lyapunov_rank(quadrant)
-        2
-       sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
-        sage: lyapunov_rank(octant)
-        3
-
-    The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
-    [Orlitzky]_::
-
-        sage: R5 = VectorSpace(QQ, 5)
-        sage: gs = R5.basis() + [ -r for r in R5.basis() ]
-        sage: K = Cone(gs)
-        sage: lyapunov_rank(K)
-        25
-
-    The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
-    [Rudolf]_::
-
-        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
-        sage: lyapunov_rank(L31)
-        1
-
-    Likewise for the `L^{3}_{\infty}` cone [Rudolf]_::
-
-        sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
-        sage: lyapunov_rank(L3infty)
-        1
-
-    A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
-    + 1` [Orlitzky]_::
-
-        sage: K = Cone([(1,0,0,0,0)])
-        sage: lyapunov_rank(K)
-        21
-        sage: K.lattice_dim()**2 - K.lattice_dim() + 1
-        21
-
-    A subspace (of dimension `m`) in `n` dimensions should have a
-    Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky]_::
-
-        sage: e1 = (1,0,0,0,0)
-        sage: neg_e1 = (-1,0,0,0,0)
-        sage: e2 = (0,1,0,0,0)
-        sage: neg_e2 = (0,-1,0,0,0)
-        sage: z = (0,0,0,0,0)
-        sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
-        sage: lyapunov_rank(K)
-        19
-        sage: K.lattice_dim()**2 - K.dim()*K.codim()
-        19
-
-    The Lyapunov rank should be additive on a product of proper cones
-    [Rudolf]_::
-
-        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
-        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
-        sage: K = L31.cartesian_product(octant)
-        sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
-        True
-
-    Two isomorphic cones should have the same Lyapunov rank [Rudolf]_.
-    The cone ``K`` in the following example is isomorphic to the nonnegative
-    octant in `\mathbb{R}^{3}`::
-
-        sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
-        sage: lyapunov_rank(K)
-        3
-
-    The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
-    itself [Rudolf]_::
-
-        sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
-        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
-        True
-
-    TESTS:
-
-    The Lyapunov rank should be additive on a product of proper cones
-    [Rudolf]_::
+    The identity is always cross-positive in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K1 = random_cone(max_ambient_dim=8,
-        ....:                  strictly_convex=True,
-        ....:                  solid=True)
-        sage: K2 = random_cone(max_ambient_dim=8,
-        ....:                  strictly_convex=True,
-        ....:                  solid=True)
-        sage: K = K1.cartesian_product(K2)
-        sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
-        True
-
-    The Lyapunov rank is invariant under a linear isomorphism
-    [Orlitzky]_::
-
-        sage: K1 = random_cone(max_ambient_dim = 8)
-        sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
-        sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
-        sage: lyapunov_rank(K1) == lyapunov_rank(K2)
-        True
-
-    The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
-    itself [Rudolf]_::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_cross_positive_on(L,K)

         True
 
         True
 

-    The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
-    be any number between `1` and `n` inclusive, excluding `n-1`
-    [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
-    trivial cone in a trivial space as well. However, in zero dimensions,
-    the Lyapunov rank of the trivial cone will be zero::
+    As is the "zero" transformation::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8,
-        ....:                 strictly_convex=True,
-        ....:                 solid=True)
-        sage: b = lyapunov_rank(K)
-        sage: n = K.lattice_dim()
-        sage: (n == 0 or 1 <= b) and b <= n
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_cross_positive_on(L,K)

         True
         True

-        sage: b == n-1
-        False

 
 

-    In fact [Orlitzky]_, no closed convex polyhedral cone can have
-    Lyapunov rank `n-1` in `n` dimensions::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: b = lyapunov_rank(K)
-        sage: n = K.lattice_dim()
-        sage: b == n-1
-        False
-
-    The calculation of the Lyapunov rank of an improper cone can be
-    reduced to that of a proper cone [Orlitzky]_::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: actual = lyapunov_rank(K)
-        sage: K_S = _restrict_to_space(K, K.span())
-        sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual()
-        sage: l = K.lineality()
-        sage: c = K.codim()
-        sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
-        sage: actual == expected
-        True
+    TESTS:

 
 

-    The Lyapunov rank of a cone is the size of a :meth:`lyapunov_like_basis`::
+    Everything in ``K.cross_positive_operators_gens()`` should be
+    cross-positive on ``K``::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_cross_positive_on(L,K)
+        ....:       for L in K.cross_positive_operators_gens() ])

         True
         True

-
-    We can make an imperfect cone perfect by adding a slack variable
-    (a Theorem in [Orlitzky]_)::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8,
-        ....:                 strictly_convex=True,
-        ....:                 solid=True)
-        sage: L = ToricLattice(K.lattice_dim() + 1)
-        sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
-        sage: lyapunov_rank(K) >= K.lattice_dim()
+        sage: all([ is_cross_positive_on(L.change_ring(SR),K)
+        ....:       for L in K.cross_positive_operators_gens() ])

         True
 
     """
         True
 
     """

-    beta = 0 # running tally of the Lyapunov rank
+    if L.base_ring().is_exact() or L.base_ring() is SR:
+        return all([ s*(L*x) >= 0
+                     for (x,s) in K.discrete_complementarity_set() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')

 
 

-    m = K.dim()
-    n = K.lattice_dim()
-    l = K.lineality()

 
 

-    if m < n:
-        # K is not solid, restrict to its span.
-        K = _restrict_to_space(K, K.span())
-
-        # Non-solid reduction lemma.
-        beta += (n - m)*n
-
-    if l > 0:
-        # K is not pointed, restrict to the span of its dual. Uses a
-        # proposition from our paper, i.e. this is equivalent to K =
-        # _rho(K.dual()).dual().
-        K = _restrict_to_space(K, K.dual().span())
-
-        # Non-pointed reduction lemma.
-        beta += l * m
-
-    beta += len(K.lyapunov_like_basis())
-    return beta
-
-
-
-def is_lyapunov_like(L,K):
+def is_Z_on(L,K):

     r"""
     r"""

-    Determine whether or not ``L`` is Lyapunov-like on ``K``.
+    Determine whether or not ``L`` is a Z-operator on ``K``.

 
 

-    We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
-    L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
+    We say that ``L`` is a Z-operator on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle \le 0` for all pairs

     `\left\langle x,s \right\rangle` in the complementarity set of
     `\left\langle x,s \right\rangle` in the complementarity set of

-    ``K``. It is known [Orlitzky]_ that this property need only be
+    ``K``. It is known that this property need only be

     checked for generators of ``K`` and its dual.
 
     checked for generators of ``K`` and its dual.
 

+    A matrix is a Z-operator on ``K`` if and only if its negation is a
+    cross-positive operator on ``K``.
+

     INPUT:
 
     - ``L`` -- A linear transformation or matrix.
     INPUT:
 
     - ``L`` -- A linear transformation or matrix.


@@ -474,296 +175,152 @@ def is_lyapunov_like(L,K):
 
     OUTPUT:
 
 
     OUTPUT:
 

-    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
+    ``True`` if it can be proven that ``L`` is a Z-operator on ``K``,

     and ``False`` otherwise.
 
     .. WARNING::
 
     and ``False`` otherwise.
 
     .. WARNING::
 

-        If this function returns ``True``, then ``L`` is Lyapunov-like
+        If this function returns ``True``, then ``L`` is a Z-operator

         on ``K``. However, if ``False`` is returned, that could mean one
         of two things. The first is that ``L`` is definitely not
         on ``K``. However, if ``False`` is returned, that could mean one
         of two things. The first is that ``L`` is definitely not

-        Lyapunov-like on ``K``. The second is more of an "I don't know"
+        a Z-operator on ``K``. The second is more of an "I don't know"

         answer, returned (for example) if we cannot prove that an inner
         answer, returned (for example) if we cannot prove that an inner

-        product is zero.
-
-    REFERENCES:
-
-    M. Orlitzky. The Lyapunov rank of an improper cone.
-    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
+        product is nonnegative.

 
     EXAMPLES:
 
 
     EXAMPLES:
 

-    The identity is always Lyapunov-like in a nontrivial space::
+    The identity is always a Z-operator in a nontrivial space::

 
         sage: set_random_seed()
 
         sage: set_random_seed()

-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 8)
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)

         sage: L = identity_matrix(K.lattice_dim())
         sage: L = identity_matrix(K.lattice_dim())

-        sage: is_lyapunov_like(L,K)
+        sage: is_Z_on(L,K)

         True
 
     As is the "zero" transformation::
 
         True
 
     As is the "zero" transformation::
 

-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)

         sage: R = K.lattice().vector_space().base_ring()
         sage: L = zero_matrix(R, K.lattice_dim())
         sage: R = K.lattice().vector_space().base_ring()
         sage: L = zero_matrix(R, K.lattice_dim())

-        sage: is_lyapunov_like(L,K)
+        sage: is_Z_on(L,K)

         True
 
         True
 

-        Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
-        on ``K``::
-
-        sage: K = random_cone(min_ambient_dim = 1, max_rays = 5)
-        sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
-        True
-
-    """
-    return all([(L*x).inner_product(s) == 0
-                for (x,s) in K.discrete_complementarity_set()])
-
-
-def random_element(K):
-    r"""
-    Return a random element of ``K`` from its ambient vector space.
-
-    ALGORITHM:
-
-    The cone ``K`` is specified in terms of its generators, so that
-    ``K`` is equal to the convex conic combination of those generators.
-    To choose a random element of ``K``, we assign random nonnegative
-    coefficients to each generator of ``K`` and construct a new vector
-    from the scaled rays.
-
-    A vector, rather than a ray, is returned so that the element may
-    have non-integer coordinates. Thus the element may have an
-    arbitrarily small norm.
-
-    EXAMPLES:
-
-    A random element of the trivial cone is zero::
-
-        sage: set_random_seed()
-        sage: K = Cone([], ToricLattice(0))
-        sage: random_element(K)
-        ()
-        sage: K = Cone([(0,)])
-        sage: random_element(K)
-        (0)
-        sage: K = Cone([(0,0)])
-        sage: random_element(K)
-        (0, 0)
-        sage: K = Cone([(0,0,0)])
-        sage: random_element(K)
-        (0, 0, 0)
-

     TESTS:
 
     TESTS:
 

-    Any cone should contain an element of itself::
+    Everything in ``K.Z_operators_gens()`` should be a Z-operator
+    on ``K``::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_rays = 8)
-        sage: K.contains(random_element(K))
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_Z_on(L,K)
+        ....:       for L in K.Z_operators_gens() ])
+        True
+        sage: all([ is_Z_on(L.change_ring(SR),K)
+        ....:       for L in K.Z_operators_gens() ])

         True
 
     """
         True
 
     """

-    V = K.lattice().vector_space()
-    F = V.base_ring()
-    coefficients = [ F.random_element().abs() for i in range(K.nrays()) ]
-    vector_gens  = map(V, K.rays())
-    scaled_gens  = [ coefficients[i]*vector_gens[i]
-                         for i in range(len(vector_gens)) ]
+    return is_cross_positive_on(-L,K)

 
 

-    # Make sure we return a vector. Without the coercion, we might
-    # return ``0`` when ``K`` has no rays.
-    v = V(sum(scaled_gens))
-    return v

 
 

-
-def positive_operators(K):
+def is_lyapunov_like_on(L,K):

     r"""
     r"""

-    Compute generators of the cone of positive operators on this cone.
+    Determine whether or not ``L`` is Lyapunov-like on ``K``.

 
 

-    OUTPUT:
+    We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
+    L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
+    `\left\langle x,s \right\rangle` in the complementarity set of
+    ``K``. This property need only be checked for generators of
+    ``K`` and its dual.

 
 

-    A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
-    Each matrix ``P`` in the list should have the property that ``P*x``
-    is an element of ``K`` whenever ``x`` is an element of
-    ``K``. Moreover, any nonnegative linear combination of these
-    matrices shares the same property.
+    INPUT:

 
 

-    EXAMPLES:
+    - ``L`` -- A linear transformation or matrix.

 
 

-    The trivial cone in a trivial space has no positive operators::
+    - ``K`` -- A polyhedral closed convex cone.

 
 

-        sage: K = Cone([], ToricLattice(0))
-        sage: positive_operators(K)
-        []
+    OUTPUT:

 
 

-    Positive operators on the nonnegative orthant are nonnegative matrices::
+    ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
+    and ``False`` otherwise.

 
 

-        sage: K = Cone([(1,)])
-        sage: positive_operators(K)
-        [[1]]
+    .. WARNING::

 
 

-        sage: K = Cone([(1,0),(0,1)])
-        sage: positive_operators(K)
-        [
-        [1 0]  [0 1]  [0 0]  [0 0]
-        [0 0], [0 0], [1 0], [0 1]
-        ]
+        If this function returns ``True``, then ``L`` is Lyapunov-like
+        on ``K``. However, if ``False`` is returned, that could mean one
+        of two things. The first is that ``L`` is definitely not
+        Lyapunov-like on ``K``. The second is more of an "I don't know"
+        answer, returned (for example) if we cannot prove that an inner
+        product is zero.

 
 

-    Every operator is positive on the ambient vector space::
+    EXAMPLES:

 
 

-        sage: K = Cone([(1,),(-1,)])
-        sage: K.is_full_space()
-        True
-        sage: positive_operators(K)
-        [[1], [-1]]
+    Lyapunov-like operators on the nonnegative orthant are diagonal
+    matrices::

 
 

-        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
-        sage: K.is_full_space()
+        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
+        sage: L = diagonal_matrix(random_vector(QQ,3))
+        sage: is_lyapunov_like_on(L,K)

         True
         True

-        sage: positive_operators(K)
-        [
-        [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
-        [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
-        ]

 
     TESTS:
 
 
     TESTS:
 

-    A positive operator on a cone should send its generators into the cone::
-
-        sage: K = random_cone(max_ambient_dim = 6)
-        sage: pi_of_K = positive_operators(K)
-        sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()])
-        True
-
-    """
-    # Sage doesn't think matrices are vectors, so we have to convert
-    # our matrices to vectors explicitly before we can figure out how
-    # many are linearly-indepenedent.
-    #
-    # The space W has the same base ring as V, but dimension
-    # dim(V)^2. So it has the same dimension as the space of linear
-    # transformations on V. In other words, it's just the right size
-    # to create an isomorphism between it and our matrices.
-    V = K.lattice().vector_space()
-    W = VectorSpace(V.base_ring(), V.dimension()**2)
-
-    tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
-
-    # Turn our matrices into long vectors...
-    vectors = [ W(m.list()) for m in tensor_products ]
-
-    # Create the *dual* cone of the positive operators, expressed as
-    # long vectors..
-    L = ToricLattice(W.dimension())
-    pi_dual = Cone(vectors, lattice=L)
-
-    # Now compute the desired cone from its dual...
-    pi_cone = pi_dual.dual()
-
-    # And finally convert its rays back to matrix representations.
-    M = MatrixSpace(V.base_ring(), V.dimension())
-
-    return [ M(v.list()) for v in pi_cone.rays() ]
-
-
-def Z_transformations(K):
-    r"""
-    Compute generators of the cone of Z-transformations on this cone.
-
-    OUTPUT:
-
-    A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
-    Each matrix ``L`` in the list should have the property that
-    ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
-    discrete complementarity set of ``K``. Moreover, any nonnegative
-    linear combination of these matrices shares the same property.
-
-    EXAMPLES:
+    The identity is always Lyapunov-like in a nontrivial space::

 
 

-    Z-transformations on the nonnegative orthant are just Z-matrices.
-    That is, matrices whose off-diagonal elements are nonnegative::
-
-        sage: K = Cone([(1,0),(0,1)])
-        sage: Z_transformations(K)
-        [
-        [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
-        [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
-        ]
-        sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
-        sage: all([ z[i][j] <= 0 for z in Z_transformations(K)
-        ....:                    for i in range(z.nrows())
-        ....:                    for j in range(z.ncols())
-        ....:                    if i != j ])
+        sage: set_random_seed()
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: L = identity_matrix(K.lattice_dim())
+        sage: is_lyapunov_like_on(L,K)

         True
 
         True
 

-    The trivial cone in a trivial space has no Z-transformations::
-
-        sage: K = Cone([], ToricLattice(0))
-        sage: Z_transformations(K)
-        []
-
-    Z-transformations on a subspace are Lyapunov-like and vice-versa::
+    As is the "zero" transformation::

 
 

-        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
-        sage: K.is_full_space()
-        True
-        sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
-        sage: zs  = span([ vector(z.list()) for z in Z_transformations(K) ])
-        sage: zs == lls
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
+        sage: R = K.lattice().vector_space().base_ring()
+        sage: L = zero_matrix(R, K.lattice_dim())
+        sage: is_lyapunov_like_on(L,K)

         True
 
         True
 

-    TESTS:
-
-    The Z-property is possessed by every Z-transformation::
+    Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
+    on ``K``::

 
 

-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim = 6)
-        sage: Z_of_K = Z_transformations(K)
-        sage: dcs = K.discrete_complementarity_set()
-        sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
-        ....:                                  for (x,s) in dcs])
+        sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
+        sage: all([ is_lyapunov_like_on(L,K)
+        ....:       for L in K.lyapunov_like_basis() ])

         True
         True

-
-    The lineality space of Z is LL::
-
-        sage: set_random_seed()
-        sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6)
-        sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
-        sage: z_cone  = Cone([ z.list() for z in Z_transformations(K) ])
-        sage: z_cone.linear_subspace() == lls
+        sage: all([ is_lyapunov_like_on(L.change_ring(SR),K)
+        ....:       for L in K.lyapunov_like_basis() ])

         True
 
     """
         True
 
     """

-    # Sage doesn't think matrices are vectors, so we have to convert
-    # our matrices to vectors explicitly before we can figure out how
-    # many are linearly-indepenedent.
-    #
-    # The space W has the same base ring as V, but dimension
-    # dim(V)^2. So it has the same dimension as the space of linear
-    # transformations on V. In other words, it's just the right size
-    # to create an isomorphism between it and our matrices.
-    V = K.lattice().vector_space()
-    W = VectorSpace(V.base_ring(), V.dimension()**2)
-
-    C_of_K = K.discrete_complementarity_set()
-    tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
-
-    # Turn our matrices into long vectors...
-    vectors = [ W(m.list()) for m in tensor_products ]
-
-    # Create the *dual* cone of the cross-positive operators,
-    # expressed as long vectors..
-    L = ToricLattice(W.dimension())
-    Sigma_dual = Cone(vectors, lattice=L)
-
-    # Now compute the desired cone from its dual...
-    Sigma_cone = Sigma_dual.dual()
-
-    # And finally convert its rays back to matrix representations.
-    # But first, make them negative, so we get Z-transformations and
-    # not cross-positive ones.
-    M = MatrixSpace(V.base_ring(), V.dimension())
-
-    return [ -M(v.list()) for v in Sigma_cone.rays() ]
+    if L.base_ring().is_exact() or L.base_ring() is SR:
+        # The "fast method" of creating a vector space based on a
+        # ``lyapunov_like_basis`` is actually slower than this.
+        return all([ s*(L*x) == 0
+                     for (x,s) in K.discrete_complementarity_set() ])
+    else:
+        # The only inexact ring that we're willing to work with is SR,
+        # since it can still be exact when working with symbolic
+        # constants like pi and e.
+        raise ValueError('base ring of operator L is neither SR nor exact')
+
+def LL_cone(K):
+    gens = K.lyapunov_like_basis()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def Sigma_cone(K):
+    gens = K.cross_positive_operators_gens()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def Z_cone(K):
+    gens = K.Z_operators_gens()
+    L = ToricLattice(K.lattice_dim()**2)
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)
+
+def pi_cone(K1, K2=None):
+    if K2 is None:
+        K2 = K1
+    gens = K1.positive_operators_gens(K2)
+    L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
+    return Cone([ g.list() for g in gens ], lattice=L, check=False)