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 lineality(K):
  64     r"""
  65     Compute the lineality of this cone.
  66
  67     The lineality of a cone is the dimension of the largest linear
  68     subspace contained in that cone.
  69
  70     OUTPUT:
  71
  72     A nonnegative integer; the dimension of the largest subspace
  73     contained within this cone.
  74
  75     REFERENCES:
  76
  77     .. [Rockafellar] R.T. Rockafellar. Convex Analysis. Princeton
  78        University Press, Princeton, 1970.
  79
  80     EXAMPLES:
  81
  82     The lineality of the nonnegative orthant is zero, since it clearly
  83     contains no lines::
  84
  85         sage: K = Cone([(1,0,0), (0,1,0), (0,0,1)])
  86         sage: lineality(K)
  87         0
  88
  89     However, if we add another ray so that the entire `x`-axis belongs
  90     to the cone, then the resulting cone will have lineality one::
  91
  92         sage: K = Cone([(1,0,0), (-1,0,0), (0,1,0), (0,0,1)])
  93         sage: lineality(K)
  94         1
  95
  96     If our cone is all of `\mathbb{R}^{2}`, then its lineality is equal
  97     to the dimension of the ambient space (i.e. two)::
  98
  99         sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
 100         sage: lineality(K)
 101         2
 102
 103     Per the definition, the lineality of the trivial cone in a trivial
 104     space is zero::
 105
 106         sage: K = Cone([], lattice=ToricLattice(0))
 107         sage: lineality(K)
 108         0
 109
 110     TESTS:
 111
 112     The lineality of a cone should be an integer between zero and the
 113     dimension of the ambient space, inclusive::
 114
 115         sage: K = random_cone(max_dim = 10)
 116         sage: l = lineality(K)
 117         sage: l in ZZ
 118         True
 119         sage: (0 <= l) and (l <= K.lattice_dim())
 120         True
 121
 122     A strictly cone should have lineality zero::
 123
 124         sage: K = random_cone(max_dim = 10, strictly_convex = True)
 125         sage: lineality(K)
 126         0
 127
 128     """
 129     return K.linear_subspace().dimension()
 130
 131
 132 def discrete_complementarity_set(K):
 133     r"""
 134     Compute the discrete complementarity set of this cone.
 135
 136     The complementarity set of this cone is the set of all orthogonal
 137     pairs `(x,s)` such that `x` is in this cone, and `s` is in its
 138     dual. The discrete complementarity set restricts `x` and `s` to be
 139     generators of their respective cones.
 140
 141     OUTPUT:
 142
 143     A list of pairs `(x,s)` such that,
 144
 145       * `x` is in this cone.
 146       * `x` is a generator of this cone.
 147       * `s` is in this cone's dual.
 148       * `s` is a generator of this cone's dual.
 149       * `x` and `s` are orthogonal.
 150
 151     EXAMPLES:
 152
 153     The discrete complementarity set of the nonnegative orthant consists
 154     of pairs of standard basis vectors::
 155
 156         sage: K = Cone([(1,0),(0,1)])
 157         sage: discrete_complementarity_set(K)
 158         [((1, 0), (0, 1)), ((0, 1), (1, 0))]
 159
 160     If the cone consists of a single ray, the second components of the
 161     discrete complementarity set should generate the orthogonal
 162     complement of that ray::
 163
 164         sage: K = Cone([(1,0)])
 165         sage: discrete_complementarity_set(K)
 166         [((1, 0), (0, 1)), ((1, 0), (0, -1))]
 167         sage: K = Cone([(1,0,0)])
 168         sage: discrete_complementarity_set(K)
 169         [((1, 0, 0), (0, 1, 0)),
 170           ((1, 0, 0), (0, -1, 0)),
 171           ((1, 0, 0), (0, 0, 1)),
 172           ((1, 0, 0), (0, 0, -1))]
 173
 174     When the cone is the entire space, its dual is the trivial cone, so
 175     the discrete complementarity set is empty::
 176
 177         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 178         sage: discrete_complementarity_set(K)
 179         []
 180
 181     TESTS:
 182
 183     The complementarity set of the dual can be obtained by switching the
 184     components of the complementarity set of the original cone::
 185
 186         sage: K1 = random_cone(max_dim=10, max_rays=10)
 187         sage: K2 = K1.dual()
 188         sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
 189         sage: actual = discrete_complementarity_set(K1)
 190         sage: actual == expected
 191         True
 192
 193     """
 194     V = K.lattice().vector_space()
 195
 196     # Convert the rays to vectors so that we can compute inner
 197     # products.
 198     xs = [V(x) for x in K.rays()]
 199     ss = [V(s) for s in K.dual().rays()]
 200
 201     return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0]
 202
 203
 204 def LL(K):
 205     r"""
 206     Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
 207     on this cone.
 208
 209     OUTPUT:
 210
 211     A list of matrices forming a basis for the space of all
 212     Lyapunov-like transformations on the given cone.
 213
 214     EXAMPLES:
 215
 216     The trivial cone has no Lyapunov-like transformations::
 217
 218         sage: L = ToricLattice(0)
 219         sage: K = Cone([], lattice=L)
 220         sage: LL(K)
 221         []
 222
 223     The Lyapunov-like transformations on the nonnegative orthant are
 224     simply diagonal matrices::
 225
 226         sage: K = Cone([(1,)])
 227         sage: LL(K)
 228         [[1]]
 229
 230         sage: K = Cone([(1,0),(0,1)])
 231         sage: LL(K)
 232         [
 233         [1 0]  [0 0]
 234         [0 0], [0 1]
 235         ]
 236
 237         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 238         sage: LL(K)
 239         [
 240         [1 0 0]  [0 0 0]  [0 0 0]
 241         [0 0 0]  [0 1 0]  [0 0 0]
 242         [0 0 0], [0 0 0], [0 0 1]
 243         ]
 244
 245     Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
 246     `L^{3}_{\infty}` cones [Rudolf et al.]_::
 247
 248         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 249         sage: LL(L31)
 250         [
 251         [1 0 0]
 252         [0 1 0]
 253         [0 0 1]
 254         ]
 255
 256         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 257         sage: LL(L3infty)
 258         [
 259         [1 0 0]
 260         [0 1 0]
 261         [0 0 1]
 262         ]
 263
 264     TESTS:
 265
 266     The inner product `\left< L\left(x\right), s \right>` is zero for
 267     every pair `\left( x,s \right)` in the discrete complementarity set
 268     of the cone::
 269
 270         sage: K = random_cone(max_dim=8, max_rays=10)
 271         sage: C_of_K = discrete_complementarity_set(K)
 272         sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
 273         sage: sum(map(abs, l))
 274         0
 275
 276     """
 277     V = K.lattice().vector_space()
 278
 279     C_of_K = discrete_complementarity_set(K)
 280
 281     tensor_products = [s.tensor_product(x) for (x,s) in C_of_K]
 282
 283     # Sage doesn't think matrices are vectors, so we have to convert
 284     # our matrices to vectors explicitly before we can figure out how
 285     # many are linearly-indepenedent.
 286     #
 287     # The space W has the same base ring as V, but dimension
 288     # dim(V)^2. So it has the same dimension as the space of linear
 289     # transformations on V. In other words, it's just the right size
 290     # to create an isomorphism between it and our matrices.
 291     W = VectorSpace(V.base_ring(), V.dimension()**2)
 292
 293     # Turn our matrices into long vectors...
 294     vectors = [ W(m.list()) for m in tensor_products ]
 295
 296     # Vector space representation of Lyapunov-like matrices
 297     # (i.e. vec(L) where L is Luapunov-like).
 298     LL_vector = W.span(vectors).complement()
 299
 300     # Now construct an ambient MatrixSpace in which to stick our
 301     # transformations.
 302     M = MatrixSpace(V.base_ring(), V.dimension())
 303
 304     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 305
 306     return matrix_basis
 307
 308
 309
 310 def lyapunov_rank(K):
 311     r"""
 312     Compute the Lyapunov (or bilinearity) rank of this cone.
 313
 314     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 315
 316     1. The dimension of the Lie algebra of the automorphism group of the
 317        cone.
 318
 319     2. The dimension of the linear space of all Lyapunov-like
 320        transformations on the cone.
 321
 322     INPUT:
 323
 324     A closed, convex polyhedral cone.
 325
 326     OUTPUT:
 327
 328     An integer representing the Lyapunov rank of the cone. If the
 329     dimension of the ambient vector space is `n`, then the Lyapunov rank
 330     will be between `1` and `n` inclusive; however a rank of `n-1` is
 331     not possible (see the first reference).
 332
 333     .. note::
 334
 335         In the references, the cones are always assumed to be proper. We
 336         do not impose this restriction.
 337
 338     .. seealso::
 339
 340         :meth:`is_proper`
 341
 342     ALGORITHM:
 343
 344     The codimension formula from the second reference is used. We find
 345     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 346     and `s` are rays of our cone. It is known that these vectors are
 347     sufficient to apply the codimension formula. Once we have all such
 348     pairs, we "brute force" the codimension formula by finding all
 349     linearly-independent `xs^{T}`.
 350
 351     REFERENCES:
 352
 353     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 354        cone and Lyapunov-like transformations, Mathematical Programming, 147
 355        (2014) 155-170.
 356
 357     .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
 358        Improper Cone. Work in-progress.
 359
 360     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 361        optimality constraints for the cone of positive polynomials,
 362        Mathematical Programming, Series B, 129 (2011) 5-31.
 363
 364     EXAMPLES:
 365
 366     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 367     [Rudolf et al.]_::
 368
 369         sage: positives = Cone([(1,)])
 370         sage: lyapunov_rank(positives)
 371         1
 372         sage: quadrant = Cone([(1,0), (0,1)])
 373         sage: lyapunov_rank(quadrant)
 374         2
 375        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 376         sage: lyapunov_rank(octant)
 377         3
 378
 379     The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
 380     [Orlitzky/Gowda]_::
 381
 382         sage: R5 = VectorSpace(QQ, 5)
 383         sage: gens = R5.basis() + [ -r for r in R5.basis() ]
 384         sage: K = Cone(gens)
 385         sage: lyapunov_rank(K)
 386         25
 387
 388     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 389     [Rudolf et al.]_::
 390
 391         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 392         sage: lyapunov_rank(L31)
 393         1
 394
 395     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 396
 397         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 398         sage: lyapunov_rank(L3infty)
 399         1
 400
 401     A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
 402     + 1` [Orlitzky/Gowda]_::
 403
 404         sage: K = Cone([(1,0,0,0,0)])
 405         sage: lyapunov_rank(K)
 406         21
 407         sage: K.lattice_dim()**2 - K.lattice_dim() + 1
 408         21
 409
 410     A subspace (of dimension `m`) in `n` dimensions should have a
 411     Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
 412
 413         sage: e1 = (1,0,0,0,0)
 414         sage: neg_e1 = (-1,0,0,0,0)
 415         sage: e2 = (0,1,0,0,0)
 416         sage: neg_e2 = (0,-1,0,0,0)
 417         sage: zero = (0,0,0,0,0)
 418         sage: K = Cone([e1, neg_e1, e2, neg_e2, zero, zero, zero])
 419         sage: lyapunov_rank(K)
 420         19
 421         sage: K.lattice_dim()**2 - K.dim()*(K.lattice_dim() - K.dim())
 422         19
 423
 424     The Lyapunov rank should be additive on a product of proper cones
 425     [Rudolf et al.]_::
 426
 427         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 428         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 429         sage: K = L31.cartesian_product(octant)
 430         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 431         True
 432
 433     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 434     The cone ``K`` in the following example is isomorphic to the nonnegative
 435     octant in `\mathbb{R}^{3}`::
 436
 437         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 438         sage: lyapunov_rank(K)
 439         3
 440
 441     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 442     itself [Rudolf et al.]_::
 443
 444         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 445         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 446         True
 447
 448     TESTS:
 449
 450     The Lyapunov rank should be additive on a product of proper cones
 451     [Rudolf et al.]_::
 452
 453         sage: K1 = random_cone(max_dim=10, strictly_convex=True, solid=True)
 454         sage: K2 = random_cone(max_dim=10, strictly_convex=True, solid=True)
 455         sage: K = K1.cartesian_product(K2)
 456         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 457         True
 458
 459     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 460     itself [Rudolf et al.]_::
 461
 462         sage: K = random_cone(max_dim=10, max_rays=10)
 463         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 464         True
 465
 466     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 467     be any number between `1` and `n` inclusive, excluding `n-1`
 468     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 469     trivial cone in a trivial space as well. However, in zero dimensions,
 470     the Lyapunov rank of the trivial cone will be zero::
 471
 472         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 473         sage: b = lyapunov_rank(K)
 474         sage: n = K.lattice_dim()
 475         sage: (n == 0 or 1 <= b) and b <= n
 476         True
 477         sage: b == n-1
 478         False
 479
 480     In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
 481     Lyapunov rank `n-1` in `n` dimensions::
 482
 483         sage: K = random_cone(max_dim=10)
 484         sage: b = lyapunov_rank(K)
 485         sage: n = K.lattice_dim()
 486         sage: b == n-1
 487         False
 488
 489     The calculation of the Lyapunov rank of an improper cone can be
 490     reduced to that of a proper cone [Orlitzky/Gowda]_::
 491
 492         sage: K = random_cone(max_dim=10)
 493         sage: actual = lyapunov_rank(K)
 494         sage: K_S = project_span(K)
 495         sage: P = project_span(K_S.dual()).dual()
 496         sage: l = lineality(K)
 497         sage: codim = K.lattice_dim() - K.dim()
 498         sage: expected = lyapunov_rank(P) + K.dim()*(l + codim) + codim**2
 499         sage: actual == expected
 500         True
 501
 502     The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
 503
 504         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 505         sage: lyapunov_rank(K) == len(LL(K))
 506         True
 507
 508     """
 509     beta = 0
 510
 511     m = K.dim()
 512     n = K.lattice_dim()
 513     l = lineality(K)
 514
 515     if m < n:
 516         # K is not solid, project onto its span.
 517         K = project_span(K)
 518
 519         # Lemma 2
 520         beta += m*(n - m) + (n - m)**2
 521
 522     if l > 0:
 523         # K is not pointed, project its dual onto its span.
 524         K = project_span(K.dual()).dual()
 525
 526         # Lemma 3
 527         beta += m * l
 528
 529     beta += len(LL(K))
 530     return beta