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 project_span(K, K2 = None):
  12     r"""
  13     Return a "copy" of ``K`` embeded in a lower-dimensional space.
  14
  15     By default, we will project ``K`` into the subspace spanned by its
  16     rays. However, if ``K2`` is not ``None``, we will project into the
  17     space spanned by the rays of ``K2`` instead.
  18
  19     EXAMPLES::
  20
  21         sage: K = Cone([(1,0,0), (0,1,0)])
  22         sage: project_span(K)
  23         2-d cone in 2-d lattice N
  24         sage: project_span(K).rays()
  25         N(1, 0),
  26         N(0, 1)
  27         in 2-d lattice N
  28
  29         sage: K = Cone([(1,0,0), (0,1,0)])
  30         sage: K2 = Cone([(0,1)])
  31         sage: project_span(K, K2).rays()
  32         N(1)
  33         in 1-d lattice N
  34
  35     """
  36     # Allow us to use a second cone to generate the subspace into
  37     # which we're "projecting."
  38     if K2 is None:
  39         K2 = K
  40
  41     # Use these to generate the new cone.
  42     cs1 = K.rays().matrix().columns()
  43
  44     # And use these to figure out which indices to drop.
  45     cs2 = K2.rays().matrix().columns()
  46
  47     perp_idxs = []
  48
  49     for idx in range(0, len(cs2)):
  50         if cs2[idx].is_zero():
  51             perp_idxs.append(idx)
  52
  53     solid_cols = [ cs1[idx] for idx in range(0,len(cs1))
  54                             if not idx in perp_idxs
  55                             and not idx >= len(cs2) ]
  56
  57     m = matrix(solid_cols)
  58     L = ToricLattice(len(m.rows()))
  59     J = Cone(m.transpose(), lattice=L)
  60     return J
  61
  62
  63 def discrete_complementarity_set(K):
  64     r"""
  65     Compute the discrete complementarity set of this cone.
  66
  67     The complementarity set of this cone is the set of all orthogonal
  68     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
  69     dual. The discrete complementarity set restricts `x` and `s` to be
  70     generators of their respective cones.
  71
  72     OUTPUT:
  73
  74     A list of pairs `(x,s)` such that,
  75
  76       * `x` is in this cone.
  77       * `x` is a generator of this cone.
  78       * `s` is in this cone's dual.
  79       * `s` is a generator of this cone's dual.
  80       * `x` and `s` are orthogonal.
  81
  82     EXAMPLES:
  83
  84     The discrete complementarity set of the nonnegative orthant consists
  85     of pairs of standard basis vectors::
  86
  87         sage: K = Cone([(1,0),(0,1)])
  88         sage: discrete_complementarity_set(K)
  89         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
  90
  91     If the cone consists of a single ray, the second components of the
  92     discrete complementarity set should generate the orthogonal
  93     complement of that ray::
  94
  95         sage: K = Cone([(1,0)])
  96         sage: discrete_complementarity_set(K)
  97         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
  98         sage: K = Cone([(1,0,0)])
  99         sage: discrete_complementarity_set(K)
 100         [((1, 0, 0), (0, 1, 0)),
 101           ((1, 0, 0), (0, -1, 0)),
 102           ((1, 0, 0), (0, 0, 1)),
 103           ((1, 0, 0), (0, 0, -1))]
 104
 105     When the cone is the entire space, its dual is the trivial cone, so
 106     the discrete complementarity set is empty::
 107
 108         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 109         sage: discrete_complementarity_set(K)
 110         []
 111
 112     TESTS:
 113
 114     The complementarity set of the dual can be obtained by switching the
 115     components of the complementarity set of the original cone::
 116
 117         sage: K1 = random_cone(max_dim=10, max_rays=10)
 118         sage: K2 = K1.dual()
 119         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 120         sage: actual = discrete_complementarity_set(K1)
 121         sage: actual == expected
 122         True
 123
 124     """
 125     V = K.lattice().vector_space()
 126
 127     # Convert the rays to vectors so that we can compute inner
 128     # products.
 129     xs = [V(x) for x in K.rays()]
 130     ss = [V(s) for s in K.dual().rays()]
 131
 132     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 133
 134
 135 def LL(K):
 136     r"""
 137     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 138     on this cone.
 139
 140     OUTPUT:
 141
 142     A list of matrices forming a basis for the space of all
 143     Lyapunov-like transformations on the given cone.
 144
 145     EXAMPLES:
 146
 147     The trivial cone has no Lyapunov-like transformations::
 148
 149         sage: L = ToricLattice(0)
 150         sage: K = Cone([], lattice=L)
 151         sage: LL(K)
 152         []
 153
 154     The Lyapunov-like transformations on the nonnegative orthant are
 155     simply diagonal matrices::
 156
 157         sage: K = Cone([(1,)])
 158         sage: LL(K)
 159         [[1]]
 160
 161         sage: K = Cone([(1,0),(0,1)])
 162         sage: LL(K)
 163         [
 164         [1 0]  [0 0]
 165         [0 0], [0 1]
 166         ]
 167
 168         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 169         sage: LL(K)
 170         [
 171         [1 0 0]  [0 0 0]  [0 0 0]
 172         [0 0 0]  [0 1 0]  [0 0 0]
 173         [0 0 0], [0 0 0], [0 0 1]
 174         ]
 175
 176     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 177     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 178
 179         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 180         sage: LL(L31)
 181         [
 182         [1 0 0]
 183         [0 1 0]
 184         [0 0 1]
 185         ]
 186
 187         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 188         sage: LL(L3infty)
 189         [
 190         [1 0 0]
 191         [0 1 0]
 192         [0 0 1]
 193         ]
 194
 195     TESTS:
 196
 197     The inner product `\left< L\left(x\right), s \right>` is zero for
 198     every pair `\left( x,s \right)` in the discrete complementarity set
 199     of the cone::
 200
 201         sage: K = random_cone(max_dim=8, max_rays=10)
 202         sage: C_of_K = discrete_complementarity_set(K)
 203         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 204         sage: sum(map(abs, l))
 205         0
 206
 207     """
 208     V = K.lattice().vector_space()
 209
 210     C_of_K = discrete_complementarity_set(K)
 211
 212     tensor_products = [s.tensor_product(x) for (x,s) in C_of_K]
 213
 214     # Sage doesn't think matrices are vectors, so we have to convert
 215     # our matrices to vectors explicitly before we can figure out how
 216     # many are linearly-indepenedent.
 217     #
 218     # The space W has the same base ring as V, but dimension
 219     # dim(V)^2. So it has the same dimension as the space of linear
 220     # transformations on V. In other words, it's just the right size
 221     # to create an isomorphism between it and our matrices.
 222     W = VectorSpace(V.base_ring(), V.dimension()**2)
 223
 224     # Turn our matrices into long vectors...
 225     vectors = [ W(m.list()) for m in tensor_products ]
 226
 227     # Vector space representation of Lyapunov-like matrices
 228     # (i.e. vec(L) where L is Luapunov-like).
 229     LL_vector = W.span(vectors).complement()
 230
 231     # Now construct an ambient MatrixSpace in which to stick our
 232     # transformations.
 233     M = MatrixSpace(V.base_ring(), V.dimension())
 234
 235     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 236
 237     return matrix_basis
 238
 239
 240
 241 def lyapunov_rank(K):
 242     r"""
 243     Compute the Lyapunov (or bilinearity) rank of this cone.
 244
 245     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 246
 247     1. The dimension of the Lie algebra of the automorphism group of the
 248        cone.
 249
 250     2. The dimension of the linear space of all Lyapunov-like
 251        transformations on the cone.
 252
 253     INPUT:
 254
 255     A closed, convex polyhedral cone.
 256
 257     OUTPUT:
 258
 259     An integer representing the Lyapunov rank of the cone. If the
 260     dimension of the ambient vector space is `n`, then the Lyapunov rank
 261     will be between `1` and `n` inclusive; however a rank of `n-1` is
 262     not possible for any cone.
 263
 264     .. note::
 265
 266         In the references, the cones are always assumed to be proper. We
 267         do not impose this restriction.
 268
 269     .. seealso::
 270
 271         :meth:`is_proper`
 272
 273     ALGORITHM:
 274
 275     The codimension formula from the second reference is used. We find
 276     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 277     and `s` are rays of our cone. It is known that these vectors are
 278     sufficient to apply the codimension formula. Once we have all such
 279     pairs, we "brute force" the codimension formula by finding all
 280     linearly-independent `xs^{T}`.
 281
 282     REFERENCES:
 283
 284     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 285        cone and Lyapunov-like transformations, Mathematical Programming, 147
 286        (2014) 155-170.
 287
 288     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 289        Improper Cone. Work in-progress.
 290
 291     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 292        optimality constraints for the cone of positive polynomials,
 293        Mathematical Programming, Series B, 129 (2011) 5-31.
 294
 295     EXAMPLES:
 296
 297     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 298     [Rudolf et al.]_::
 299
 300         sage: positives = Cone([(1,)])
 301         sage: lyapunov_rank(positives)
 302         1
 303         sage: quadrant = Cone([(1,0), (0,1)])
 304         sage: lyapunov_rank(quadrant)
 305         2
 306        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 307         sage: lyapunov_rank(octant)
 308         3
 309
 310     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 311     [Rudolf et al.]_::
 312
 313         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 314         sage: lyapunov_rank(L31)
 315         1
 316
 317     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 318
 319         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 320         sage: lyapunov_rank(L3infty)
 321         1
 322
 323     The Lyapunov rank should be additive on a product of cones
 324     [Rudolf et al.]_::
 325
 326         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 327         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 328         sage: K = L31.cartesian_product(octant)
 329         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 330         True
 331
 332     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 333     The cone ``K`` in the following example is isomorphic to the nonnegative
 334     octant in `\mathbb{R}^{3}`::
 335
 336         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 337         sage: lyapunov_rank(K)
 338         3
 339
 340     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 341     itself [Rudolf et al.]_::
 342
 343         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 344         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 345         True
 346
 347     TESTS:
 348
 349     The Lyapunov rank should be additive on a product of cones
 350     [Rudolf et al.]_::
 351
 352         sage: K1 = random_cone(max_dim=10, max_rays=10)
 353         sage: K2 = random_cone(max_dim=10, max_rays=10)
 354         sage: K = K1.cartesian_product(K2)
 355         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 356         True
 357
 358     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 359     itself [Rudolf et al.]_::
 360
 361         sage: K = random_cone(max_dim=10, max_rays=10)
 362         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 363         True
 364
 365     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 366     be any number between `1` and `n` inclusive, excluding `n-1`
 367     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 368     trivial cone in a trivial space as well. However, in zero dimensions,
 369     the Lyapunov rank of the trivial cone will be zero::
 370
 371         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 372         sage: b = lyapunov_rank(K)
 373         sage: n = K.lattice_dim()
 374         sage: (n == 0 or 1 <= b) and b <= n
 375         True
 376         sage: b == n-1
 377         False
 378
 379     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 380     Lyapunov rank `n-1` in `n` dimensions::
 381
 382         sage: K = random_cone(max_dim=10, max_rays=16)
 383         sage: b = lyapunov_rank(K)
 384         sage: n = K.lattice_dim()
 385         sage: b == n-1
 386         False
 387
 388     The calculation of the Lyapunov rank of an improper cone can be
 389     reduced to that of a proper cone [Orlitzky/Gowda]_::
 390
 391         sage: K = random_cone(max_dim=15, max_rays=25)
 392         sage: actual = lyapunov_rank(K)
 393         sage: K_S = project_span(K)
 394         sage: J_T1 = project_span(K_S.dual()).dual()
 395         sage: J_T2 = project_span(K, K_S.dual())
 396         sage: J_T2 = Cone(J_T2.rays(), lattice=J_T1.lattice())
 397         sage: J_T1 == J_T2
 398         True
 399         sage: J_T = J_T1
 400         sage: l = K.linear_subspace().dimension()
 401         sage: codim = K.lattice_dim() - K.dim()
 402         sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
 403         sage: actual == expected
 404         True
 405
 406     """
 407     return len(LL(K))