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