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 discrete_complementarity_set(K):
  12     r"""
  13     Compute the discrete complementarity set of this cone.
  14
  15     The complementarity set of this cone is the set of all orthogonal
  16     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
  17     dual. The discrete complementarity set restricts `x` and `s` to be
  18     generators of their respective cones.
  19
  20     OUTPUT:
  21
  22     A list of pairs `(x,s)` such that,
  23
  24       * `x` is in this cone.
  25       * `x` is a generator of this cone.
  26       * `s` is in this cone's dual.
  27       * `s` is a generator of this cone's dual.
  28       * `x` and `s` are orthogonal.
  29
  30     EXAMPLES:
  31
  32     The discrete complementarity set of the nonnegative orthant consists
  33     of pairs of standard basis vectors::
  34
  35         sage: K = Cone([(1,0),(0,1)])
  36         sage: discrete_complementarity_set(K)
  37         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
  38
  39     If the cone consists of a single ray, the second components of the
  40     discrete complementarity set should generate the orthogonal
  41     complement of that ray::
  42
  43         sage: K = Cone([(1,0)])
  44         sage: discrete_complementarity_set(K)
  45         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
  46         sage: K = Cone([(1,0,0)])
  47         sage: discrete_complementarity_set(K)
  48         [((1, 0, 0), (0, 1, 0)),
  49           ((1, 0, 0), (0, -1, 0)),
  50           ((1, 0, 0), (0, 0, 1)),
  51           ((1, 0, 0), (0, 0, -1))]
  52
  53     When the cone is the entire space, its dual is the trivial cone, so
  54     the discrete complementarity set is empty::
  55
  56         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
  57         sage: discrete_complementarity_set(K)
  58         []
  59
  60     TESTS:
  61
  62     The complementarity set of the dual can be obtained by switching the
  63     components of the complementarity set of the original cone::
  64
  65         sage: K1 = random_cone(max_dim=10, max_rays=10)
  66         sage: K2 = K1.dual()
  67         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
  68         sage: actual = discrete_complementarity_set(K1)
  69         sage: actual == expected
  70         True
  71
  72     """
  73     V = K.lattice().vector_space()
  74
  75     # Convert the rays to vectors so that we can compute inner
  76     # products.
  77     xs = [V(x) for x in K.rays()]
  78     ss = [V(s) for s in K.dual().rays()]
  79
  80     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
  81
  82
  83 def lyapunov_rank(K):
  84     r"""
  85     Compute the Lyapunov (or bilinearity) rank of this cone.
  86
  87     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
  88
  89     1. The dimension of the Lie algebra of the automorphism group of the
  90        cone.
  91
  92     2. The dimension of the linear space of all Lyapunov-like
  93        transformations on the cone.
  94
  95     INPUT:
  96
  97     A closed, convex polyhedral cone.
  98
  99     OUTPUT:
 100
 101     An integer representing the Lyapunov rank of the cone. If the
 102     dimension of the ambient vector space is `n`, then the Lyapunov rank
 103     will be between `1` and `n` inclusive; however a rank of `n-1` is
 104     not possible (see the first reference).
 105
 106     .. note::
 107
 108         In the references, the cones are always assumed to be proper. We
 109         do not impose this restriction.
 110
 111     .. seealso::
 112
 113         :meth:`is_proper`
 114
 115     ALGORITHM:
 116
 117     The codimension formula from the second reference is used. We find
 118     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 119     and `s` are rays of our cone. It is known that these vectors are
 120     sufficient to apply the codimension formula. Once we have all such
 121     pairs, we "brute force" the codimension formula by finding all
 122     linearly-independent `xs^{T}`.
 123
 124     REFERENCES:
 125
 126     1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone
 127        and Lyapunov-like transformations, Mathematical Programming, 147
 128        (2014) 155-170.
 129
 130     2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 131        optimality constraints for the cone of positive polynomials,
 132        Mathematical Programming, Series B, 129 (2011) 5-31.
 133
 134     EXAMPLES:
 135
 136     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`::
 137
 138         sage: positives = Cone([(1,)])
 139         sage: lyapunov_rank(positives)
 140         1
 141         sage: quadrant = Cone([(1,0), (0,1)])
 142         sage: lyapunov_rank(quadrant)
 143         2
 144         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 145         sage: lyapunov_rank(octant)
 146         3
 147
 148     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one::
 149
 150         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 151         sage: lyapunov_rank(L31)
 152         1
 153
 154     Likewise for the `L^{3}_{\infty}` cone::
 155
 156         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 157         sage: lyapunov_rank(L3infty)
 158         1
 159
 160     The Lyapunov rank should be additive on a product of cones::
 161
 162         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 163         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 164         sage: K = L31.cartesian_product(octant)
 165         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 166         True
 167
 168     Two isomorphic cones should have the same Lyapunov rank. The cone
 169     ``K`` in the following example is isomorphic to the nonnegative
 170     octant in `\mathbb{R}^{3}`::
 171
 172         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 173         sage: lyapunov_rank(K)
 174         3
 175
 176     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 177     itself::
 178
 179         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 180         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 181         True
 182
 183     TESTS:
 184
 185     The Lyapunov rank should be additive on a product of cones::
 186
 187         sage: K1 = random_cone(max_dim=10, max_rays=10)
 188         sage: K2 = random_cone(max_dim=10, max_rays=10)
 189         sage: K = K1.cartesian_product(K2)
 190         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 191         True
 192
 193     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 194     itself::
 195
 196         sage: K = random_cone(max_dim=10, max_rays=10)
 197         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 198         True
 199
 200     """
 201     V = K.lattice().vector_space()
 202
 203     C_of_K = discrete_complementarity_set(K)
 204
 205     matrices = [x.tensor_product(s) for (x,s) in C_of_K]
 206
 207     # Sage doesn't think matrices are vectors, so we have to convert
 208     # our matrices to vectors explicitly before we can figure out how
 209     # many are linearly-indepenedent.
 210     #
 211     # The space W has the same base ring as V, but dimension
 212     # dim(V)^2. So it has the same dimension as the space of linear
 213     # transformations on V. In other words, it's just the right size
 214     # to create an isomorphism between it and our matrices.
 215     W = VectorSpace(V.base_ring(), V.dimension()**2)
 216
 217     def phi(m):
 218         r"""
 219         Convert a matrix to a vector isomorphically.
 220         """
 221         return W(m.list())
 222
 223     vectors = [phi(m) for m in matrices]
 224
 225     return (W.dimension() - W.span(vectors).rank())