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