]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/cone/cone.py
Rename LL() to lyapunov_like_basis().
[sage.d.git] / mjo / cone / cone.py
index 1a3ee3b2fad16e3a05938cb7b9c6a0beae9518d1..cd835a764f9baf7fd8dd2082aa948501610698f5 100644 (file)
@@ -76,11 +76,11 @@ def _basically_the_same(K1, K2):
     if K1.is_strictly_convex() != K2.is_strictly_convex():
         return False
 
-    if len(LL(K1)) != len(LL(K2)):
+    if len(K1.lyapunov_like_basis()) != len(K2.lyapunov_like_basis()):
         return False
 
-    C_of_K1 = discrete_complementarity_set(K1)
-    C_of_K2 = discrete_complementarity_set(K2)
+    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
 
@@ -211,233 +211,6 @@ def _restrict_to_space(K, W):
     return Cone(K_W_rays, lattice=L)
 
 
-
-def discrete_complementarity_set(K):
-    r"""
-    Compute a discrete complementarity set of this cone.
-
-    A discrete complementarity set of `K` is the set of all orthogonal
-    pairs `(x,s)` such that `x \in G_{1}` and `s \in G_{2}` for some
-    generating sets `G_{1}` of `K` and `G_{2}` of its dual. Polyhedral
-    convex cones are input in terms of their generators, so "the" (this
-    particular) discrete complementarity set corresponds to ``G1
-    == K.rays()`` and ``G2 == K.dual().rays()``.
-
-    OUTPUT:
-
-    A list of pairs `(x,s)` such that,
-
-      * Both `x` and `s` are vectors (not rays).
-      * `x` is one of ``K.rays()``.
-      * `s` is one of ``K.dual().rays()``.
-      * `x` and `s` are orthogonal.
-
-    REFERENCES:
-
-    .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
-       Improper Cone. Work in-progress.
-
-    EXAMPLES:
-
-    The discrete complementarity set of the nonnegative orthant consists
-    of pairs of standard basis vectors::
-
-        sage: K = Cone([(1,0),(0,1)])
-        sage: discrete_complementarity_set(K)
-        [((1, 0), (0, 1)), ((0, 1), (1, 0))]
-
-    If the cone consists of a single ray, the second components of the
-    discrete complementarity set should generate the orthogonal
-    complement of that ray::
-
-        sage: K = Cone([(1,0)])
-        sage: discrete_complementarity_set(K)
-        [((1, 0), (0, 1)), ((1, 0), (0, -1))]
-        sage: K = Cone([(1,0,0)])
-        sage: discrete_complementarity_set(K)
-        [((1, 0, 0), (0, 1, 0)),
-          ((1, 0, 0), (0, -1, 0)),
-          ((1, 0, 0), (0, 0, 1)),
-          ((1, 0, 0), (0, 0, -1))]
-
-    When the cone is the entire space, its dual is the trivial cone, so
-    the discrete complementarity set is empty::
-
-        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
-        sage: discrete_complementarity_set(K)
-        []
-
-    Likewise when this cone is trivial (its dual is the entire space)::
-
-        sage: L = ToricLattice(0)
-        sage: K = Cone([], ToricLattice(0))
-        sage: discrete_complementarity_set(K)
-        []
-
-    TESTS:
-
-    The complementarity set of the dual can be obtained by switching the
-    components of the complementarity set of the original cone::
-
-        sage: set_random_seed()
-        sage: K1 = random_cone(max_ambient_dim=6)
-        sage: K2 = K1.dual()
-        sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
-        sage: actual = discrete_complementarity_set(K1)
-        sage: sorted(actual) == sorted(expected)
-        True
-
-    The pairs in the discrete complementarity set are in fact
-    complementary::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=6)
-        sage: dcs = discrete_complementarity_set(K)
-        sage: sum([x.inner_product(s).abs() for (x,s) in dcs])
-        0
-
-    """
-    V = K.lattice().vector_space()
-
-    # Convert rays to vectors so that we can compute inner products.
-    xs = [V(x) for x in K.rays()]
-
-    # We also convert the generators of the dual cone so that we
-    # return pairs of vectors and not (vector, ray) pairs.
-    ss = [V(s) for s in K.dual().rays()]
-
-    return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
-
-
-def LL(K):
-    r"""
-    Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
-    on this cone.
-
-    OUTPUT:
-
-    A list of matrices forming a basis for the space of all
-    Lyapunov-like transformations on the given cone.
-
-    EXAMPLES:
-
-    The trivial cone has no Lyapunov-like transformations::
-
-        sage: L = ToricLattice(0)
-        sage: K = Cone([], lattice=L)
-        sage: LL(K)
-        []
-
-    The Lyapunov-like transformations on the nonnegative orthant are
-    simply diagonal matrices::
-
-        sage: K = Cone([(1,)])
-        sage: LL(K)
-        [[1]]
-
-        sage: K = Cone([(1,0),(0,1)])
-        sage: LL(K)
-        [
-        [1 0]  [0 0]
-        [0 0], [0 1]
-        ]
-
-        sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
-        sage: LL(K)
-        [
-        [1 0 0]  [0 0 0]  [0 0 0]
-        [0 0 0]  [0 1 0]  [0 0 0]
-        [0 0 0], [0 0 0], [0 0 1]
-        ]
-
-    Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
-    `L^{3}_{\infty}` cones [Rudolf et al.]_::
-
-        sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
-        sage: LL(L31)
-        [
-        [1 0 0]
-        [0 1 0]
-        [0 0 1]
-        ]
-
-        sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
-        sage: LL(L3infty)
-        [
-        [1 0 0]
-        [0 1 0]
-        [0 0 1]
-        ]
-
-    If our cone is the entire space, then every transformation on it is
-    Lyapunov-like::
-
-        sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
-        sage: M = MatrixSpace(QQ,2)
-        sage: M.basis() == LL(K)
-        True
-
-    TESTS:
-
-    The inner product `\left< L\left(x\right), s \right>` is zero for
-    every pair `\left( x,s \right)` in the discrete complementarity set
-    of the cone::
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: C_of_K = discrete_complementarity_set(K)
-        sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
-        sage: sum(map(abs, l))
-        0
-
-    The Lyapunov-like transformations on a cone and its dual are related
-    by transposition, but we're not guaranteed to compute transposed
-    elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
-    \right)`
-
-        sage: set_random_seed()
-        sage: K = random_cone(max_ambient_dim=8)
-        sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
-        sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
-        sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
-        sage: LL2_vecs = [ V(m.list()) for m in LL2   ]
-        sage: V.span(LL1_vecs) == V.span(LL2_vecs)
-        True
-
-    """
-    V = K.lattice().vector_space()
-
-    C_of_K = discrete_complementarity_set(K)
-
-    tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ]
-
-    # 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.
-    W = VectorSpace(V.base_ring(), V.dimension()**2)
-
-    # Turn our matrices into long vectors...
-    vectors = [ W(m.list()) for m in tensor_products ]
-
-    # Vector space representation of Lyapunov-like matrices
-    # (i.e. vec(L) where L is Luapunov-like).
-    LL_vector = W.span(vectors).complement()
-
-    # Now construct an ambient MatrixSpace in which to stick our
-    # transformations.
-    M = MatrixSpace(V.base_ring(), V.dimension())
-
-    matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
-
-    return matrix_basis
-
-
-
 def lyapunov_rank(K):
     r"""
     Compute the Lyapunov rank (or bilinearity rank) of this cone.
@@ -641,11 +414,12 @@ def lyapunov_rank(K):
         sage: actual == expected
         True
 
-    The Lyapunov rank of any cone is just the dimension of ``LL(K)``::
+        The Lyapunov rank of any cone is just the dimension of
+        ``K.lyapunov_like_basis()``::
 
         sage: set_random_seed()
         sage: K = random_cone(max_ambient_dim=8)
-        sage: lyapunov_rank(K) == len(LL(K))
+        sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
         True
 
     We can make an imperfect cone perfect by adding a slack variable
@@ -683,7 +457,7 @@ def lyapunov_rank(K):
         # Non-pointed reduction lemma.
         beta += l * m
 
-    beta += len(LL(K))
+    beta += len(K.lyapunov_like_basis())
     return beta
 
 
@@ -741,15 +515,16 @@ def is_lyapunov_like(L,K):
         sage: is_lyapunov_like(L,K)
         True
 
-    Everything in ``LL(K)`` should be Lyapunov-like on ``K``::
+        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 LL(K)])
+        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 discrete_complementarity_set(K)])
+                for (x,s) in K.discrete_complementarity_set()])
 
 
 def random_element(K):
@@ -807,3 +582,194 @@ def random_element(K):
     # return ``0`` when ``K`` has no rays.
     v = V(sum(scaled_gens))
     return v
+
+
+def positive_operators(K):
+    r"""
+    Compute generators of the cone of positive operators on this cone.
+
+    OUTPUT:
+
+    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.
+
+    EXAMPLES:
+
+    The trivial cone in a trivial space has no positive operators::
+
+        sage: K = Cone([], ToricLattice(0))
+        sage: positive_operators(K)
+        []
+
+    Positive operators on the nonnegative orthant are nonnegative matrices::
+
+        sage: K = Cone([(1,)])
+        sage: positive_operators(K)
+        [[1]]
+
+        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]
+        ]
+
+    Every operator is positive on the ambient vector space::
+
+        sage: K = Cone([(1,),(-1,)])
+        sage: K.is_full_space()
+        True
+        sage: positive_operators(K)
+        [[1], [-1]]
+
+        sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
+        sage: K.is_full_space()
+        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:
+
+    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:
+
+    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 ])
+        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::
+
+        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
+        True
+
+    TESTS:
+
+    The Z-property is possessed by every Z-transformation::
+
+        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])
+        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
+        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() ]