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):
  12     r"""
  13     Project ``K`` into its own span.
  14
  15     EXAMPLES::
  16
  17         sage: K = Cone([(1,)])
  18         sage: project_span(K) == K
  19         True
  20
  21         sage: K2 = Cone([(1,0)])
  22         sage: project_span(K2).rays()
  23         N(1)
  24         in 1-d lattice N
  25         sage: K3 = Cone([(1,0,0)])
  26         sage: project_span(K3).rays()
  27         N(1)
  28         in 1-d lattice N
  29         sage: project_span(K2) == project_span(K3)
  30         True
  31
  32     TESTS:
  33
  34     The projected cone should always be solid::
  35
  36         sage: K = random_cone(max_dim = 10)
  37         sage: K_S = project_span(K)
  38         sage: K_S.is_solid()
  39         True
  40
  41     If we do this according to our paper, then the result is proper::
  42
  43         sage: K = random_cone(max_dim = 10)
  44         sage: K_S = project_span(K)
  45         sage: P = project_span(K_S.dual()).dual()
  46         sage: P.is_proper()
  47         True
  48
  49     """
  50     L = K.lattice()
  51     F = L.base_field()
  52     Q = L.quotient(K.sublattice_complement())
  53     vecs = [ vector(F, reversed(list(Q(r)))) for r in K.rays() ]
  54
  55     newL = None
  56     if len(vecs) == 0:
  57         newL = ToricLattice(0)
  58
  59     return Cone(vecs, lattice=newL)
  60
  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 (see the first reference).
 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 full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 311     [Orlitzky/Gowda]_::
 312
 313         sage: R5 = VectorSpace(QQ, 5)
 314         sage: gens = R5.basis() + [ -r for r in R5.basis() ]
 315         sage: K = Cone(gens)
 316         sage: lyapunov_rank(K)
 317         25
 318
 319     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 320     [Rudolf et al.]_::
 321
 322         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 323         sage: lyapunov_rank(L31)
 324         1
 325
 326     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 327
 328         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 329         sage: lyapunov_rank(L3infty)
 330         1
 331
 332     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 333     + 1` [Orlitzky/Gowda]_::
 334
 335         sage: K = Cone([(1,0,0,0,0)])
 336         sage: lyapunov_rank(K)
 337         21
 338         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 339         21
 340
 341     A subspace (of dimension `m`) in `n` dimensions should have a
 342     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 343
 344         sage: e1 = (1,0,0,0,0)
 345         sage: neg_e1 = (-1,0,0,0,0)
 346         sage: e2 = (0,1,0,0,0)
 347         sage: neg_e2 = (0,-1,0,0,0)
 348         sage: zero = (0,0,0,0,0)
 349         sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
 350         sage: lyapunov_rank(K)
 351         19
 352         sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
 353         19
 354
 355     The Lyapunov rank should be additive on a product of proper cones
 356     [Rudolf et al.]_::
 357
 358         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 359         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 360         sage: K = L31.cartesian_product(octant)
 361         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 362         True
 363
 364     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 365     The cone ``K`` in the following example is isomorphic to the nonnegative
 366     octant in `\mathbb{R}^{3}`::
 367
 368         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 369         sage: lyapunov_rank(K)
 370         3
 371
 372     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 373     itself [Rudolf et al.]_::
 374
 375         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 376         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 377         True
 378
 379     TESTS:
 380
 381     The Lyapunov rank should be additive on a product of proper cones
 382     [Rudolf et al.]_::
 383
 384         sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
 385         sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
 386         sage: K = K1.cartesian_product(K2)
 387         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 388         True
 389
 390     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 391     itself [Rudolf et al.]_::
 392
 393         sage: K = random_cone(max_dim=10, max_rays=10)
 394         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 395         True
 396
 397     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 398     be any number between `1` and `n` inclusive, excluding `n-1`
 399     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 400     trivial cone in a trivial space as well. However, in zero dimensions,
 401     the Lyapunov rank of the trivial cone will be zero::
 402
 403         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 404         sage: b = lyapunov_rank(K)
 405         sage: n = K.lattice_dim()
 406         sage: (n == 0 or 1 <= b) and b <= n
 407         True
 408         sage: b == n-1
 409         False
 410
 411     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 412     Lyapunov rank `n-1` in `n` dimensions::
 413
 414         sage: K = random_cone(max_dim=10)
 415         sage: b = lyapunov_rank(K)
 416         sage: n = K.lattice_dim()
 417         sage: b == n-1
 418         False
 419
 420     The calculation of the Lyapunov rank of an improper cone can be
 421     reduced to that of a proper cone [Orlitzky/Gowda]_::
 422
 423         sage: K = random_cone(max_dim=10)
 424         sage: actual = lyapunov_rank(K)
 425         sage: K_S = project_span(K)
 426         sage: P = project_span(K_S.dual()).dual()
 427         sage: l = K.linear_subspace().dimension()
 428         sage: codim = K.lattice_dim() - K.dim()
 429         sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
 430         sage: actual == expected
 431         True
 432
 433     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 434
 435         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 436         sage: lyapunov_rank(K) == len(LL(K))
 437         True
 438
 439     """
 440     beta = 0
 441
 442     m = K.dim()
 443     n = K.lattice_dim()
 444     l = K.linear_subspace().dimension()
 445
 446     if m < n:
 447         # K is not solid, project onto its span.
 448         K = project_span(K)
 449
 450         # Lemma 2
 451         beta += m*(n - m) + (n - m)**2
 452
 453     if l > 0:
 454         # K is not pointed, project its dual onto its span.
 455         K = project_span(K.dual()).dual()
 456
 457         # Lemma 3
 458         beta += m * l
 459
 460     beta += len(LL(K))
 461     return beta