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 LL(K):
  84     r"""
  85     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
  86     on this cone.
  87
  88     OUTPUT:
  89
  90     A list of matrices forming a basis for the space of all
  91     Lyapunov-like transformations on the given cone.
  92
  93     EXAMPLES:
  94
  95     The trivial cone has no Lyapunov-like transformations::
  96
  97         sage: L = ToricLattice(0)
  98         sage: K = Cone([], lattice=L)
  99         sage: LL(K)
 100         []
 101
 102     The Lyapunov-like transformations on the nonnegative orthant are
 103     simply diagonal matrices::
 104
 105         sage: K = Cone([(1,)])
 106         sage: LL(K)
 107         [[1]]
 108
 109         sage: K = Cone([(1,0),(0,1)])
 110         sage: LL(K)
 111         [
 112         [1 0]  [0 0]
 113         [0 0], [0 1]
 114         ]
 115
 116         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 117         sage: LL(K)
 118         [
 119         [1 0 0]  [0 0 0]  [0 0 0]
 120         [0 0 0]  [0 1 0]  [0 0 0]
 121         [0 0 0], [0 0 0], [0 0 1]
 122         ]
 123
 124     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 125     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 126
 127         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 128         sage: LL(L31)
 129         [
 130         [1 0 0]
 131         [0 1 0]
 132         [0 0 1]
 133         ]
 134
 135         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 136         sage: LL(L3infty)
 137         [
 138         [1 0 0]
 139         [0 1 0]
 140         [0 0 1]
 141         ]
 142
 143     TESTS:
 144
 145     The inner product `\left< L\left(x\right), s \right>` is zero for
 146     every pair `\left( x,s \right)` in the discrete complementarity set
 147     of the cone::
 148
 149         sage: K = random_cone(max_dim=8, max_rays=10)
 150         sage: C_of_K = discrete_complementarity_set(K)
 151         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 152         sage: sum(map(abs, l))
 153         0
 154
 155     """
 156     V = K.lattice().vector_space()
 157
 158     C_of_K = discrete_complementarity_set(K)
 159
 160     tensor_products = [s.tensor_product(x) for (x,s) in C_of_K]
 161
 162     # Sage doesn't think matrices are vectors, so we have to convert
 163     # our matrices to vectors explicitly before we can figure out how
 164     # many are linearly-indepenedent.
 165     #
 166     # The space W has the same base ring as V, but dimension
 167     # dim(V)^2. So it has the same dimension as the space of linear
 168     # transformations on V. In other words, it's just the right size
 169     # to create an isomorphism between it and our matrices.
 170     W = VectorSpace(V.base_ring(), V.dimension()**2)
 171
 172     # Turn our matrices into long vectors...
 173     vectors = [ W(m.list()) for m in tensor_products ]
 174
 175     # Vector space representation of Lyapunov-like matrices
 176     # (i.e. vec(L) where L is Luapunov-like).
 177     LL_vector = W.span(vectors).complement()
 178
 179     # Now construct an ambient MatrixSpace in which to stick our
 180     # transformations.
 181     M = MatrixSpace(V.base_ring(), V.dimension())
 182
 183     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 184
 185     return matrix_basis
 186
 187
 188
 189 def lyapunov_rank(K):
 190     r"""
 191     Compute the Lyapunov (or bilinearity) rank of this cone.
 192
 193     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 194
 195     1. The dimension of the Lie algebra of the automorphism group of the
 196        cone.
 197
 198     2. The dimension of the linear space of all Lyapunov-like
 199        transformations on the cone.
 200
 201     INPUT:
 202
 203     A closed, convex polyhedral cone.
 204
 205     OUTPUT:
 206
 207     An integer representing the Lyapunov rank of the cone. If the
 208     dimension of the ambient vector space is `n`, then the Lyapunov rank
 209     will be between `1` and `n` inclusive; however a rank of `n-1` is
 210     not possible (see the first reference).
 211
 212     .. note::
 213
 214         In the references, the cones are always assumed to be proper. We
 215         do not impose this restriction.
 216
 217     .. seealso::
 218
 219         :meth:`is_proper`
 220
 221     ALGORITHM:
 222
 223     The codimension formula from the second reference is used. We find
 224     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 225     and `s` are rays of our cone. It is known that these vectors are
 226     sufficient to apply the codimension formula. Once we have all such
 227     pairs, we "brute force" the codimension formula by finding all
 228     linearly-independent `xs^{T}`.
 229
 230     REFERENCES:
 231
 232     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 233        cone and Lyapunov-like transformations, Mathematical Programming, 147
 234        (2014) 155-170.
 235
 236     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 237        optimality constraints for the cone of positive polynomials,
 238        Mathematical Programming, Series B, 129 (2011) 5-31.
 239
 240     EXAMPLES:
 241
 242     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 243     [Rudolf et al.]_::
 244
 245         sage: positives = Cone([(1,)])
 246         sage: lyapunov_rank(positives)
 247         1
 248         sage: quadrant = Cone([(1,0), (0,1)])
 249         sage: lyapunov_rank(quadrant)
 250         2
 251        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 252         sage: lyapunov_rank(octant)
 253         3
 254
 255     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 256     [Rudolf et al.]_::
 257
 258         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 259         sage: lyapunov_rank(L31)
 260         1
 261
 262     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 263
 264         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 265         sage: lyapunov_rank(L3infty)
 266         1
 267
 268     The Lyapunov rank should be additive on a product of cones
 269     [Rudolf et al.]_::
 270
 271         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 272         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 273         sage: K = L31.cartesian_product(octant)
 274         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 275         True
 276
 277     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 278     The cone ``K`` in the following example is isomorphic to the nonnegative
 279     octant in `\mathbb{R}^{3}`::
 280
 281         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 282         sage: lyapunov_rank(K)
 283         3
 284
 285     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 286     itself [Rudolf et al.]_::
 287
 288         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 289         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 290         True
 291
 292     TESTS:
 293
 294     The Lyapunov rank should be additive on a product of cones
 295     [Rudolf et al.]_::
 296
 297         sage: K1 = random_cone(max_dim=10, max_rays=10)
 298         sage: K2 = random_cone(max_dim=10, max_rays=10)
 299         sage: K = K1.cartesian_product(K2)
 300         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 301         True
 302
 303     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 304     itself [Rudolf et al.]_::
 305
 306         sage: K = random_cone(max_dim=10, max_rays=10)
 307         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 308         True
 309
 310     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 311     be any number between `1` and `n` inclusive, excluding `n-1`
 312     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 313     trivial cone in a trivial space as well. However, in zero dimensions,
 314     the Lyapunov rank of the trivial cone will be zero::
 315
 316         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 317         sage: b = lyapunov_rank(K)
 318         sage: n = K.lattice_dim()
 319         sage: (n == 0 or 1 <= b) and b <= n
 320         True
 321         sage: b == n-1
 322         False
 323
 324     """
 325     return len(LL(K))