Semidefinite programming bounds for the average kissing number
Maria Dostert Alexander Kolpakov Fernando Mário de Oliveira Filho
math.MG math.OC 52C17, 90C22, 90C34



  The average kissing number of ^n is the supremum of the average
  degrees of contact graphs of packings of finitely many balls (of any
  radii) in ^n. We provide an upper bound for the average kissing
  number based on semidefinite programming that improves previous
  bounds in dimensions 3, …, 9. A very simple upper bound for
  the average kissing number is twice the kissing number; in
  dimensions 6, …, 9 our new bound is the first to improve on
  this simple upper bound.


P ≈ NP, at least in Visual Question Answering
This work was supported by NVIDIA's NVAIL program and the BMBF project DeFuseNN (Grant 01IW17002)
    

1 Authors contributed equally
Shailza Jolly1,1,2,
Sebastian Palacio1,1,2,
Joachim Folz1,2,
Federico Raue2,
Jörn Hees2,
Andreas Dengel1,2


1DFKI GmbH

2TU Kaiserslautern 

Kaiserslautern, Germany 

<firstname.lastname@dfki.de>


        ^1 Department of Statistical Sciences, University College London, London WC1E 6BT, UK

    ^2 School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK

^3 The Alan Turing Institute, The British Library, London NW1 2DB, UK 


===================================================================================================================================================================================================================================================================

M. Dostert, A. Kolpakov, and F.M. de Oliveira
  FilhoSemidefinite programming bounds for the average kissing
  number





§ INTRODUCTION


A packing of balls in ^n is a finite set of
interior-disjoint closed balls. The contact graph of a
packing  is the graph with vertex set  in which two
balls X and Y are adjacent if they intersect, that is, if they are
tangent to each other.

Contact graphs of packings of disks on the plane are characterized by
the Koebe-Andreev-Thurston theorem <cit.>: they are
precisely the (simple) planar graphs. In higher dimensions, no such
simple characterization is known (see the paper by
Glazyrin <cit.> for a nice discussion), and therefore
research has been focused on understanding the behavior of some
specific parameters of contact graphs.

In this paper, we consider the average degree of contact graphs. More
precisely, we are interested in the average kissing number
of ^n, namely

    _n = sup{ δ(G) : G is the
        contact graph of a packing of balls in ^n },

where δ(G) denotes the average degree of G.

Lower bounds for _n can be obtained by constructions; a
simple idea is to consider lattice packings. Given a
lattice Λ⊆^n with shortest vectors of length d,
we consider the set of all balls of radius d/2 centered on the
lattice points. These balls have disjoint interiors and so we have a
packing of infinitely many balls. Each ball in this packing has the
same number of tangent balls, called the kissing number of
the lattice Λ. The lattice kissing number of ^n,
denoted by τ^*_n, is the largest kissing number of any lattice
in ^n; immediately we have _n ≥τ^*_n.  Conway and
Sloane <cit.> list lower bounds for τ^*_n,
and hence for _n, for n up to 128. For n = 3, a
construction of Eppstein, Kuperberg, and Ziegler <cit.>
gives _3 ≥ 12.612, while τ_3^* = 12.

On the side of upper bounds, it is easy to see
that _n ≤ 2 τ_n, where τ_n is the kissing
  number of ^n, that is, the maximum number of interior-disjoint
unit balls that can simultaneously touch a central unit ball. Indeed,
say  is a packing of balls and let r(X) be the radius of the
ball X ∈; let G = (, E) be the contact graph
of . In G, the number of neighbors of a ball X ∈
that have radius at least r(X) is at most the kissing
number τ_n. So

    |E| ≤∑_X ∈ |{ {X, Y}∈ E : r(X) ≤
      r(Y) }| ≤τ_n ||,

whence the average degree of G is 2|E|/|| ≤ 2 τ_n. Though
simple, this bound is still the best known for all n ≥ 10.

Kuperberg and Schramm <cit.> gave the first nontrivial
upper bound for the average kissing number in dimension 3, proving
that _3 ≤ 8 + 4 √(3) =
14.928…. Glazyrin <cit.> refined their approach and
showed that _3 ≤ 13.955; he also extended their result to
higher dimensions and managed to beat the upper bound of 2τ_n
for n = 4 and 5. In this paper, we use semidefinite programming to
refine Glazyrin's approach (see <ref>
and <ref> below), obtaining better upper bounds
for n = 3, …, 9; see
Table <ref>. In <ref> we discuss an
alternative approach related to the linear programming bound of Cohn
and Elkies <cit.> for the sphere packing density.







 §.§ Notation and preliminaries


The Euclidean inner product on ^n is denoted by
x · y = x_1 y_1 + ⋯ + x_n y_n for x, y ∈^n. The
(n-1)-dimensional unit sphere is
S^n-1 = {  x ∈^n : x = 1 }; the distance between
points x, y ∈ S^n-1 is arccos x · y. The surface
measure on the (n-1)-dimensional sphere of radius ρ is denoted
by ω_ρ; we write ω = ω_1.

A spherical cap in S^n-1 of center x ∈ S^n-1 and
radius α is the set of all points in S^n-1 at distance at
most α from x, namely

    {  y ∈ S^n-1 : x · y ≥cosα }.

The normalized area of this cap is

    ω(S^n-2)/ω(S^n-1)∫_cosα^1 (1 -
      u^2)^(n-3)/2  du.

Spherical caps are defined similarly for spheres of radius other
than 1. Of course, the normalized area of a cap of radius α is
the same irrespective of the radius of the sphere, and is given by the
formula above. The area of a spherical cap can be computed by a
recurrence; see Appendix <ref>.

Let V be a measure space. A kernel is a real-valued 
square-integrable function on V × V. If f V → is
square integrable, then f⊗ f^* is the kernel that
maps (x, y) to f(x) f(y).





§ GLAZYRIN'S UPPER BOUND



Glazyrin <cit.> refines and extends previous work by
Kuperberg and Schramm <cit.> and obtains as a result
the best upper bounds on the average kissing number in dimensions 3,
…, 5. Here is a short description of Glazyrin's approach,
following his presentation. In <ref> we will see how
Glazyrin's bounds can be improved with the help of semidefinite
programming.

Fix ρ > 1 and the dimension n ≥ 3. For r > 0, let B_r be
a ball of radius r tangent to the ball of radius 1 centered at the
origin. The intersection of B_r with the sphere ρ S^n-1 of
radius ρ centered at the origin, if nonempty, is a spherical cap
on ρ S^n-1. The normalized area of this spherical cap is
denoted by A_n,ρ(r), that is,

    A_n,ρ(r) = ω_ρ(B_r ∩ρ
        S^n-1)/ω_ρ(ρ S^n-1),

which as a function of r is monotonically increasing.


  
  If n ≥ 3, ρ > 1, and r > 0,
  then A_n,ρ(r) + A_n,ρ(1/r) ≥ 2 A_n,ρ(1).



If ρ≥ 3, then A_n,ρ(1) = 0 and the result follows, so we
assume ρ < 3. For s > 0, let B_s be a ball of radius s
tangent to the ball of radius 1 centered at the origin. Let us assume
first that both intersections B_r ∩ρ S^n-1
and B_1/r∩ρ S^n-1 are nonempty and hence that both are
spherical caps in ρ S^n-1; let α and β denote
their radii.



Using the law of cosines (see Figure <ref>) we can
determine both cosα and cosβ and as a consequence get

    cosα + cosβ = ρ^2 + 3/2ρ.

Let C denote the right-hand side above, so C ∈ (1, 2)
since 1 < ρ < 3. Write x = cosα;
then cosβ = C - x. Use the formula for the area of a spherical
cap to get

    A_n,ρ(r) + A_n,ρ(1/r) = K ∫_x^1 (1 -
      u^2)^(n-3)/2  du + K ∫_C-x^1 (1 - u^2)^(n-3)/2  du,

where K is a positive constant depending only on the
dimension n.

From the above expression, if f(x) is such
that A_n,ρ(r) + A_n,ρ(1/r) = K f(x), then

    f'(x) = -(1 - x^2)^(n-3)/2 + (1 - (C - x)^2)^(n-3)/2.

So f is monotonically decreasing when 1 - x^2 ≥ 1 - (C-x)^2,
that is, when x ∈ [C-1, C/2], and monotonically increasing
when x ∈ [C/2, 1]. A global minimum of f is therefore attained
at x = C/2, which implies that α = β and so r = 1/r and
therefore r = 1, proving the theorem when
both B_r ∩ρ S^n-1 and B_1/r∩ρ S^n-1 are
nonempty.

Now say B_r ∩ρ S^n-1 is empty;
then B_1/r∩ρ S^n-1 is not empty. Note r < r',
where r' = (ρ - 1) / 2; moreover B_r'∩ρ S^n-1 is a
single point and so A_n,ρ(r') = 0.  Since
also B_1/r'∩ρ S^n-1≠∅, we know
that A_n,ρ(r') + A_n,ρ(1 / r') ≥ 2 A_n,ρ(1) and
hence, since A_n,ρ is monotonically increasing, we get

    A_n,ρ(r) + A_n,ρ(1/r) ≥ A_n,ρ(r') + A_n,ρ(1/r') ≥
      2 A_n,ρ(1),

as we wanted.


Fix ρ > 1 and consider a unit ball centered at the origin. Any
configuration of pairwise interior-disjoint balls (of any radii)
tangent to the central unit ball covers a certain fraction of the
sphere ρ S^n-1 of radius ρ centered at the origin. The
supremum of this covered fraction taken over all possible
configurations is denoted by _n(ρ).


  
  If n ≥ 3 and 1 < ρ < 3, then
  
    _n ≤_n(ρ)/A_n,ρ(1).




Let G = (, E) be the contact graph of a packing  of
balls in ^n. Denote by r(X) the radius of a ball X ∈. On
the one hand, applying Lemma <ref> we get

    ∑_{X, Y}∈ E A_n,ρ(r(X) / r(Y)) + A_n,ρ(r(Y) / r(X))
      ≥ 2 A_n,ρ(1) |E|.

On the other hand, writing N(X) for the set of neighbors
of X ∈, we get

    ∑_{X, Y}∈ E A_n,ρ(r(X) / r(Y)) + A_n,ρ(r(Y) / r(X))
         = ∑_X ∈∑_Y ∈ N(X) A_n,ρ(r(Y) / r(X))
       ≤_n(ρ) ||.


Since ρ < 3, we have A_n,ρ(1) > 0. Putting it all together
we then get

    2|E|/||≤_n(ρ)/A_n,ρ(1),

finishing the proof.


Note that A_n,ρ(1) is simply the normalized area of a spherical
cap of radius α such that

    cosα = ρ^2+3/4ρ,

and so we can compute A_n,ρ(1) explicitly for all n ≥
3. By using in Theorem <ref> the trivial inequality
_n(ρ) ≤ 1 and taking ρ = √(3), we obtain upper
bounds for _n: for n = 3 we get the upper bound
of 14.928… from Kuperberg and Schramm <cit.>;
for n = 4 and 5 we get the upper bounds 34.680…
and 77.756… of Glazyrin <cit.>. The
choice ρ = √(3) is optimal when using the upper
bound _n(ρ) ≤ 1.

Extending techniques of Florian <cit.>,
Glazyrin <cit.> gives a better upper bound
for _3(√(3)), and so obtains _3 ≤ 13.955.





§ REFINING GLAZYRIN'S APPROACH USING SEMIDEFINITE PROGRAMMING



From Theorem <ref> we see that better upper bounds
for _n(ρ) have the potential to give us better upper bounds
for the average kissing number. We will see now how semidefinite
programming can be used to provide upper bounds for _n(ρ);
these upper bounds lead to improved upper bounds for the average
kissing number for n = 3, …, 9 (see Table <ref>).

For fixed 1 < ρ < 3, the function A_n,ρ(r) is increasing
in r and has a limit at infinity, which we denote
by A_n,ρ(∞). It is actually easy to compute this limit:
as r increases, the ball of radius r tangent to the central ball
of radius 1 resembles more and more a hyperplane tangent to the
central ball, so A_n,ρ(∞) is the normalized area
of a spherical cap of radius α such
that cosα = 1 / ρ.



Say two interior-disjoint balls of radii r and s touch a central
unit ball and let x and y be the contact points between each of
the balls and the central ball. Apply the law of cosines (see
Figure <ref>) to get

    x · y ≤1 + r + s - rs/1 + r + s + rs.

Denote the right-hand side above by (r, s).

If F [0,1]^2 → is a kernel and U ⊆ [0, 1] is a
finite set, then the matrix (F(u, v))_u, v ∈ U is a
principal submatrix of F. We denote by P_k^n the Jacobi
polynomial of degree k and parameters α = β = (n-3)/2,
normalized so P_k^n(1) = 1 (see the book by Szegö <cit.>
for background on Jacobi polynomials).


The following theorem is our basic tool to find upper bounds
for _n(ρ).


  
  Let n ≥ 3 be an integer, ρ be such that 1 < ρ < 3,
  and R be such that R > (ρ - 1) / 2. Let r be an increasing
  bijection from [0,1] to [(ρ-1)/2, R] and let
  a [0, 1] → be such that
  a(u) ≥ A_n,ρ(r(u))^1/2 for all u ∈ [0, 1] and
  a(1) ≥ A_n,ρ(∞)^1/2.

  Fix an integer d > 0 and for every k = 0, …, d
  let F_k [0,1]^2 → be a kernel; write
  
    f(t, u, v) = ∑_k=0^d F_k(u, v) P_k^n(t)    for t ∈
          [-1, 1] and u, v ∈ [0,1].

  If f and the kernels F_k are such that
  
  
  (i) every principal submatrix of F_0 - a ⊗ a^* is
    positive semidefinite,

  
  (ii) every principal submatrix of F_k is positive
    semidefinite for k = 0, …, d, and

  
  (iii) f(t, u, v) ≤ 0 whenever -1 ≤ t ≤(r(u), r(v)),
  
  then _n(ρ) ≤max{  f(1, u, u) : u ∈ [0, 1] }.


This theorem is very similar to Theorem 1.2 of de Laat, Oliveira, and
Vallentin <cit.>. They consider configurations of spherical
caps of different radii, but the radii are taken from a finite list of
possibilities; then the function f is matrix valued. Here we work
with configurations of balls of different radii, and the list of
possible radii is infinite, namely the interval [(ρ - 1) / 2,
R]. For this reason we work with the function f as defined in the
theorem above; f can be seen as a kernel-valued function that
assigns to each t ∈ [-1, 1] a kernel on [0,1]^2.


Consider any configuration  of interior-disjoint balls of any
radii tangent to the unit ball centered at the origin. Let Δ be
the normalized area of ρ S^n-1 covered by this configuration
and assume Δ > 0. Since a ball of radius less than (ρ-1) /
2 tangent to the central unit ball does not intersect ρ S^n-1,
we assume that each ball in  has radius at least (ρ - 1) /
2.

Given a ball B ∈, consider the point x ∈ S^n-1 where
it touches the central unit ball. If the radius of B is in the
interval [(ρ-1) / 2, R], then let u ∈ [0, 1] be such
that r(u) is the radius of B; otherwise, set u = 1. Now B will
be represented by the pair (x, u);
let I ⊆ S^n-1× [0, 1] be the set of pairs
representing each ball in .

We will need the following claim: if every principal submatrix
of F [0,1]^2 → is positive semidefinite and if k ≥ 0
is an integer, then the matrix

    ( F(u, v) P_k^n(x · y) )_(x, u), (y, v) ∈ I

is positive semidefinite.

The proof of this claim is as follows. Write

    S    = {  x ∈ S^n-1 : (x, u) ∈ I for some u ∈
          [0,1] }  and
    
        T    = {  u ∈ [0,1] : (x, u) ∈ I for some x ∈
          S^n-1 }.

The addition theorem for Gegenbauer
polynomials <cit.> implies that there is
a real finite-dimensional Hilbert space H and vectors p(x) ∈ H
for x ∈ S such that P_k^n(x · y) = ⟨ p(x), p(y)⟩
for all x, y ∈ S, where ⟨·, ·⟩ denotes the
inner product in H. Similarly, since every principal submatrix
of F is positive semidefinite, there is a real finite-dimensional
Hilbert space, which we may assume to be H as well, and
vectors q(u) ∈ H for u ∈ T such
that F(u, v) = ⟨ q(u), q(v)⟩ for all u, v ∈ T. But
then

    F(u, v) P_k^n(x · y) = ⟨ p(x) ⊗ q(u), p(y) ⊗
      q(v)⟩

for all (x, u), (y, v) ∈ I, and the claim follows.

Since P_0^n is the constant one polynomial, the claim just proved
together with (i), (ii), and the definition of f implies that the
matrix

    (f(x· y, u, v) - a(u) a(v))_(x, u), (y, v) ∈ I

is positive semidefinite. Hence

    ∑_(x, u), (y, v) ∈ I f(x · y, u, v) a(u) a(v) - a(u)^2
      a(v)^2 ≥ 0.


Since  is a configuration of interior-disjoint balls,
if (x, u), (y, v) ∈ I then
-1 ≤ x · y ≤(r(u), r(v)). Now use (iii) and split the
sum on the left-hand side above into the diagonal and off-diagonal
terms to get the inequality

    0    ≤∑_(x, u), (y, v) ∈ I f(x · y, u, v) a(u) a(v)-
        a(u)^2 a(v)^2
       ≤∑_(x, u) ∈ I f(1, u, u) a(u)^2 -
        (∑_(x, u) ∈ I a(u)^2)^2,

so

    ∑_(x, u) ∈ I a(u)^2 ≤max{  f(1, u, u) : u ∈ [0,
      1] }.

Finally, by the construction of I and the properties satisfied
by a we know that Δ is at most the left-hand side above, and
so the theorem follows.


To use Theorem <ref> we need to specify the
kernels F_k. One way to do so is to fix an integer N > 0
and functions p_0, …, p_N [0,1] →. Then, given a
matrix A ∈^(N+1) × (N+1), set

    F(u, v) = ∑_i,j=0^N A_ij p_i(u) p_j(v).

It is easy to check that, if A is positive semidefinite, then every
principal submatrix of F is positive semidefinite. Similarly,
if a = α_0 p_0 + ⋯ + α_N p_N and the matrix
(A_ij - α_i α_j)_i,j=0^N is positive
semidefinite, then every principal submatrix of F - a ⊗ a^* is
positive semidefinite.

In this way, Theorem <ref> can be rephrased as a
semidefinite program. By choosing different
functions p_0, …, p_N, one obtains different optimization
problems, and there is an interplay between the functions chosen to
specify the kernels and the quality of the approximation a
of u ↦ A_n,ρ(r(u))^1/2 that one can obtain. In the next
two sections we will use this approach to construct optimization
problems that give bounds for _n(ρ) that lead to the new
upper bounds for the average kissing number in Table <ref>;
for the functions p we will take alternately step functions and
polynomials.





 §.§ Step functions


Let us first set the p_i to be step functions.
Fix R > (ρ - 1)/2 and let r [0, 1] → [(ρ-1)/2, R] be
such that

    r(u) = (R - (ρ - 1) / 2) u + (ρ-1)/2.

Note that r is an increasing bijection between [0, 1]
and [(ρ - 1) / 2, R].

Now fix an integer N > 0 and
points 0 = s_0 < s_1 < ⋯ < s_N < s_N+1 = 1.
Let S_i = [s_i, s_i+1) for i = 0, …, N - 1
and S_N = [s_N, s_N+1]. Let p_i be the function that is 1
on S_i and 0 everywhere else.

The function a is now simple to specify: for u ∈ [0, 1] set

    a(u) = 
        A_n,ρ(r(s_i+1))^1/2   if u ∈ S_i for some i <
          N;
    
        A_n,ρ(∞)^1/2   if u ∈ S_N.

Then a is an upper approximation of the
function u ↦ A_n,ρ(r(u))^1/2 (since this is a
monotonically increasing function) as needed in
Theorem <ref>, and moreover a is a linear
combination of the p_i functions.

Each kernel F_k is parameterized by an (N+1) × (N+1)
positive-semidefinite matrix A_k as follows:

    F_k(u, v) = ∑_i,j=0^N A_k,ij p_i(u) p_j(v)    for
        all u, v ∈ [0, 1],

where A_k,ij = (A_k)_ij.  So the kernels F_k are constant on
the sets S_i × S_j, and therefore F_k can be quite naturally
identified with A_k. For fixed u, v ∈ [0, 1], the
function t↦ f(t, u, v) defined in (<ref>) is a
polynomial on t; for i, j = 0, …, N we write f_ij(t)
for the common value that f(t, u, v) assumes on S_i × S_j,
that is,

    f_ij(t) = f(t, s_i, s_j) = ∑_k=0^d ∑_i,j=0^N A_k,ij
      p_i(s_i) p_j(s_j) P_k^n(t) = ∑_k=0^d A_k,ij P_k^n(t).


Say that u ∈ S_i and v ∈ S_j; from (<ref>) we
have (r(u), r(v)) ≤(r(s_i), r(s_j)). So to ensure that f
satisfies item (iii) of Theorem <ref> we have to
ensure that, for all i, j = 0, …, N,

    f_ij(t) ≤ 0    whenever -1 ≤ t ≤(r(s_i),
        r(s_j)).


Summarizing, if a = α_0 p_0 + ⋯ + α_N p_N, then any
feasible solution of the following optimization problem gives an upper
bound for _n(ρ):

    min   2lmax{ f_ii(1) : i = 0, …, N }
       2lf_ij(t) = ∑_k=0^d A_k,ij P_k^n(t),
       f_ij(t) ≤ 0   whenever -1 ≤ t ≤(r(s_i),
          r(s_j)),
       (A_0,ij - α_i α_j)_i,j=0^N   is
          positive semidefinite,
       A_k ∈^(N+1) × (N+1)   is positive semidefinite for k =
          0, …, d.


To model the nonpositivity constraints on the functions f_ij we
ensure nonpositivity on a finite sample of points
in [-1, (r(s_i), r(s_j))]. For the complete approach and a
description of how the solutions found by a solver can be rigorously
verified, see Appendix <ref>.

Problem <ref> can be used to give better bounds for the
average kissing number in dimensions 5, …, 9; see
Table <ref>. In dimension 3, we could not beat Glazyrin's
bound using this problem; in dimension 4, the bound provided is
better than Glazyrin's bound, but worse than the bound
of <ref> below. To obtain the bounds of
Table <ref>, we used ρ = 2, N = 30,
and R ≈ 184.25; see the verification script
(cf. Appendix <ref>) for precise information.





 §.§ Polynomials



We now take the functions p_i to be polynomials.  Fix an
integer N > 0 and let p_i(u) = u^i for i = 0, …, N.
Fix R > (ρ - 1) / 2 and
let r [0, 1] → [(ρ - 1) / 2, R] be such that

    r(u) = (R - (ρ - 1) / 2) u^2 + (ρ - 1) / 2.

Note that r is an increasing bijection between [0, 1]
and [(ρ - 1) / 2, R]. Note also that, in comparison
with (<ref>), we have u^2 instead of u above; we could
use u, but we noticed that this leads to a worse upper
approximation a.

To get the function a we solve a simple linear program
that is set up as follows. We have variables α_0,
…, α_N for the coefficients of p_0, …, p_N. We fix
some ϵ≥ 0 and a finite sample S of points in [0, 1]
and consider the constraints

    [ α_0 p_0(u) + ⋯ + α_N p_N(u) ≥
        A_n,ρ(r(u))^1/2 + ϵ                                            for u ∈ S,;       α_0 p_0(1) + ⋯ + α_N p_N(1) ≥
        A_n,ρ(∞)^1/2. ]

The objective is to minimize

    max{ α_0 p_0(u) + ⋯ + α_N p_N(u) -
      A_n,ρ(r(u))^1/2 : u ∈ S }.




By taking a fine enough sample of points and a small
positive ϵ≈ 10^-5, the
function a = α_0 p_0 + ⋯ + α_N p_N will satisfy the
properties in Theorem <ref>, though this has to
be verified (see Appendix <ref>). If N is large
enough, the function a should also be close enough
to u ↦ A_n,ρ(r(u))^1/2 (see
Figure <ref>).

Each kernel F_k is parameterized by an (N+1) × (N+1)
positive-semidefinite matrix A_k as follows:

    F_k(u, v) = ∑_i,j=0^N A_k,ij p_i(u) p_j(v)    for
        all u, v ∈ [0, 1].

So the function f defined in (<ref>) is a polynomial
on t, u, and v. Recalling (<ref>), item (iii) in
Theorem <ref> then asks that f(t, u, v) should
be nonpositive on the semialgebraic set { (t, u, v) : s_i(t, u, v)
≥ 0 for i = 1, …, 4 }, where the s_i are the
following polynomials:

    s_1    = t + 1,
    
        s_2    = 1 + r(u) + r(v) - r(u) r(v) - t (1 + r(u) + r(v) + r(u) r(v)),
    
        s_3    = u(1-u) + v(1-v),  and
    
        s_4    = u(1-u) v(1-v);

note that s_3(t, u, v) ≥ 0 and s_4(t, u, v) ≥ 0 if and only
if 0 ≤ u ≤ 1 and 0 ≤ v ≤ 1.[Why note
  take s_3 = u(1-u) and s_4 = v(1-v)? Then it is also true
  that s_3(t, u, v) ≥ 0 and s_4(t, u, v) ≥ 0 if and only
  if 0 ≤ u ≤ 1 and 0 ≤ v ≤ 1. The reason for our
  choice is that the polynomials s_3 and s_4
  in (<ref>) are symmetric in u and v, and this
  helps us reduce the size of the corresponding semidefinite program;
  see Appendix <ref>.]

The nonpositivity condition can then be restricted to a sum-of-squares
condition: we require that there exist polynomials q_1,
…, q_5, each a sum of squares, such that

    f = -s_1 q_1 - s_2 q_2 - s_3 q_3 - s_4 q_4 - q_5,

since then f is clearly nonpositive in the required domain.

Finally, the upper bound on _n(ρ) is given by the maximum
value of the function u ↦ f(1, u, u) on [0, 1]; note that
this function is a univariate polynomial on u of degree 2N. A
theorem of Lukács <cit.> says that this
maximum is equal to the minimum λ for which there are
univariate polynomials l_1 and l_2, each a sum of squares, such
that

    λ - f(1, u, u) = l_1(u) + u (1 - u) l_2(u).


Putting it all together, we get the following optimization problem, any
feasible solution of which gives an upper bound for _n(ρ):

    min   λ
       f(t, u, v) = ∑_k=0^d ∑_i,j=0^N A_k,ij p_i(u) p_j(v)
        P_k^n(t),
       q_1, …, q_5 are sum-of-squares polynomials satisfying (<ref>),
       l_1 and l_2 are sum-of-squares polynomials
          satisfying (<ref>),
       (A_0,ij - α_i α_j)_i,j=0^N is
          positive semidefinite,
       A_k ∈^(N+1) × (N+1) is positive semidefinite for k =
          0, …, d.


Problem (<ref>) can be rewritten as a
semidefinite program, where each sum-of-squares polynomial
is parameterized by a positive-semidefinite
matrix. Appendix <ref> gives a detailed description
of the semidefinite program we solve and an overview of
how the solution found by the solver can be verified to be
feasible.

The approach of this section provides better bounds for the average
kissing number in dimensions 3 and 4; see
Table <ref>. In higher dimensions we could not manage to
obtain any bounds using this approach, since the problems are
infeasible when polynomials of low degree are used, and too large when
polynomials of high degree are used.





§ A DIRECT LINEAR PROGRAMMING BOUND



Suppose we want to find the average kissing number, but that we
restrict ourselves to packings having balls of a few prescribed radii,
say r_1 < ⋯ < r_N. How many neighbors can a vertex in the
contact graph of such a packing have? Or, in other words, how many
balls can touch a given ball in the packing? Certainly, the largest
number of balls touching a central ball is attained when the central
ball has the largest possible radius, r_N, and every ball touching
it has the smallest possible radius, r_1. Hence the maximum degree
of the contact graph, and by consequence its average degree, is at
most the maximum number of interior-disjoint balls of radius r_1
that can simultaneously touch a central ball of radius r_N.

This is a simple upper bound for this restricted average kissing
number, but one could object it is too local: the bound does
not take the whole packing into account, ignoring the interaction
between different balls. In particular, it is usually impossible for
every vertex in the packing to have the maximum possible degree, since
not every vertex can be the largest ball surrounded by several small
balls.

The bound for the average kissing number given by
Theorem <ref> appears to be similarly local. It is
based on the parameter _n(ρ), which is not defined in terms
of a packing of balls and therefore cannot take into account the
interaction between different balls in a packing. We discuss now an
alternative idea, based on the linear programming bound of Cohn and
Elkies <cit.> for the sphere packing density, which seems to
overcome this issue.

A continuous (matrix-valued)
function f^n →^N × N is of positive
  type if for every finite set U ⊆^n the block matrix

    (f(x-y))_x,y ∈ U

is positive semidefinite. Matrix-valued functions of positive type are
straightforward extensions of functions of positive type; see e.g. the
paper by de Laat, Oliveira, and Vallentin <cit.> for
more on such functions.


  
  Let r_1, …, r_N be any positive numbers. If f^n
  →^N × N is a continuous function of positive type such
  that
  
  
  (i) f(x)_ij≤ 0 if x≥ r_i + r_j and
  
  (ii) f(x)_ij≤ -1 if x = r_i + r_j,
  
  then the average degree of the contact graph of a packing of balls
  of radii r_1, …, r_N is at most max{  f(0)_ii :
  i = 1, …, N }.



Let  be a packing of balls of radii r_1, …, r_N and
let I ⊆{1, …, N}×^n be such
that (i, x) ∈ I if and only if  has a ball of radius r_i
centered at x. Since f is of positive type we know that

    ∑_(i, x), (j, y) ∈ I f(x-y)_ij≥ 0.


Let G = (, E) be the contact graph of the packing .
Split the sum above into three parts: the diagonal terms, the terms
corresponding to pairs of balls that do not touch, and the terms
corresponding to pairs of balls that do touch. Since f satisfies (i)
and (ii) we get

    0    ≤∑_(i, x) ∈ I f(0)_ii + ∑_(i, x),
            (j, y) ∈ I
    r_i+r_j ≠x-y  f(x-y)_ij
        + ∑_(i, x),
            (j, y) ∈ I
    r_i+r_j = x-y  f(x-y)_ij
       ≤∑_(i, x) ∈ I f(0)_ii - 2|E|
       ≤ |I| max{ f(0)_ii : i = 1, …, N } - 2|E|,

whence 2|E| / || ≤max{ f(0)_ii : i = 1,
  …, N }, and the theorem follows.


This theorem gives a direct bound for the average kissing number
instead of the rather indirect bound of Theorem <ref>
via the parameter _n(ρ). Moreover, we have a
two-point bound, that takes into account interactions between
pairs of balls in a packing.

Though Theorem <ref> is stated for a finite number of
possible radii, it can be easily extended to account for radii in any
bounded interval [a, b] if one uses kernel-valued functions. It is
not immediately clear, however, how to adapt the theorem for packings
of balls of arbitrary radii, and so Theorem <ref>
cannot be directly used to compute upper bounds for the average
kissing number.

For any given radii r_1, …, r_N, it is possible to reduce the
problem of finding functions satisfying the conditions in
Theorem <ref> to a semidefinite program;
such a reduction was employed by de Laat, Oliveira, and
Vallentin <cit.> for a similar problem. In this way,
concrete bounds can be computed.

The bound of Theorem <ref> can be adapted to packings
of balls of finitely many possible radii r_1, …, r_N, namely
by changing the definition of _n(ρ). Fix k = 1, …, N
and ρ > 1; consider the unit ball centered at the origin. Any
configuration of pairwise interior-disjoint balls of
radii r_i / r_k, for i = 1, …, N, covers a certain fraction
of the sphere ρ S^n-1 of radius ρ centered at the
origin. Denote the supremum of this covered fraction, taken over all
such configurations, by _n^k(ρ), and let

    _n(ρ) = max{_n^1(ρ), …, _n^N(ρ)}.

Following the proof of Theorem <ref>, we see that the
average degree of the contact graph of any packing of balls of
radii r_1, …, r_N is at most

    _n(ρ)/A_n,ρ(1).

Upper bounds for _n^k(ρ) can be computed using the approach
of de Laat, Oliveira, and Vallentin <cit.>.

So it is possible to compare numerically, for different choices of
radii, bounds given by Theorem <ref> with bounds given
by Theorem <ref>. The case N = 1 is particularly
simple, since then _n(ρ) is the kissing number τ_n
of ^n times the area covered by the spherical cap, which
is A_n,ρ(1). So

    _n(ρ)/A_n,ρ(1) = τ_n,

and moreover any upper bound for the kissing number, like the
linear programming bound of Delsarte, Goethals, and
Seidel <cit.>, gives an upper bound for the average
degree of the contact graph. As for Theorem <ref>, we
have observed numerically that it provides a worse upper bound than
the linear programming bound for the kissing number. Surprisingly,
when more than one radius is considered, the bound of
Theorem <ref> becomes even worse;
Table <ref> contains some results.







§ ACKNOWLEDGEMENTS


We would like to thank the Complex Systems and Big Data Competence
Centre at the University of Neuchâtel for providing access to their
computational cluster “Cervino”.



ain-bib.tex











§ COMPUTING THE AREA OF A SPHERICAL CAP



To use Theorem <ref> we need to be able to
compute the normalized area of a spherical cap of radius α
in S^n-1, which is given by

    ω(S^n-2)/ω(S^n-1)∫_cosα^1 (1 -
      u^2)^(n-3)/2  du.

The factor before the integral can be computed, to any desired
precision, by means of a recurrence. We will now derive a recurrence
relation for the integral above, making it possible to compute the
normalized area to any desired precision.

Let s be a real number. The Taylor series of u ↦ (1 - u)^s
around 0 is

    ∑_k=0^∞ (-1)^k s (s-1) ⋯ (s - k + 1) u^k/k!
      = ∑_k=0^∞ (-s)_k u^k/k!,

where for a real number a and an integer k ≥ 0 we denote
by (a)_k the shifted factorial:

    (a)_k = 
        1,   if k = 0;
    
        a(a+1)⋯(a+k-1),   otherwise.


Substitute u by u^2 and integrate term-by-term to get

    ∫ (1 - u^2)^s  du = ∑_k=0^∞ (-s)_k
      u^2k+1/(2k+1) k! = u _2F_1(1/2, -s; 3/2; u^2),

where _2F_1 is the hypergeometric
series <cit.>.

So we want to compute

    F_n(u) = _2F_1(1/2, -(n - 3) / 2; 3 / 2; u^2).

Equation (15.2.11) in the book by Abramovitz and
Stegun <cit.> gives us the relation

    (c - b) _2F_1(a, b-1; c; z) + (2b - c - bz + az) _2F_1(a, b;
      c; z)
     + b(z-1) _2F_1(a, b+1; c; z) = 0.

Take a = 1/2, b = -(n - 3) / 2, and c = 3/2 to get an
expression for F_n+2(u) in terms of F_n(u) and F_n-2(u). It
is now easy to obtain a recurrence; the base cases are

    F_3(u)    = 1,   
        F_4(u)    = u (1 - u^2)^1/2 + arcsin u/2u,
    
        F_5(u)    = 1 - u^2 / 3,    and   
        F_6(u)    = u (5 - 2u^2) (1 - u^2)^1/2 + 3arcsin u/8u.






§ THE SEMIDEFINITE PROGRAM FOR STEP FUNCTIONS AND RIGOROUS
  VERIFICATION



For each i, j = 0, …, N, to implement the nonpositivity
constraint for f_ij in problem (<ref>) we select a
finite sample _ij of points in [-1, (r(s_i),
r(s_j))]. Consider the matrix W such
that W_ij = α_i α_j and fix ϵ > 0. To obtain an
upper bound we solve the following semidefinite program,
in which the role of the A_0 variable changes in comparison
with (<ref>):

    min   max{ f_ii(1) : i = 0, …, N }
       f_ij(t) = W_ij + ∑_k=0^d A_k,ij P_k^n(t),
       f_ij(t) ≤ -ϵ   for t ∈_ij,
       A_k ∈^(N+1) × (N+1) is positive semidefinite for k =
          0, …, d.


In practice, we select samples of 50 points for each i, j and
set ϵ = 10^-5. We solve the resulting problem with standard
solvers and obtain a tentative optimal value z^*.  The next step is
to remove the objective function and add it as a constraint, requiring
that

    max{ f_ii(1) : i = 0, …, N }≤ z^* +
      η,

where η≈ 10^-3. When we solve this feasibility problem,
the solver returns a strictly feasible solution, that is, a solution
in which every matrix A_k is positive definite. We observed that
this solution immediately satisfies the original nonpositivity
constraints of (<ref>).

To verify that we have indeed a feasible solution, we only have to
verify that each A_k is positive semidefinite and compute an upper
bound on the value of f_ij on [-1, (r(s_i), r(s_j))]. Since
each A_k is actually positive definite, we use high-precision
floating-point arithmetic to compute for each A_k its Cholesky
decomposition L_k. Then we replace A_k by L_k L_k^, so A_k
becomes positive semidefinite by construction.

To get an upper bound for the value of f_ij on the corresponding
interval, we use interval arithmetic. We split the original interval
into subintervals and evaluate f_ij on each subinterval, obtaining
for each subinterval an upper bound on the value of f_ij. In this
way, we obtain an upper bound u_ij on the value of f_ij on the
original interval.

Since the definition of f_ij uses W, which in practice is
computed numerically, it is not enough to have u_ij≤
0. Indeed, if we use the exact value of W_ij instead of an
approximation, then f_ij could change to a positive number. To
prevent this from happening, we need to ensure that u_ij is
negative enough compared to the absolute error in the computation
of W_ij. If we use the formulas of Appendix <ref> to
compute W_ij using high-precision interval arithmetic, then we
have a rigorous bound on the absolute error of each W_ij. This
whole verification approach is implemented in a Sage <cit.>
script included with the arXiv version of this paper.





§ THE SEMIDEFINITE PROGRAM FOR POLYNOMIAL INTERPOLATION AND 
  RIGOROUS VERIFICATION


Two steps are required in order to obtain rigorous upper bounds for
the average kissing number via the approach
of <ref>. First, we must find a polynomial that
approximates the spherical-cap-area function from above, and we must
prove that this polynomial is really an upper bound for this
function. Second, we must rewrite problem (<ref>) as a
semidefinite program, find good solutions for it, and prove that they
are feasible. In this section, we will see how both steps can be
carried out.





 §.§ Verifying the approximation for A_n,ρ



In <ref> we have seen how a polynomial a
satisfying the conditions described in
Theorem <ref> can be found. Here we quickly
describe how it can be rigorously verified that a given polynomial a
satisfies these conditions; this verification approach is implemented
in a Sage <cit.> script included with the arXiv version of this
paper.

Say ρ and R are fixed and let r(u) be defined as
in (<ref>). We want to prove
that a(u)^2 - A_n,ρ(r(u)) ≥ 0 for all u ∈ [0, 1] and
moreover that a(1)^2 - A_n,ρ(∞) ≥ 0; the difficulty
lies in testing the validity of the first set of conditions.

Say we have an upper bound M on the absolute value of the derivative
of the function f(u) = a(u)^2 - A_n,ρ(r(u)) on the
interval [0, 1]. For an integer N ≥ 1, write ϵ = 1 / (N
+ 1) and consider the points u_k = kϵ for k = 0,
…, N+1. If η is the minimum of f on the points u_k,
then the mean-value theorem implies that

    f(u) ≥η - ϵ M / 2    for all u ∈ [0, 1].

So, as long as η - ϵ M / 2 ≥ 0, the function f is
nonnegative on [0, 1]. Computing the function f on lots of points
yields a guess for η; then we find N such
that ϵ = 1 / (N+1) ≤ 2η/M and try to test the function
on the points u_k. This is the approach implemented in the
verification script.

It remains to see how to compute the upper bound M on the absolute
value of the derivative of f. Since a(u), and hence a(u)^2, is a
polynomial, it is easy to compute an upper bound on the absolute value
of the derivative of a(u)^2 rigorously using interval
arithmetic.

Bounding the derivative of A_n,ρ(r(u)) is also simple. Indeed,
recall that A_n,ρ(s) equals

    ω(S^n-2)/ω(S^n-1)∫_x^1 (1 -
      z^2)^(n-3)/2  dz,

where

    x = ρ^2 + 2s + 1/2ρ(1+s).

So

    A_n,ρ'(s) = -ω(S^n-2)/ω(S^n-1)
      (1-x^2)^(n-3)/2ρ^2-1/2ρ(1+s)^2.

Note that 0 ≤ x ≤ 1, so |(1-x^2)^(n-3)/2| ≤ 1. The
rightmost fraction above is largest when s = 0 and
|r'(u)| = |2 (R - (ρ-1)/2) u| ≤ 2(R - (ρ-1)/2)
for u ∈ [0,1]. Finally, apply the chain rule to get

    |d A_n,ρ(r(u))/du| ≤ω(S^n-2)/ω(S^n-1)|ρ^2-1|/ρ (R -
      (ρ-1)/2)

for all u ∈ [0, 1].





 §.§ The semidefinite program and how to verify feasibility



We quickly discuss how problem (<ref>) is transformed
into a semidefinite program and then how a provably
feasible solution can be found for it.

To transform (<ref>) into a semidefinite program, one has
to encode each sum-of-squares polynomial in terms of
positive-semidefinite matrices; here is the well-known
recipe. Let P_1, …, P_m be a basis of the
space [x_1, …, x_n]_≤ k of n-variable real polynomials
of degree at most k and for x = (x_1, …, x_n) let v_k(x) be
the vector such that

    (v_k(x))_i = P_i(x)    for i = 1, …, m.

We can see v_k as a “vector” whose entries are the
polynomials P_i.

A polynomial p of degree 2k is a sum of squares if and only
if there is a positive-semidefinite matrix X such that

    p(x) = ⟨ v_k(x) v_k(x)^, X⟩,

where ⟨ A, B⟩ =  AB is the trace inner product
between symmetric matrices A and B. Note that above we have a
polynomial identity: both the left and right-hand sides are
polynomials that we require to be equal. This polynomial identity can
be rewritten as a set of linear constraints on the entries of the
matrix X: for each monomial of degree at most 2k we have one
constraint relating the coefficient of the monomial on both left and
right-hand sides.

The polynomials in (<ref>) are univariate and have
degree 2N. So to model this constraint we
use v_N(u) = (1, u, u^2, …, u^N).

Since the polynomial f(t, u, v) is symmetric on u and v, that
is, f(t, u, v) = f(t, v, u), and since the same holds for the
polynomials s_1, …, s_4 in (<ref>), we are able
to further reduce the sizes of the positive-semidefinite matrices
needed to encode (<ref>).

Indeed, say that a sum-of-squares polynomial p ∈[t, u, v] is
symmetric on u and v. It is not necessarily true that there is a
sum-of-squares decomposition p = q_1^2 + ⋯ + q_m^2 of p in
which every q_i is symmetric. However, we have that

    [t, u, v] = [t, u+v, uv] ⊕ (u - v) [t,
      u+v, uv],

where [t, u+v, uv] is the ring of polynomials symmetric on u
and v. So say q_i = a_i + (u-v) b_i, where a_i and b_i are
symmetric. Then

    p(t, u, v)    = (1/2)(p(t, u, v) + p(t, v, u))
       =1/2∑_i=1^m a_i(t, u, v)^2 + 2(u-v)a_i(t, u,
        v)b_i(t, u, v) + (u - v)^2 b_i(t, u, v)^2
       =1/2∑_i=1^m   + a_i(t, v, u)^2 + 2(v-u)a_i(t, v,
        u)b_i(t, v, u) + (v - u)^2 b_i(t, v, u)^2
       =∑_i=1^m a_i(t, u, v)^2 + (u-v)^2 b_i(t, u, v)^2.

For any given k, let v_k(t, u, v) be obtained from a basis
of [t, u+v, uv]_≤ k as before. The above discussion implies
that p ∈[t, u, v] of degree 2k is a symmetric sum of squares
if and only if there are positive-semidefinite matrices X and X'
such that

    p = ⟨ v_k v_k^, X⟩ + ⟨ (u-v)^2 v_k-1
      v_k-1^, X'⟩.


Now it should be clear how (<ref>) can be rewritten as a
semidefinite program. Only one technical detail remains,
namely how to determine the degrees of the sum-of-squares polynomials
in (<ref>). Here we use for each polynomial the smallest
possible degree such that no term in the right-hand side
of (<ref>) has degree larger than d + 2N, which is the
degree of f.

The bounds for n = 3 and 4 in Table <ref> were obtained
by solving the semidefinite program described above with d = 10
and N = 6, 8, respectively. To solve the problem we use the
high-precision solver SDPA-GMP <cit.>, but even so the
solution found by the solver is not truly feasible. To extract a
rational feasible solution from it we use the Julia library developed
by Dostert, de Laat, and Moustrou <cit.>. This allows us
to provide a rigorous upper bound for the average kissing number. The
script to generate the semidefinite program and to obtain an exact
rational solution is included with the arXiv version of this paper.


