]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/cone/cone.py
Remove _basically_the_same from cone.py.
[sage.d.git] / mjo / cone / cone.py
index 3f915ed66fb4ea33629bc1c30ff4e021f59f8607..5f058ca9a8c055f972ffffaa216975517e0dc723 100644 (file)
@@ -8,89 +8,6 @@ addsitedir(abspath('../../'))
 from sage.all import *
 
 
-def _basically_the_same(K1, K2):
-    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
-
-    if K1.is_solid() != K2.is_solid():
-        return False
-
-    if K1.is_strictly_convex() != K2.is_strictly_convex():
-        return False
-
-    if len(K1.LL()) != len(K2.LL()):
-        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
-
-
-
 def _restrict_to_space(K, W):
     r"""
     Restrict this cone a subspace of its ambient space.
@@ -213,26 +130,21 @@ def _restrict_to_space(K, W):
 
 def lyapunov_rank(K):
     r"""
-    Compute the Lyapunov rank (or bilinearity rank) of this cone.
-
-    The Lyapunov rank of a cone can be thought of in (mainly) two ways:
+    Compute the Lyapunov rank of this cone.
 
-    1. The dimension of the Lie algebra of the automorphism group of the
-       cone.
-
-    2. The dimension of the linear space of all Lyapunov-like
-       transformations on the cone.
-
-    INPUT:
-
-    A closed, convex polyhedral cone.
+    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:
 
-    An integer representing the Lyapunov rank of the cone. If the
-    dimension of the ambient vector space is `n`, then the Lyapunov rank
-    will be between `1` and `n` inclusive; however a rank of `n-1` is
-    not possible (see [Orlitzky/Gowda]_).
+    A nonnegative integer representing the Lyapunov rank of this cone.
+
+    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]_.
 
     ALGORITHM:
 
@@ -245,21 +157,21 @@ def lyapunov_rank(K):
 
     REFERENCES:
 
-    .. [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.
+    .. [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.
 
-    .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
-       Improper Cone. Work in-progress.
+    M. Orlitzky. The Lyapunov rank of an improper cone.
+    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 
-    .. [Rudolf et al.] 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.
+    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.
 
     EXAMPLES:
 
     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
-    [Rudolf et al.]_::
+    [Rudolf]_::
 
         sage: positives = Cone([(1,)])
         sage: lyapunov_rank(positives)
@@ -272,7 +184,7 @@ def lyapunov_rank(K):
         3
 
     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
-    [Orlitzky/Gowda]_::
+    [Orlitzky]_::
 
         sage: R5 = VectorSpace(QQ, 5)
         sage: gs = R5.basis() + [ -r for r in R5.basis() ]
@@ -281,20 +193,20 @@ def lyapunov_rank(K):
         25
 
     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
-    [Rudolf et al.]_::
+    [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 et al.]_::
+    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/Gowda]_::
+    + 1` [Orlitzky]_::
 
         sage: K = Cone([(1,0,0,0,0)])
         sage: lyapunov_rank(K)
@@ -303,7 +215,7 @@ def lyapunov_rank(K):
         21
 
     A subspace (of dimension `m`) in `n` dimensions should have a
-    Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
+    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)
@@ -317,7 +229,7 @@ def lyapunov_rank(K):
         19
 
     The Lyapunov rank should be additive on a product of proper cones
-    [Rudolf et al.]_::
+    [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)])
@@ -325,7 +237,7 @@ def lyapunov_rank(K):
         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
         True
 
-    Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
+    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}`::
 
@@ -334,7 +246,7 @@ def lyapunov_rank(K):
         3
 
     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
-    itself [Rudolf et al.]_::
+    itself [Rudolf]_::
 
         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
@@ -343,7 +255,7 @@ def lyapunov_rank(K):
     TESTS:
 
     The Lyapunov rank should be additive on a product of proper cones
-    [Rudolf et al.]_::
+    [Rudolf]_::
 
         sage: set_random_seed()
         sage: K1 = random_cone(max_ambient_dim=8,
@@ -357,7 +269,7 @@ def lyapunov_rank(K):
         True
 
     The Lyapunov rank is invariant under a linear isomorphism
-    [Orlitzky/Gowda]_::
+    [Orlitzky]_::
 
         sage: K1 = random_cone(max_ambient_dim = 8)
         sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
@@ -366,7 +278,7 @@ def lyapunov_rank(K):
         True
 
     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
-    itself [Rudolf et al.]_::
+    itself [Rudolf]_::
 
         sage: set_random_seed()
         sage: K = random_cone(max_ambient_dim=8)
@@ -390,7 +302,7 @@ def lyapunov_rank(K):
         sage: b == n-1
         False
 
-    In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
+    In fact [Orlitzky]_, no closed convex polyhedral cone can have
     Lyapunov rank `n-1` in `n` dimensions::
 
         sage: set_random_seed()
@@ -401,7 +313,7 @@ def lyapunov_rank(K):
         False
 
     The calculation of the Lyapunov rank of an improper cone can be
-    reduced to that of a proper cone [Orlitzky/Gowda]_::
+    reduced to that of a proper cone [Orlitzky]_::
 
         sage: set_random_seed()
         sage: K = random_cone(max_ambient_dim=8)
@@ -414,15 +326,15 @@ def lyapunov_rank(K):
         sage: actual == expected
         True
 
-    The Lyapunov rank of any cone is just the dimension of ``K.LL()``::
+    The Lyapunov rank of a cone is the size of a :meth:`lyapunov_like_basis`::
 
         sage: set_random_seed()
         sage: K = random_cone(max_ambient_dim=8)
-        sage: lyapunov_rank(K) == len(K.LL())
+        sage: lyapunov_rank(K) == len(K.lyapunov_like_basis())
         True
 
     We can make an imperfect cone perfect by adding a slack variable
-    (a Theorem in [Orlitzky/Gowda]_)::
+    (a Theorem in [Orlitzky]_)::
 
         sage: set_random_seed()
         sage: K = random_cone(max_ambient_dim=8,
@@ -434,7 +346,7 @@ def lyapunov_rank(K):
         True
 
     """
-    beta = 0
+    beta = 0 # running tally of the Lyapunov rank
 
     m = K.dim()
     n = K.lattice_dim()
@@ -456,7 +368,7 @@ def lyapunov_rank(K):
         # Non-pointed reduction lemma.
         beta += l * m
 
-    beta += len(K.LL())
+    beta += len(K.lyapunov_like_basis())
     return beta
 
 
@@ -493,8 +405,8 @@ def is_lyapunov_like(L,K):
 
     REFERENCES:
 
-    .. [Orlitzky] M. Orlitzky. The Lyapunov rank of an
-       improper cone (preprint).
+    M. Orlitzky. The Lyapunov rank of an improper cone.
+    http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 
     EXAMPLES:
 
@@ -514,10 +426,11 @@ def is_lyapunov_like(L,K):
         sage: is_lyapunov_like(L,K)
         True
 
-    Everything in ``K.LL()`` 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 K.LL()])
+        sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
         True
 
     """
@@ -670,3 +583,104 @@ def positive_operators(K):
     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() ]