mjo/cone/cone.py

   1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
   2 # have to explicitly mangle our sitedir here so that "mjo.cone"
   3 # resolves.
   4 from os.path import abspath
   5 from site import addsitedir
   6 addsitedir(abspath('../../'))
   7
   8 from sage.all import *
   9
  10
  11 def lyapunov_rank(K):
  12     r"""
  13     Compute the Lyapunov (or bilinearity) rank of this cone.
  14
  15     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
  16
  17     1. The dimension of the Lie algebra of the automorphism group of the
  18        cone.
  19
  20     2. The dimension of the linear space of all Lyapunov-like
  21        transformations on the cone.
  22
  23     INPUT:
  24
  25     A closed, convex polyhedral cone.
  26
  27     OUTPUT:
  28
  29     An integer representing the Lyapunov rank of the cone. If the
  30     dimension of the ambient vector space is `n`, then the Lyapunov rank
  31     will be between `1` and `n` inclusive; however a rank of `n-1` is
  32     not possible (see the first reference).
  33
  34     .. note::
  35
  36         In the references, the cones are always assumed to be proper. We
  37         do not impose this restriction.
  38
  39     .. seealso::
  40
  41         :meth:`is_proper`
  42
  43     ALGORITHM:
  44
  45     The codimension formula from the second reference is used. We find
  46     all pairs `(x,s)` in the complementarity set of `K` such that `x`
  47     and `s` are rays of our cone. It is known that these vectors are
  48     sufficient to apply the codimension formula. Once we have all such
  49     pairs, we "brute force" the codimension formula by finding all
  50     linearly-independent `xs^{T}`.
  51
  52     REFERENCES:
  53
  54     1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
  55        and Lyapunov-like transformations, Mathematical Programming, 147
  56        (2014) 155-170.
  57
  58     2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
  59        optimality constraints for the cone of positive polynomials,
  60        Mathematical Programming, Series B, 129 (2011) 5-31.
  61
  62     EXAMPLES:
  63
  64     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
  65
  66         sage: positives = Cone([(1,)])
  67         sage: lyapunov_rank(positives)
  68         1
  69         sage: quadrant = Cone([(1,0), (0,1)])
  70         sage: lyapunov_rank(quadrant)
  71         2
  72         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
  73         sage: lyapunov_rank(octant)
  74         3
  75
  76     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
  77
  78         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
  79         sage: lyapunov_rank(L31)
  80         1
  81
  82     Likewise for the `L^{3}_{\infty}` cone::
  83
  84         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
  85         sage: lyapunov_rank(L3infty)
  86         1
  87
  88     The Lyapunov rank should be additive on a product of cones::
  89
  90         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
  91         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
  92         sage: K = L31.cartesian_product(octant)
  93         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
  94         True
  95
  96     Two isomorphic cones should have the same Lyapunov rank. The cone
  97     ``K`` in the following example is isomorphic to the nonnegative
  98     octant in `\mathbb{R}^{3}`::
  99
 100         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 101         sage: lyapunov_rank(K)
 102         3
 103
 104     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 105     itself::
 106
 107         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 108         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 109         True
 110
 111     """
 112     V = K.lattice().vector_space()
 113
 114     xs = [V(x) for x in K.rays()]
 115     ss = [V(s) for s in K.dual().rays()]
 116
 117     # WARNING: This isn't really C(K), it only contains the pairs
 118     # (x,s) in C(K) where x,s are extreme in their respective cones.
 119     C_of_K = [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 120
 121     matrices = [x.column() * s.row() for (x,s) in C_of_K]
 122
 123     # Sage doesn't think matrices are vectors, so we have to convert
 124     # our matrices to vectors explicitly before we can figure out how
 125     # many are linearly-indepenedent.
 126     #
 127     # The space W has the same base ring as V, but dimension
 128     # dim(V)^2. So it has the same dimension as the space of linear
 129     # transformations on V. In other words, it's just the right size
 130     # to create an isomorphism between it and our matrices.
 131     W = VectorSpace(V.base_ring(), V.dimension()**2)
 132
 133     def phi(m):
 134         r"""
 135         Convert a matrix to a vector isomorphically.
 136         """
 137         return W(m.list())
 138
 139     vectors = [phi(m) for m in matrices]
 140
 141     return (W.dimension() - W.span(vectors).rank())