Add basic Lyapunov rank implementation for polyhedral cones.

diff --git a/mjo/all.py b/mjo/all.py
index c0676a2430b4249d16b52624f5f70cb42c9f58d0..cf758a85865e6f9cf86dca89e5ce8bb408f5d635 100644 (file)
--- a/mjo/all.py
+++ b/mjo/all.py
@@ -3,6 +3,7 @@ Import all of the other code, so that the user doesn't have to do it
 in his script. Instead, he can just `from mjo.all import *`.
 """
 
+from cone.cone import *
 from cone.symmetric_psd import *
 from cone.doubly_nonnegative import *
 from cone.completely_positive import *

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
new file mode 100644 (file)
index 0000000..2304338
--- /dev/null
+++ b/mjo/cone/cone.py
@@ -0,0 +1,141 @@
+# 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 *
+
+
+def lyapunov_rank(K):
+    r"""
+    Compute the Lyapunov (or bilinearity) rank of this cone.
+
+    The Lyapunov rank of a cone can be thought of in (mainly) two ways:
+
+    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.
+
+    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 the first reference).
+
+    .. note::
+
+        In the references, the cones are always assumed to be proper. We
+        do not impose this restriction.
+
+    .. seealso::
+
+        :meth:`is_proper`
+
+    ALGORITHM:
+
+    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}`.
+
+    REFERENCES:
+
+    1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
+       and Lyapunov-like transformations, Mathematical Programming, 147
+       (2014) 155-170.
+
+    2. 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`::
+
+        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 `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
+
+        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::
+
+        sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
+        sage: lyapunov_rank(L3infty)
+        1
+
+    The Lyapunov rank should be additive on a product of cones::
+
+        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. 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::
+
+        sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
+        sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
+        True
+
+    """
+    V = K.lattice().vector_space()
+
+    xs = [V(x) for x in K.rays()]
+    ss = [V(s) for s in K.dual().rays()]
+
+    # WARNING: This isn't really C(K), it only contains the pairs
+    # (x,s) in C(K) where x,s are extreme in their respective cones.
+    C_of_K = [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
+
+    matrices = [x.column() * s.row() 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)
+
+    def phi(m):
+        r"""
+        Convert a matrix to a vector isomorphically.
+        """
+        return W(m.list())
+
+    vectors = [phi(m) for m in matrices]
+
+    return (W.dimension() - W.span(vectors).rank())