mjo/cone/rearrangement.py

   1 from sage.all import *
   2
   3 def rearrangement_cone(p,n,lattice=None):
   4     r"""
   5     Return the rearrangement cone of order ``p`` in ``n`` dimensions.
   6
   7     The rearrangement cone in ``n`` dimensions has as its elements
   8     vectors of length ``n``. For inclusion in the cone, the smallest
   9     ``p`` components of a vector must sum to a nonnegative number.
  10
  11     For example, the rearrangement cone of order ``p == 1`` has its
  12     single smallest component nonnegative. This implies that all
  13     components are nonnegative, and that therefore the rearrangement
  14     cone of order one is the nonnegative orthant.
  15
  16     When ``p == n``, the sum of all components of a vector must be
  17     nonnegative for inclusion in the cone. That is, the cone is a
  18     half-space in ``n`` dimensions.
  19
  20     INPUT:
  21
  22       - ``p`` -- The number of components to "rearrange."
  23
  24       - ``n`` -- The dimension of the ambient space for the resulting cone.
  25
  26       - ``lattice`` -- (default: ``None``) an ambient lattice of rank ``n``
  27                        to be passed to the :func:`Cone` constructor.
  28
  29     OUTPUT:
  30
  31     A polyhedral closed convex cone object representing a rearrangement
  32     cone of order ``p`` in ``n`` dimensions. Each generating ray will
  33     have the integer ring as its base ring.
  34
  35     If a ``lattice`` was specified, then the resulting cone will live in
  36     that lattice unless its rank is incompatible with the dimension
  37     ``n`` (in which case a ``ValueError`` is raised).
  38
  39     ALGORITHM:
  40
  41     The generators for the rearrangement cone are given by [Jeong]_
  42     Theorem 5.2.3.
  43
  44     REFERENCES:
  45
  46     .. [GowdaJeong] Muddappa Seetharama Gowda and Juyoung Jeong.
  47        Spectral cones in Euclidean Jordan algebras.
  48        Linear Algebra and its Applications, 509, 286-305.
  49        doi:10.1016/j.laa.2016.08.004
  50
  51     .. [HenrionSeeger] Rene Henrion and Alberto Seeger.
  52        Inradius and Circumradius of Various Convex Cones Arising in
  53        Applications. Set-Valued and Variational Analysis, 18(3-4),
  54        483-511, 2010. doi:10.1007/s11228-010-0150-z
  55
  56     .. [Jeong] Juyoung Jeong.
  57        Spectral sets and functions on Euclidean Jordan algebras.
  58        University of Maryland, Baltimore County, Ph.D. thesis, 2017.
  59
  60     SETUP::
  61
  62         sage: from mjo.cone.rearrangement import rearrangement_cone
  63
  64     EXAMPLES:
  65
  66     The rearrangement cones of order one are nonnegative orthants::
  67
  68         sage: rearrangement_cone(1,1) == Cone([(1,)])
  69         True
  70         sage: rearrangement_cone(1,2) == Cone([(0,1),(1,0)])
  71         True
  72         sage: rearrangement_cone(1,3) == Cone([(0,0,1),(0,1,0),(1,0,0)])
  73         True
  74
  75     When ``p == n``, the resulting cone will be a half-space, so we
  76     expect its lineality to be one less than ``n`` because it will
  77     contain a hyperplane but is not the entire space::
  78
  79         sage: rearrangement_cone(5,5).lineality()
  80         4
  81
  82     All rearrangement cones are proper when ``p`` is less than ``n`` by
  83     [Jeong]_ Proposition 5.2.1::
  84
  85         sage: all( rearrangement_cone(p,n).is_proper()
  86         ....:              for n in range(10)
  87         ....:              for p in range(1, n) )
  88         True
  89
  90     The Lyapunov rank of the rearrangement cone of order ``p`` in ``n``
  91     dimensions is ``n`` for ``p == 1`` or ``p == n`` and one otherwise,
  92     by [Jeong]_ Corollary 5.2.4::
  93
  94         sage: all( rearrangement_cone(p,n).lyapunov_rank() == n
  95         ....:              for n in range(2, 10)
  96         ....:              for p in [1, n-1] )
  97         True
  98         sage: all( rearrangement_cone(p,n).lyapunov_rank() == 1
  99         ....:              for n in range(3, 10)
 100         ....:              for p in range(2, n-1) )
 101         True
 102
 103     TESTS:
 104
 105     All rearrangement cones are permutation-invariant by [Jeong]_
 106     Proposition 5.2.1::
 107
 108         sage: n = ZZ.random_element(2,10).abs()
 109         sage: p = ZZ.random_element(1,n)
 110         sage: K = rearrangement_cone(p,n)
 111         sage: P = SymmetricGroup(n).random_element().matrix()
 112         sage: all( K.contains(P*r) for r in K )
 113         True
 114
 115     The smallest ``p`` components of every element of the rearrangement
 116     cone should sum to a nonnegative number (this tests that the
 117     generators really are what we think they are)::
 118
 119         sage: set_random_seed()
 120         sage: def _has_rearrangement_property(v,p):
 121         ....:     return sum( sorted(v)[0:p] ) >= 0
 122         sage: all( _has_rearrangement_property(
 123         ....:      rearrangement_cone(p,n).random_element(),
 124         ....:      p
 125         ....:    )
 126         ....:    for n in range(2, 10)
 127         ....:    for p in range(1, n-1)
 128         ....: )
 129         True
 130
 131     The rearrangenent cone of order ``p`` is contained in the rearrangement
 132     cone of order ``p + 1`` by [Jeong]_ Proposition 5.2.1::
 133
 134         sage: set_random_seed()
 135         sage: n = ZZ.random_element(2,10)
 136         sage: p = ZZ.random_element(1,n)
 137         sage: K1 = rearrangement_cone(p,n)
 138         sage: K2 = rearrangement_cone(p+1,n)
 139         sage: all( x in K2 for x in K1 )
 140         True
 141
 142     The rearrangement cone of order ``p`` is linearly isomorphic to the
 143     rearrangement cone of order ``n - p`` when ``p`` is less than ``n``,
 144     by [Jeong]_ Proposition 5.2.1::
 145
 146         sage: set_random_seed()
 147         sage: n = ZZ.random_element(2,10)
 148         sage: p = ZZ.random_element(1,n)
 149         sage: K1 = rearrangement_cone(p,n)
 150         sage: K2 = rearrangement_cone(n-p, n)
 151         sage: Mp = (1/p)*matrix.ones(QQ,n) - identity_matrix(QQ,n)
 152         sage: Cone( (Mp*K2.rays()).columns() ).is_equivalent(K1)
 153         True
 154
 155     The order ``p`` should be between ``1`` and ``n``, inclusive::
 156
 157         sage: rearrangement_cone(0,3)
 158         Traceback (most recent call last):
 159         ...
 160         ValueError: order p=0 should be between 1 and n=3, inclusive
 161         sage: rearrangement_cone(5,3)
 162         Traceback (most recent call last):
 163         ...
 164         ValueError: order p=5 should be between 1 and n=3, inclusive
 165
 166     If a ``lattice`` was given, it is actually used::
 167
 168         sage: L = ToricLattice(3, 'M')
 169         sage: rearrangement_cone(2, 3, lattice=L)
 170         3-d cone in 3-d lattice M
 171
 172     Unless the rank of the lattice disagrees with ``n``::
 173
 174         sage: L = ToricLattice(1, 'M')
 175         sage: rearrangement_cone(2, 3, lattice=L)
 176         Traceback (most recent call last):
 177         ...
 178         ValueError: lattice rank=1 and dimension n=3 are incompatible
 179
 180     """
 181     if p < 1 or p > n:
 182         raise ValueError('order p=%d should be between 1 and n=%d, inclusive'
 183                          %
 184                          (p,n))
 185
 186     if lattice is None:
 187         lattice = ToricLattice(n)
 188
 189     if lattice.rank() != n:
 190         raise ValueError('lattice rank=%d and dimension n=%d are incompatible'
 191                          %
 192                          (lattice.rank(), n))
 193
 194     I = identity_matrix(ZZ,n)
 195     M = matrix.ones(ZZ,n) - p*I
 196     G = identity_matrix(ZZ,n).rows() + M.rows()
 197     return Cone(G, lattice=lattice)