mjo/cone/schur.py

   1 """
   2 The Schur cone, as described in the "Critical angles..." papers by
   3 Iusem, Seeger, and Sossa. It defines the Schur ordering on `R^{n}`.
   4 """
   5
   6 from sage.all import *
   7
   8 def schur_cone(n):
   9     r"""
  10     Return the Schur cone in ``n`` dimensions.
  11
  12     INPUT:
  13
  14     - ``n`` -- the dimension the ambient space.
  15
  16     OUTPUT:
  17
  18     A rational closed convex Schur cone of dimension ``n``.
  19
  20     REFERENCES:
  21
  22     .. [GourionSeeger] Daniel Gourion and Alberto Seeger.
  23        Critical angles in polyhedral convex cones: numerical and
  24        statistical considerations. Mathematical Programming, 123:173-198,
  25        2010, doi:10.1007/s10107-009-0317-2.
  26
  27     .. [IusemSeegerOnPairs] Alfredo Iusem and Alberto Seeger.
  28        On pairs of vectors achieving the maximal angle of a convex cone.
  29        Mathematical Programming, 104(2-3):501-523, 2005,
  30        doi:10.1007/s10107-005-0626-z.
  31
  32     .. [SeegerSossaI] Alberto Seeger and David Sossa.
  33        Critical angles between two convex cones I. General theory.
  34        TOP, 24(1):44-65, 2016, doi:10.1007/s11750-015-0375-y.
  35
  36     SETUP::
  37
  38         sage: from mjo.cone.nonnegative_orthant import nonnegative_orthant
  39         sage: from mjo.cone.schur import schur_cone
  40
  41     EXAMPLES:
  42
  43     Verify the claim that the maximal angle between any two generators
  44     of the Schur cone and the nonnegative quintant is ``3*pi/4``::
  45
  46         sage: P = schur_cone(5)
  47         sage: Q = nonnegative_orthant(5)
  48         sage: G = ( g.change_ring(QQbar).normalized() for g in P )
  49         sage: H = ( h.change_ring(QQbar).normalized() for h in Q )
  50         sage: actual = max(arccos(u.inner_product(v)) for u in G for v in H)
  51         sage: expected = 3*pi/4
  52         sage: abs(actual - expected).n() < 1e-12
  53         True
  54
  55     The dual of the Schur cone is the "downward monotonic cone"
  56     [GourionSeeger]_, whose elements' entries are in non-increasing
  57     order::
  58
  59         sage: set_random_seed()
  60         sage: n = ZZ.random_element(10)
  61         sage: K = schur_cone(n).dual()
  62         sage: x = K.random_element()
  63         sage: all( x[i] >= x[i+1] for i in xrange(n-1) )
  64         True
  65
  66     TESTS:
  67
  68     We get the trivial cone when ``n`` is zero::
  69
  70         sage: schur_cone(0).is_trivial()
  71         True
  72
  73     The Schur cone induces the majorization ordering::
  74
  75         sage: set_random_seed()
  76         sage: def majorized_by(x,y):
  77         ....:     return (all(sum(x[0:i]) <= sum(y[0:i])
  78         ....:                 for i in xrange(x.degree()-1))
  79         ....:             and sum(x) == sum(y))
  80         sage: n = ZZ.random_element(10)
  81         sage: V = VectorSpace(QQ, n)
  82         sage: S = schur_cone(n)
  83         sage: majorized_by(V.zero(), S.random_element())
  84         True
  85         sage: x = V.random_element()
  86         sage: y = V.random_element()
  87         sage: majorized_by(x,y) == ( (y-x) in S )
  88         True
  89
  90     """
  91
  92     def _f(i,j):
  93         if i == j:
  94             return 1
  95         elif j - i == 1:
  96             return -1
  97         else:
  98             return 0
  99
 100     # The "max" below catches the trivial case where n == 0.
 101     S = matrix(ZZ, max(0,n-1), n, _f)
 102
 103     # Likewise, when n == 0, we need to specify the lattice.
 104     return Cone(S.rows(), ToricLattice(n))