From 012175f3a2591586099b4e955bb440958f2d63d7 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 11 May 2015 10:59:20 -0400 Subject: [PATCH] Add basic Lyapunov rank implementation for polyhedral cones. --- mjo/all.py | 1 + mjo/cone/cone.py | 141 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+) create mode 100644 mjo/cone/cone.py diff --git a/mjo/all.py b/mjo/all.py index c0676a2..cf758a8 100644 --- 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 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()) -- 2.43.2