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     Try the formula in my paper::
 208
 209         sage: K = random_cone(max_dim=15, max_rays=25)
 210         sage: actual = lyapunov_rank(K)
 211         sage: K_S = project_span(K)
 212         sage: J_T1 = project_span(K, K_S.dual())
 213         sage: J_T2 = project_span(K_S.dual()).dual()
 214         sage: J_T2 = Cone(J_T2.rays(), lattice=J_T1.lattice())
 215         sage: J_T1 == J_T2
 216         True
 217         sage: J_T = J_T1
 218         sage: l = K.linear_subspace().dimension()
 219         sage: codim = K.lattice_dim() - K.dim()
 220         sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2
 221         sage: actual == expected
 222         True
 223
 224     """
 225     V = K.lattice().vector_space()
 226
 227     C_of_K = discrete_complementarity_set(K)
 228
 229     tensor_products = [s.tensor_product(x) for (x,s) in C_of_K]
 230
 231     # Sage doesn't think matrices are vectors, so we have to convert
 232     # our matrices to vectors explicitly before we can figure out how
 233     # many are linearly-indepenedent.
 234     #
 235     # The space W has the same base ring as V, but dimension
 236     # dim(V)^2. So it has the same dimension as the space of linear
 237     # transformations on V. In other words, it's just the right size
 238     # to create an isomorphism between it and our matrices.
 239     W = VectorSpace(V.base_ring(), V.dimension()**2)
 240
 241     # Turn our matrices into long vectors...
 242     vectors = [ W(m.list()) for m in tensor_products ]
 243
 244     # Vector space representation of Lyapunov-like matrices
 245     # (i.e. vec(L) where L is Luapunov-like).
 246     LL_vector = W.span(vectors).complement()
 247
 248     # Now construct an ambient MatrixSpace in which to stick our
 249     # transformations.
 250     M = MatrixSpace(V.base_ring(), V.dimension())
 251
 252     matrix_basis = [ M(v.list()) for v in LL_vector.basis() ]
 253
 254     return matrix_basis
 255
 256
 257
 258 def lyapunov_rank(K):
 259     r"""
 260     Compute the Lyapunov (or bilinearity) rank of this cone.
 261
 262     The Lyapunov rank of a cone can be thought of in (mainly) two ways:
 263
 264     1. The dimension of the Lie algebra of the automorphism group of the
 265        cone.
 266
 267     2. The dimension of the linear space of all Lyapunov-like
 268        transformations on the cone.
 269
 270     INPUT:
 271
 272     A closed, convex polyhedral cone.
 273
 274     OUTPUT:
 275
 276     An integer representing the Lyapunov rank of the cone. If the
 277     dimension of the ambient vector space is `n`, then the Lyapunov rank
 278     will be between `1` and `n` inclusive; however a rank of `n-1` is
 279     not possible (see the first reference).
 280
 281     .. note::
 282
 283         In the references, the cones are always assumed to be proper. We
 284         do not impose this restriction.
 285
 286     .. seealso::
 287
 288         :meth:`is_proper`
 289
 290     ALGORITHM:
 291
 292     The codimension formula from the second reference is used. We find
 293     all pairs `(x,s)` in the complementarity set of `K` such that `x`
 294     and `s` are rays of our cone. It is known that these vectors are
 295     sufficient to apply the codimension formula. Once we have all such
 296     pairs, we "brute force" the codimension formula by finding all
 297     linearly-independent `xs^{T}`.
 298
 299     REFERENCES:
 300
 301     .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
 302        cone and Lyapunov-like transformations, Mathematical Programming, 147
 303        (2014) 155-170.
 304
 305     .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
 306        optimality constraints for the cone of positive polynomials,
 307        Mathematical Programming, Series B, 129 (2011) 5-31.
 308
 309     EXAMPLES:
 310
 311     The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
 312     [Rudolf et al.]_::
 313
 314         sage: positives = Cone([(1,)])
 315         sage: lyapunov_rank(positives)
 316         1
 317         sage: quadrant = Cone([(1,0), (0,1)])
 318         sage: lyapunov_rank(quadrant)
 319         2
 320        sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 321         sage: lyapunov_rank(octant)
 322         3
 323
 324     The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
 325     [Rudolf et al.]_::
 326
 327         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 328         sage: lyapunov_rank(L31)
 329         1
 330
 331     Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
 332
 333         sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
 334         sage: lyapunov_rank(L3infty)
 335         1
 336
 337     The Lyapunov rank should be additive on a product of cones
 338     [Rudolf et al.]_::
 339
 340         sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
 341         sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
 342         sage: K = L31.cartesian_product(octant)
 343         sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
 344         True
 345
 346     Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
 347     The cone ``K`` in the following example is isomorphic to the nonnegative
 348     octant in `\mathbb{R}^{3}`::
 349
 350         sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
 351         sage: lyapunov_rank(K)
 352         3
 353
 354     The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
 355     itself [Rudolf et al.]_::
 356
 357         sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
 358         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 359         True
 360
 361     TESTS:
 362
 363     The Lyapunov rank should be additive on a product of cones
 364     [Rudolf et al.]_::
 365
 366         sage: K1 = random_cone(max_dim=10, max_rays=10)
 367         sage: K2 = random_cone(max_dim=10, max_rays=10)
 368         sage: K = K1.cartesian_product(K2)
 369         sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
 370         True
 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 = random_cone(max_dim=10, max_rays=10)
 376         sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
 377         True
 378
 379     The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
 380     be any number between `1` and `n` inclusive, excluding `n-1`
 381     [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
 382     trivial cone in a trivial space as well. However, in zero dimensions,
 383     the Lyapunov rank of the trivial cone will be zero::
 384
 385         sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True)
 386         sage: b = lyapunov_rank(K)
 387         sage: n = K.lattice_dim()
 388         sage: (n == 0 or 1 <= b) and b <= n
 389         True
 390         sage: b == n-1
 391         False
 392
 393     """
 394     return len(LL(K))