Philippe G Lefloch

Jean-Marc Mercier
email: jean-marc.mercier@mpg-partners.com

Mesh-free error integration in arbitrary dimensions: a numerical study of discrepancy functions

We are interested in mesh-free formulas based on the Monte-Carlo methodology for the approximation of multi-dimensional integrals, and we investigate their accuracy when the functions belong to a reproducing-kernel space. A kernel typically captures regularity and qualitative properties of functions "beyond" the standard Sobolev regularity class. We are interested in the issue whether quantitative error bounds can be a priori guaranteed in applications (including mathematical finance, scientific computing, and machine learning). Our main contribution is a numerical study of the error discrepancy function based on a comparison between several numerical strategies, when one varies the choice of the kernel, the number of approximation points, and the dimension of the problem. We consider two strategies in order to localize to a bounded set the standard kernels defined in the whole Euclidian space (exponential, multiquadric, Gaussian, truncated), namely, on one hand the class of periodic kernels defined via a discrete Fourier transform on a lattice and, on the other hand, a class of transport-based kernels. Relying on the Poisson formula on a lattice together with heuristic arguments, we study the derivation of theoretical bounds for the discrepancy function of periodic kernels. Then, for four kernels of particular interest we perform extensive numerical experiments and generate the optimal distributions of points and the discrepancy error functions. Our numerical results validate our theoretical observations and, importantly, provide us with quantitative estimates for the error made with a kernel-based strategy as opposed to a purely random strategy.

Introduction

An error approximation formula. We are motivated here by applications to partial differential equations arising continuum physics, including the development of mesh-free methods in fluid dynamics and material sciences [START_REF] Bessa | A meshfree unification: reproducing kernel peridynamics[END_REF][START_REF] Günther | Implementation of boundary conditions for mesh-less methods[END_REF][START_REF] Koester | Conforming window functions for mesh-free methods[END_REF][START_REF] Li | Mesh-free particle methods[END_REF][START_REF] Liu | Mesh-free methods: moving beyond the finite element method[END_REF][START_REF] Nakano | Convergence of meshfree collocation methods for fully nonlinear parabolic equations[END_REF][START_REF] Oh | Mesh-free particle methods for thin plates[END_REF][START_REF] Zhou | Solving phase field equations using a mesh-less method[END_REF]. Specifically, we are interested in approximating multi-dimensional integrals via Monte-Carlo formula and deriving error estimates, in which the dependency with respect to the dimension of the problem and other important parameters is specified in a quantitative manner. By revisiting this problem of multivariate integration, our purpose is to clarify the derivation and validity of such estimates whose importance has been highlighted in recent years in artificial intelligence, for mesh-free computations of partial differential equations, and mathematical finance. The existing literature emphasizes the role of Sobolev-type spaces, while we would like here to stress the importance of kernel-based Hilbert spaces. In many applications, one is interested in preserving certain a priori structure that are available a priori and the choice of a kernel is dictated by properties (symmetry, scaling, regularity, decay, etc.) that should be incorporated in the approximation algorithms. Therefore, it is desirable to have a flexible framework that encompasses a wide class of kernels, as we consider in the present paper.

More specifically, within a given Hilbert (or Banach) space we seek to optimize the choice of the interpolation points in an integral approximation formula and establish a sharp error estimate within the chosen class of regularity and decay. Two parameters are of primary interest, namely, the dimension D ≥ 1 of the problem and the number of interpolation points N ≥ 1, and it is essential to have quantitative estimates with a specified dependency in N, D that can be determined from the kernels of interest. In the present paper, we contribute to this general objective and provide a systematic study and comparison of several classes of kernels, which we refer to as periodic kernels and transported kernels. Our periodic framework for periodic kernels is motivated by work by Cohn and Elkies [START_REF] Cohn | New upper bounds on sphere packings[END_REF] who studied the problem of sphere packing. Our result depends upon a "kernel density" function which arises as a key factor in a quantitative bound. We build here on many earlier works on the subject, including contributions in approximation theory [START_REF] Babuska | Survey of mesh-less and generalized finite element methods: a unified approach[END_REF][START_REF] Fasshauer | Mesh-free methods[END_REF][START_REF] Fasshauer | Mesh-free approximation methods with Matlab[END_REF][START_REF] Matsumoto | Mersenne twister: a 623-dimensionally equidistributed uniform pseudorandom number generator[END_REF][START_REF] Opfer | Multiscale kernels[END_REF][START_REF] Salehi | A moving least square reproducing polynomial mesh-less method[END_REF][START_REF] Wendland | Sobolev-type error estimates for interpolation by radial basis functions[END_REF][START_REF] Wendland | Scattered data approximation[END_REF].

The discrepancy function associated with a kernel. To any kernel K : Ω × Ω → R defined on a bounded and open subset Ω ⊂ R D and satisfying a positivity condition (see Section 2.1), we associate a Banach space H s,p K (Ω) of real-valued functions defined on Ω with regularity exponent s > 0 and integrability exponent p ∈ [1, +∞). (More generally, the Lebesgue measure could be replaced by a probability measure.) Then, an "abstract" error integration estimate reads, for any N ≥ 1 and any function ϕ ∈ H s,p K (Ω), sup

ϕ H s,p K (Ω) =1 inf x 1 ,...,x N ∈Ω 1 |Ω| Ω ϕ(x) dx - 1 N 1≤n≤N ϕ(x n ) ≤ E s,p K (N, D), (1.1)
in which the discrepancy function E s,p K (N, D) is independent of ϕ. Hence, (1.1) provides us -in the class of functions under consideration-a factorization of the error in two contributions: ϕ H s,p K (Ω) measures the regularity of the function while the discrepancy function is related to the best distribution of N points in Ω. The challenge is to control E s,p K (N, D), which can be expressed in several forms: 1. In the physical space Ω, the function E s,p K (N, D) can be formulated with a pseudo-distance associated with the kernel.

2. In suitable spectral variables determined from an operator naturally associated with the kernel, the function E s,p K (N, D) takes a rather explicit form involving the eigenfunctions and eigenvectors of this operator.

However, both formulations are difficult to work with directly -except in dimension D = 1. So, we introduce below a third standpoint which is more efficient in order to control and its dependency with respect to D and N , that is, we introduce the class of "lattice-based" kernels (in a tensorial form), as we call them. In this context, we can express the function E s,p K (N, D) via: 3. a Poisson formula in dual discrete Fourier variables associated with a lattice (see next section).

Interestingly, for this latter class of kernels, quantitative estimates can be established that involve the the notion of a "lattice density" function, as we explain it in this paper, and shed some light on the problem of the curse of dimensionality. A priori and quantitative error bounds are obtained at any order of accuracy at the expense of possibly increasing the regularity of the functions under consideration. Importantly for the applications, the error function is controlled quantitatively in a given functional framework.

Evaluation of the discrepancy function. We focus attention on a selected list of kernels which we construct by a nonlinear transformation of four translation-invariant kernels, that is, of the form K(x, y) = ϕ(x -y). We choose kernels that are commonly used in the applications, namely the exponential 1 , multiquadric, Gaussian, and truncated kernels. Their Fourier transform ϕ defined on R D is known explicitly and is listed in Table 1.1. (see [START_REF] Wendland | Scattered data approximation[END_REF])

Table 1.1: Four kernels and their Fourier transforms on R D

These kernels are defined in R D and we proceed by "localizing" them to a bounded domain Ω, taken to be the unit cube [0, 1] D for simplicity in the presentation. We propose two methods for such a localization of a kernel K defined on R D × R D :

• The periodic version K per of K defined from a discrete Fourier transform.

• The transported version K tran of K defined via a nonlinear transport map. This provides us with eight kernels (listed in Table 1.2) and our main purpose is to investigate the discrepancy function associated with each of them.

In principle, we could use numerically any one of the three expressions of the error function which we derive below and attempt to minimize it over the set of Y . In most cases, this requires a computation which, in general, cannot be done explicitly and a numerical integration of this function would be very costly, especially in large dimensions. We discuss this below. In particular, due to the (non-convex) form of the kernel, minimizing the error function in the physical space is computationally challenging. By introducing a periodic version based on the discrete Fourier transform, we arrive at an expression that is computationally tractable. We are able to make comparisons between these kernels and investigate the rate of convergence while comparing with the case when the points are randomly chosen. Our numerical results confirm and support our theoretical discussion.

Applications and perspectives. The material in this paper should be useful for analyzing meshfree methods for computing solutions to partial differential equations and deriving quantitative bounds for algorithms used in pattern recognition and artificial intelligence. The estimate discussed here provides us with a key building block in order to establish an error analysis of the transported mesh-free method presented in the companion paper [START_REF] Lefloch | The Transport-based Mesh-free Method (TMM) and its applications in mathematical finance and fluid dynamics[END_REF]. The method therein can be regarded as a generalization of the Lagrangian mesh-free method use in computational fluid dynamics, but also allows to include Navier-Stokes-type diffusive terms.

Most of the literature on error integration estimates is focused on functions with Sobolev regularity while we are interested here in functions with regularity adapted to specific applications. For instance, the standard choice of radially-symmetric kernels leads to functional spaces that are variants of Sobolev spaces and, in particular, are invariant by translations. Allowing more general kernels allows one to describe local (direction-dependent) properties of functions. For instance, a kernel we discuss below is adapted to measure the regularity of functions of the form ϕ = 0<n1<...<n k ≤D ϕ n1,...,n k (x n1 , . . . , x n k ), relevant in mathematical finance. The strategy in [START_REF] Lefloch | A new method for solving Kolmogorov equations in mathematical finance[END_REF][START_REF] Lefloch | The Transport-based Mesh-free Method (TMM) and its applications in mathematical finance and fluid dynamics[END_REF] is now applied in industrial applications [START_REF] Mercier | Hedging strategies for net interest income and economic values of equity[END_REF] and its accuracy can be explained in the light of the present study. This is relevant when considering the valuations of complex financial products (including the so-called American exercising) written on a large number of underlyings, and aiming at computing rapidly complex risk measures; see [START_REF] Lefloch | A new method for solving Kolmogorov equations in mathematical finance[END_REF].

Outline of this paper. In Section 2, we present some basic material on reproducing kernel spaces. In Section 3, we discuss our methodology for constructing the two classes of kernels of main interest. In Section 4, the kernels studied in the present paper are presented and some their properties discussed. In Section 5, we derive several expressions of the discrepancy function, depending whether physical, spectral, or Fourier variables are used and, next, in Section 6 we derive estimates on the discrepancy error. In Section 7, we present and discuss our numerical results for each of the kernels of interest. 2 Functional framework based on a reproducing kernel

exponential (E) multiquadric (M) Gaussian (G) truncated (T) Tra exp(-|erf -1 (x)| 1 ) 1 + |erf -1 (x)| 2 -(D+1)/2 exp(-|erf -1 (x)| 2 /2) sup 1 -|erf -1 (x)|, 0 D Per exp(-|x + α| 1 ) (1 + |x + α| 2 ) -(D+1)/2 exp(-|x + α| 2 /2) sup(1 -|x + α|, 0) D

Discrete setup

The class of admissible kernels. Since we are primarily interested in kernels defined on a bounded set, in the present section we restrict attention to this class -although we will allow ourselves to manipulate kernels defined on the whole Euclidian space and treated as "seed data" in order to generate the kernels of actual interest. A reproducing kernel provides a convenient way to generate a broad class of Hilbert spaces (or, more generally, Banach space); cf. [START_REF] Fasshauer | Mesh-free methods[END_REF][START_REF] Wendland | Sobolev-type error estimates for interpolation by radial basis functions[END_REF]. A bounded and continuous function

K : Ω × Ω → R on a bounded open set Ω ⊂ R D is called an admissible kernel if it satisfies (1)
the symmetry property: K(x, y) = K(y, x) for all x, y ∈ Ω, and (2) the positivity property: for any collection of N distinct points Y = (y 1 , . . . , y N ) in Ω, the symmetric matrix

K(Y, Y ) = K(y m , y n ) 1≤n,m≤N is positive definite in the sense that a T K(Y, Y )a > 0 for all a ∈ R N \ {0}.
It is said to be uniformly positive if there exists a uniform constant c > 0 such that for any collection of distinct points Y one has

a T K(Y, Y )a ≥ c |a| 2 for all a ∈ R N .
Clearly, any admissible kernel also satisfies

K(x, x) ≥ 0, K(x, y) 2 ≤ K(x, x) K(y, y), x, y ∈ Ω. ( 2.1)
This implies that 2 K(x, y) ≤ K(x, x) + K(y, y) and, therefore, the non-negative function

D(x, y) := K(x, x) + K(y, y) -2K(x, y) ≥ 0, x, y ∈ Ω, (2.2)
can be interpreted as a "pseudo-distance" in view of the properties D(x, x) = 0 and D(x, y) = D(y, x).

(The triangle inequality need not hold.) Many examples of admissible kernels will be presented in the next two sections. Finite dimensional framework. Given any finite collection of points Y = (y 1 , . . . , y N ) chosen in Ω, we introduce the (finite dimensional) vector space H Y K (Ω) consisting of all linear combinations of the basis functions x → K(x, y n ). In other words, we set

H Y K (Ω) := 1≤m≤N a m K(•, y m ) / a = (a 1 , . . . , a N ) ∈ R N . (2.3) Since K is continuous, H Y K (Ω) ⊂ C(Ω)
embeds into the space C(Ω) of all continuous functions on Ω. To any two functions ϕ = 1≤m≤N a m K(•, y m ) and ψ = 1≤n≤N b n K(•, y n ), we associate the bilinear expression

ϕ, ψ H Y K (Ω) := a T K(Y, Y )b = 1≤m≤N 1≤n≤N a m b n K(y m , y n ) (2.4)
(with a = (a m ), etc.), which endows the space H Y K (Ω) with a Hilbertian structure with norm ϕ 2

H Y K (Ω) := a T K(Y, Y )a. Now, the so-called reproducing kernel property (immediate from (2.4)) K(•, y m ), K(•, y n ) H Y K (Ω) = K(y m , y n ), (2.5)
allows one to relate the coefficients of the decomposition of a function ϕ = 1≤m≤N a m K(•, y m ) to its scalar product with the basis functions, namely

ϕ, K(•, y n ) H Y K (Ω) = 1≤m≤N a m K(•, y m ), K(•, y n ) H Y K (Ω) = 1≤m≤N a m K(y m , y n ) = a T K(Y, y n ). (2.6)
Discrete spectral decomposition. Since K(Y, Y ) is a symmetric and positive definite matrix, it admits real and positive eigenvalues, denoted by λ n Y > 0, together with a basis of right-eigenvectors ζ n Y ∈ R D (with n = 1, 2, . . . , N ) chosen to be unitary, and satisfying

K(Y, Y )ζ n Y = λ n Y ζ n Y , n = 1, 2, . . . , N. (2.7)
This decomposition is useful in order to define, for any s ≥ 0 and p ≥ 1, the finite dimensional Banach space H Y,s,p K (Ω) of all functions of the form 1≤m≤N a m K(•, y m ) with finite norm

ϕ p H Y,s,p K (Ω) := 1≤n≤N (λ n Y ) -sp ϕ, ζ n Y p H Y K (Ω) .
(2.8)

Projection operator. Consider a function f ∈ C(Ω) and introduce the vector f (Y ) = f (y 1 ), . . . , f (y N ) consisting of the values of this function at the given points. We define its projection P Y (f ) into the discrete space H Y K (Ω) by setting

P Y (f ) := a T K(•, Y ), a := K(Y, Y ) -1 f (Y ). (2.9a)
Clearly, this defines a projection since P Y • P Y (f ) = P Y (f ) and, in fact, P Y (ϕ) = ϕ for any function ϕ = a T K(•, Y ) belonging to the space H Y K (Ω) (where we used that K is symmetric). Moreover, the norm of this projection reads

P Y (f ) 2 H Y K (Ω) = f (Y ) T K(Y, Y ) -1 f (Y ).
(2.9b)

The partition of unity. A basis is naturally associated with the discrete space H Y K (Ω), that is, N functions θ n Y : Ω → R taking the values 0 or 1 at the points of the set Y . Precisely, writing

θ Y := (θ 1 Y , . . . , θ N Y ), we define θ Y := K(Y, Y ) -1 K(Y, •). (2.10a) It follows that θ Y (y m ) 1≤n,m≤N = K(Y, Y ) -1 K(Y, Y ) =
Id (the identity matrix), and using the Kronecker symbol, we have θ n Y (y m ) = δ nm , while the scalar product of any two basis functions is

θ m Y , θ n Y H Y K (Ω) = K -1 (y m , y n ). (2.10b)
This partition of unity is useful in expressing the projection of a function f ∈ C(Ω), namely

P Y (f ) = 1≤n≤N f (y n )θ n Y .

Continuous setup

Functional spaces. Given an admissible kernel K : Ω × Ω → R, we now introduce the infinite dimensional space H K (Ω) consisting of all finite linear combinations of the functions K(x, •) parametrized by x ∈ Ω, that is, H K (Ω) := Span K(•, x) / x ∈ Ω , which we endow with the scalar product and norm defined in the finite dimensional setup; see (2.1) (where now Y and N are no longer fixed). By construction, the reproducing kernel property (2.5) also holds in H K (Ω), i.e.

K(•, x), K(•, y) H K (Ω) = K(x, y), x, y ∈ Ω. (2.11a)
From the Cauchy-Schwarz inequality, it follows that for any ϕ ∈ H K (Ω) and x ∈ Ω

|ϕ(x)| = | K(•, x), ϕ H K (Ω) | ≤ K(•, x) H K (Ω) ϕ H K (Ω) = K(x, x) ϕ H K (Ω) . (2.11b)
Since the kernel is continuous and bounded, the "point evaluation" ϕ → ϕ(x) is thus a linear and bounded functional on H K (Ω) (for any x ∈ Ω).

We have defined a pre-Hilbert space, that is, a vector space endowed with a scalar product and, in order to obtain a complete metric space, the completion of the pre-Hilbert space H K (Ω) is considered by taking all linear combinations based on countably many points Y = (y 1 , y 2 , . . .) in Ω. The corresponding space is denoted by H K (Ω) and is the reproducing Hilbert space generated from the kernel K.

Clearly, both properties (2.11) remain true in H K (Ω). Observe also that

|ϕ(x) -ϕ(y)| = | K(•, x) -K(•, y), ϕ H K (Ω) | K(•, x) -K(•, y) H K (Ω) ϕ H K (Ω) = D(x, y) ϕ H K (Ω) ,
x, y ∈ Ω.

(2.12)

Since K and thus D are continuous in Ω, we have the embedding H K (Ω) ⊂ C(Ω), that is, all of the functions are continuous, at least. Mercer decomposition. We now consider the linear operator T K : L2 (Ω) → L 2 (Ω) defined by T K (a) := Ω K(•, x)a(x) dx (for a ∈ L 2 (Ω)) on the Hilbert space L 2 (Ω) endowed with its standard inner product. We have

T K (a) L 2 (Ω) ≤ |Ω| 1/2 K L 2 (Ω×Ω) a L 2 (Ω) ,
so this operator is continuous and self-adjoint. It is easily checked to be compact: if a sequence a p a weakly in the L 2 norm then T K (a p ) → T K (a) strongly in the L 2 norm. The classical spectral decomposition applies and the operator T K admits an (at most countable) non-increasing sequence of eigenvalues λ i > 0 and a corresponding set of eigenfunctions ζ i such that T K (ζ i ) = λ i ζ i and the family ζ 1 , ζ 2 , . . . forms an orthonormal basis of L 2 (Ω) for the L 2 inner product. Furthermore, since the kernel is continuous and bounded, the eigenfunctions ζ i are continuous, at least.

We then introduce the kernel

L(x, y) = j=1,2,... λ j ζ j (x)ζ j (y), x, y ∈ Ω, (2.13)
in which the sum converges in the L 2 sense. This kernel is admissible since, for any collection of points Y = (y 1 , . . . , y N ) in R N , the bilinear form

L Y (x, x ) = 1≤n,m≤N L(y n , y m )x n x m (with x, x ∈ R N ), satisfies L Y (x, x) = j=1,2,... λ j 1≤n,m≤N ζ j (y n )ζ j (y m )x n x m = j=1,2,... λ j 1≤n≤N ζ j (y n )x n 2 ≥ 0.
In fact, one can check that L coincides with the given kernel K and (2.13) represents its spectral decomposition. The family of functions λ 1/2 i ζ i i≥1 is an orthonormal basis of the space H K (Ω), as follows from the defining relation

T K ζ i = λ i ζ i , namely 2 < λ 1/2 i ζ i , λ 1/2 j ζ j > H K (Ω) = λ 1/2 i λ 1/2 j 1 λ i < T K ζ i , ζ j > H K (Ω) = λ 1/2 j λ 1/2 i ζ i , ζ j L 2 (Ω) = δ ij . (2.14)
This series representation of the elements in H K (Ω) is referred to as the Mercer representation (which may be non-unique). In short, Mercer theorem states that any admissible kernel K : Ω × Ω → R can be viewed as a positive, self-adjoint and compact operator on L 2 (Ω). Banach spaces. Based on the Mercer representation, we can define the L p -based spaces at any order of differentiability. Namely, for any s ≥ 0 and p ∈ [1, +∞) we consider the Banach space

H s,p K (Ω) := ϕ ∈ L p (Ω) / i=1,2,... λ -ps/2 i ϕ, ζ i L 2 (Ω) p < +∞ (2.15)
endowed with the norm ϕ H s,p K (Ω)

p := i=1,2,... λ -ps/2 i ϕ, ζ i L 2 (Ω) p .
When p is chosen to be 2, the space H s K (Ω) := H s,2 K (Ω) is a Hilbert space endowed with the inner product

f, g H s K (Ω) = i=1,2,... (λ i ) -s f, ζ i L 2 (Ω) g, ζ i L 2 (Ω) . (2.16)
In particular,

H K (Ω) = H 1 K (Ω) = H 1,2 K (Ω)
and we recover the Hilbert space defined earlier.

3 Methodology for defining classes of kernels Translation invariant kernels. Our first task now is to provide some preliminary material and a classification involving the notions of translation-invariant kernel, periodic kernel, tensorial kernel, and radiallysymmetric kernel. In the present section we allow ourselves to introduce kernels defined on the whole of R D , although Section 2 was restricted to kernels defined on a bounded set; namely, we will use kernel defined on R D only for defining the kernels of interest defined on a bounded set. We begin with the class of the form K(x, y) = χ(x -y) (the examples in Table 1.1 being of this type), which we refer to as translation-invariant kernels

K(x, y) = χ(x -y), x, y ∈ R D . (3.1a)
This family is parametrized by a generating function χ : R D → R, which (after normalization) must satisfy (as we check below)

χ(0) = 1, χ(-x) = χ(x), x ∈ R D , χ(ξ) ≥ 0 ξ ∈ R D . ( 3
(Y, Y )(ξ, ξ) = 1≤n,m≤N K(y n , y m )ξ n ξ m for ξ ∈ R N . We obtain the positivity property 1≤n,m≤N χ, e -i<y n -y m ,•> D ,D ξ n ξ m = χ, n ξ n e -i<y n ,•> 2 D ,D ≥ 0, in which •, • D ,D
is the duality braket for distributions. We thus find

R D ×R D K(x, y)ϕ(x)ϕ(y) dxdy = R D ×R D χ(x -y)ϕ(x)ϕ(y) dxdy = R D (χ ϕ)(x)ϕ(x) dx = R D (χ ϕ)(ξ) ϕ(ξ) dξ = R D | ϕ| 2 χ dξ ≥ 0, or equivalently R D ×R D K(x, y)ϕ(x)ϕ(y) dxdy ≥ 0.
Radially-symmetric kernels. Among translation-invariant kernels, the class of radially-symmetric kernels is an important subclass and corresponds to the case where the function χ that only depends on the modulus of its argument only, that is,

K(x, y) = χ(|x -y|), x, y ∈ R D , (3.2)
where the generating function χ : R → R is now regarded as a function of a real variable.

Localization via scaling: the transported kernels

The translation-invariant fail to be sufficiently localized (say in the sense that K fails to be L 1 (R D × R D ).), and we now explain how a transportation map can be applied in order to "localize" a kernel to a bounded set. Many proofs will be postponed to the Appendix, including the one of the following result.

Proposition 3.1. Let K seed : R D × R D → R be a bounded admissible kernel and µ be a sufficiently regular, probability measure such that supp(µ) is a convex set (with, therefore, µ ≥ 0 and µ(R D ) = 1). Then the kernel

K loc (x, x ) := K seed (x, x )µ(x)µ(x ), x, x ∈ R D , (3.3a
)

is admissible and belongs to L 1 (R D × R D ). Moreover, let Ω ⊂ R D be
any open and convex subset with normalized volume |Ω| = 1, and consider a transportation map for the measure µ, that is, a one-to-one map S : Ω → supp(µ) such that S = ∇h for some convex function h : Ω → R and S # µ = dy (the Lebesgue measure). Then, K tran (y, y ) := K seed (S(y), S(y )), y, y ∈ Ω (3.3b) defines an admissible kernel, referred to as the transported kernel associated with K seed and µ. Furthermore, for any ϕ ∈ L 1 µ (R D ) (the set of functions that are integrable for the measure µ) and any choice of points x n = S(y n ) one has

R D ϕ(x) dµ(x) - 1 N 1≤n≤N ϕ(x n ) = Ω (ϕ • S)(y) dy - 1 N 1≤n≤N (ϕ • S)(y n ). (3.3c)
In the statement above, the transportation maps satisfies

R D ϕdµ = Ω (ϕ • S)dy for any continuous function ϕ ∈ L 1 µ (R D ). Provided S is sufficiently regular, one has S # µ = | det ∇S| • S -1 dy,
where | det ∇S| is the Jacobian of S. The convexity of supp(µ) ensures that S is one-to-one and continuous from Ω onto the support set supp(µ). Furthermore, thanks to the localization argument above, deriving an error estimate associated with the left-hand side of (3.3c) reduces to (1.1), and the role of the measure µ is eliminated. Within the framework of kernel spaces, the relation between the µ-weighted norm on R D and the unweighted norm on the bounded set Ω is as follows (with obvious notation):

ϕ H s,p K seed ,µ (R D ) = ϕ H K loc,s,p (R D ) = ϕ • S H K tran,s,p (Ω) .

Localization via periodization: the periodic kernels

A discrete lattice. Motivated by Cohn and Elkies's work [START_REF] Cohn | New upper bounds on sphere packings[END_REF] on the sphere packing problem, we propose to embed the support of a general kernel in a periodic lattice. By periodicity, we can always extend to R D the kernel defined to an elementary cell. Importantly, the terms arising in the spectral decomposition can be controlled almost explicitly, thanks to a Poisson decomposition formula associated with the lattice.

A family of D vectors l 1 , l 2 , . . . , l D ∈ R D being given, we consider their convex hull C ⊂ R D which serves as the fundamental cell of our discrete lattice and whose volume is denoted by |C|. By suitably translating C, we thus generate the periodic lattice L := 1≤d≤D α d l d α = (α 1 , . . . , α D ) ∈ Z D , and we denote its dual by L * := α * ∈ R D < α, α * >∈ Z for all α ∈ L . We denote by l * 1 , . . . , l * D the vectors generating the elementary cell C * ⊂ R D of the dual lattice. Next, we introduce the discrete Hilbert space l 2 (L * ) consisting of all functions defined on the dual lattice, endowed with the inner product

f, g l 2 (L * ) := 1 |C| α * ∈L * f (α * )g(α * ), f, g : L * → R. (3.4)
Poisson formula. Consider the Fourier transform of a real-valued function ϕ : R D → R defined on R D and let us restrict it to the dual lattice, that is, consider the discrete values

ϕ(α * ) := R D ϕ(x) e -2iπ<x,α * > dx. (3.5)
Then, the so-called Poisson formula reads

α∈L ϕ(α + x) = 1 |C| α * ∈L * e 2iπ<x,α * > ϕ(α * ), x ∈ R D . (3.6)
Provided the function ϕ is supported on the cell C, the sum in the left-hand side contains a single term and, therefore,

ϕ(x) = 1 |C| α * ∈L * e 2iπ<x,α * > ϕ(α * ), x ∈ C, provided supp(ϕ) ⊂ C (3.7a)
or, with our notation,

ϕ(x) = e 2iπ<x,•> , ϕ l 2 (L * ) , x ∈ C, provided supp(ϕ) ⊂ C. (3.7b)
Hence, a function defined on the cell C can be recovered (via an discrete inverse Fourier transform) from the values of its Fourier transform on the dual lattice L * . More generally, let us derive an identity that will be useful to us later on. Consider now a collection (y 1 , . . . , y N ) of points in R D and, in view of (3.7), let us write

1≤n,m≤N α∈L ϕ(α + y n -y m ) = 1 |C| 1≤n,m≤N α * ∈L * e 2iπ<y n -y m ,α * > ϕ(α * ) = 1 |C| α * ∈L * 1≤n≤N e 2iπ<y n ,α * > 2 ϕ(α * ).
We can arrange that, in the left-hand side, the sum over α ∈ L reduces to a single term obtain when α = 0. We reach the following conclusion.

Lemma 3.2. For any function ϕ supported on the elementary cell C of a lattice L, that is, supp(ϕ) ⊂ C and for any finite collection (y 1 , . . . , y N ) satisfying the "localization property" y n -y m ∈ C for all n, m = 1, 2, . . . , N , the following identity holds:

1≤n,m≤N ϕ(y n -y m ) = 1 |C| α * ∈L * 1≤n≤N e 2iπ<y n ,α * > 2 ϕ(α * ). (3.8)
Periodic kernels. The interest of the following class of kernels lies in the fact that their spectral decomposition can be determined almost explicitly, in terms of exponential functions defined on the dual lattice. Namely, Lemma 3.2 allows us to pass from the continuous physical variables on C to the discrete Fourier variables on L * . The role of the function ρ introduced below is going to be played by (the restriction to the lattice of) the Fouier transform of an arbitrary kernel. Indeed, our definition below provides a way to transform a "seed" kernel K seed defined on R D to a periodic kernel K per defined on the lattice cell C. (See the proof in the appendix.) Proposition 3.3 (Periodic kernels associated with a generating function). Consider a discrete lattice L ⊂ R D generated from an elementary cell C, and let ρ : L * → (0, +∞) be a positive, integrable, and even function, that is,

ρ 1/2 ∈ 2 (L * ), ρ(-α * ) = ρ(α * ) ≥ 0 ( with α * ∈ L * ). (3.9a)
Then, the discrete Fourier transform of ρ extended to the whole of R D , that is,

K per (x, y) := e 2iπ<x-y,•> , ρ 2 (L * ) , x, y ∈ R D (3.9b)
defines an admissible kernel on R D which is periodic with period C and its associated Hilbert space is

H K per (C) = ϕ ∈ C(R D ) C-periodic / ϕ ρ -1/2 ∈ 2 (L * ) (3.9c) endowed with the norm f, g H K per (C) = f ρ -1/2 , gρ -1/2 2 (L * ) .
More generally, provided ρ s ∈ p (L * ) for some p ∈ [1, +∞) and s ≥ 0, we can also introduce the Banach space

H s,p K per (C) = ϕ ∈ C(R D ) C-periodic ϕρ -s/p ∈ p (L * ) .
(3.10)

Further generating techniques

The tensor technique. A broad class of examples on R D can be obtained by tensor decomposition from an admissible kernel in one dimension, say

K seed : R × R → R, namely K(x, y) = 1≤d≤D K seed (x d , y d ), x, y ∈ R D . (3.11a)
This applies, particularly, to a translation-invariant kernel, in which case we choose any even function χ seed : R → R and set

K(x, y) = 1≤d≤D χ seed (x d -y d ), x, y ∈ R D . (3.11b)
An example: the tensorial truncated kernel. The kernel K T (x, y) = 1 -|x -y| + has the form above and can also be expressed as a convolution, namely

K T (x, y) = χ(x -y) = (λ * λ)(x -y), λ(x) = 1 [-1/2,1/2] (x), x ∈ R. (3.12a)
Here, we have χ(ξ) = λ 2 (ξ) and λ(ξ) = sin(2πξ) 2πξ . More generally, in dimension D we consider

K(x, y) = 1≤d≤D K seed (x d , y d ) = 1≤d≤D 1 -|x d -y d | + , (3.12b)
written also as a convolution

K(x, y) = (λ * λ)(x -y) with λ(x) = 1 [0,1] D (x). Again, we have K(x, y) = χ(x -y) with χ(ξ) = λ 2 (ξ), and now λ(ξ) = 1≤d≤D sin(2πξ d ) 2πξ d .
The normalization technique. The transformations below can also serve as building blocks in order to adapt existing examples to a particular application. If K is an admissible kernel on R D we can consider

K norm (x, y) = K(x, y) K(x, x) 1/2 K(y, y) 1/2 , x, y ∈ R D . (3.13)
Clearly, we have K norm (x, x) = 1 and the required admissibility conditions are easily checked. The coefficients a n and a norm n of the corresponding decomposition (2.3) of a function in K K (R D ) and K K norm (R D ), respectively, are related by a norm n = a n K(y n , y n ) 1/2 , so that the two spaces are quite similar from the application standpoint. We thus regard (3.13) as a normalization procedure.

Taking sums and products. If K 1 , K 2 are admissible kernels in R D , then it is also easily checked that, a, b > 0 being some given constants,

K 3 (x, y) = aK 1 (x, y) + bK 2 (x, y), K 4 (x, y) = K 1 (x, y)K 2 (x, y), (3.14)
With the notation used in Section 2 the matrices K 3 (Y, Y ) and K 4 (Y, Y ) are symmetric positive definite, as follows easily from a standard linear algebra argument. Zonal kernels. If a function φ is such that (x, x ) → φ(xx ) is a one-dimensional kernel, then the following formula

K zonal (x, y) = φ(< x, y >), x, y ∈ R D (3.15)
(where < x, y > stands for the Euclidian inner product) defines an admissible kernel. This class is often used by the artificial intelligence community. Convolution kernels. Another class of translation-invariant kernels can be generated by choosing a function

λ ∈ L 2 (R D ), R D λ(x)λ(-x) dx = 0, R D ( λ) 2 dξ = 1 (3.16a)
and then defining our generating function χ by the convolution formula

K(x, y) := χ(x -y) = (λ λ)(x -y), (λ λ)(y) = R D λ(x)λ(y -x) dx R D λ(x)λ(-x) dx , y ∈ R D . (3.16b)
Indeed, it is clear that χ(0) = 1 and χ = ( λ) 2 ≥ 0, so that the Fourier transform of χ is a probability measure. This class of kernels is used in machine learning, for instance in combination of multi-layer neural networks.

4 Designing kernels on a bounded domain

Standard kernels taken as seed data

We focus our attention to the standard examples listed in Table 1.1, that are radially-symmetric, translationinvariant kernels

K(x, y) = χ(x -y). Exponential kernel K E . The choice χ E (ξ) = (1 + |ξ| 2 ) -m (with m > D/2) leads to the standard Sobolev space H K E (R D ) = W m,2 (R D )
, and is a standard choice in the numerical analysis literature.

Multiquadric kernel K M . The choice χ M (ξ) = e -|ξ| is only Lipschitz continuous at the origin and is relevant for representing sufficiently smooth functions with polynomial decay, hence provides (slightly) more information than the Gaussian one (below).

Gaussian kernel K G . The choice χ G (ξ) = e -|ξ| 2 provides an exponential decay in both the Fourier and the physical spaces. Both functions χ G ≥ 0 and χ G ≥ 0 are globally positive. The Gaussian kernel is adapted to the description of smooth and fast decaying functions which have"almost" compact support in physical and Fourier variables. Hence, the function space H K G (R D ) is "small" and, in term, provides limited "information" on the functions.

Truncated kernel K T . A more interesting and also quite standard choice is obtained by truncation in the physical space, namely, χ T (x) = (1 -|x| 1 ) l + (with l ≥ D/2) where the notation a + := sup(a, 0) stands for the positive part. This kernel is only Lipschitz continuous in the physical space.

Periodic kernels of interest

Objective. We now present, in their rescaled form, the four periodic kernels already listed in Table 1.2 (and associated with each of the examples in Table 1.1). Periodic kernels are translation-invariant, i.e. K per (x, y) = χ per (x -y) and we plot the corresponding functions χ per in Figure 4.1. Moreover, we also plot the level set of the Fourier transform of χ in Figure 4.2 for the two-dimensional case. For the sake of simplicity in the notation, since the lattice and the dual lattice coincide we simply write α (instead of α * ) for a general element of the lattice or dual lattice.

Periodic tensorial exponential kernel K per E . Consider the one-dimensional exponential kernel given by

K E (x, y) = exp(-|x -y|) = χ E (|x -y|), x, y ∈ R, (4.1a)
with χ E (ξ) = 2 1+4π 2 ξ 2 . We make this kernel tensorial and periodic using the localization method in Section 3.3 (see also (6.8) Here, we have also introduced a parameter 0 < τ D → 0 as D → +∞. Observe that the Fourier transform of τ D e -τ D |x| is

below), that is, for x, y ∈ [0, 1] D K per E (x, y) = χ per E (x -y) = α∈Z D e 2iπ<x-y,α> D d=1 (1 + 4π 2 α 2 d /τ 2 D ) = D d=1 α d ∈Z e 2iπ(x d -y d )α d 1 + 4π 2 α 2 d /τ 2 D . ( 4
2 1+4π 2 α 2 d /τ 2 D
. Thus, thanks to the Poisson formula (3.6), our kernel coincides with

K per E (x, y) = α∈Z D τ D exp(-τ D |x -y + α| 1 ) = τ D α∈Z D D d=1 exp(-τ D |x d -y d + α d |) = D d=1 2τ D e τ D -1 exp(τ D |x d -y d |) + exp τ D (1 -|x d -y d |) . (4.1c)
We plot the corresponding function χ per E in Figure 4.3 for several dimensions. The Banach space associated to this kernel (defined in (3.10)) reads

H s,p K per E ([0, 1] D ) = ϕ is periodic on [0, 1] D D d=1 (1 + 4π 2 α 2 d /τ 2 D ) s/p ϕ p (Z D ) < +∞ . (4.1d)
In particular, the space H 1,1

K per E ([0, 1] D ) has been found to be relevant in mathematical finance. Periodic Gaussian kernel. The Gaussian kernel is translation-invariant, radial, and smooth, and reads

K G (x, y) = exp(-|x -y| 2 ) = D d=1 exp(-|x d -y d | 2 ) = χ G (x -y), (4.2)
with χ G (ξ) = 2 -D/2 exp(-|ξ| 2 /4). We define its periodic version as

K per G (x, y) := 1 τ D α∈Z D D d=1 exp(- |x d -y d + α| 2 4τ 2 ) = D d=1 ϑ 3 (iπ|x d -y d |, τ D ), (4.3)
where

ϑ 3 (z; τ D ) = n∈Z exp(-π 2 τ ) n 2
cos(2πnz) is nothing but the third Jacobi-theta function. Here, τ D is a numerical parameter. The diagonal term enjoys the following decay

K per G (x, x) = ϑ 3 (0; τ D ) D (1 + 2e -π 2 τ D ) D , and we thus choose τ D = ln(2D) π 2 to ensure K per G (x, x) ≤ e.
Let us denote (with some abuse of notation)

K per G (x, y) = χ per G (x -y) = d χ per G (x d -y d ).
We plot the function χ per G in Figure 4.4 for several dimensions.

Periodic multiquadric kernel. The one-dimensional multiquadric kernel is translation-invariant, radially-symmetric, and smooth, and reads

K M (x, y) = 1 1 + |x -y| 2 = χ M (|x -y|), x, y ∈ R, (4.4a)
with χ M (ξ) = exp(-|ξ|), so this kernel is nothing but the Fourier transform of the exponential kernel. above. We define its periodic version by setting, for all x, y ∈ [0, 1] D ,

K per M (x, y) := α∈Z D D d=1 e 2iπ<x-y,α>-τ D |α|1 = D d=1 α d ∈Z e 2iπ(x d -y d )α d -τ D |α d | , (4.4b)
K per M (x, y) = D d=1 α d ∈Z 1 1 + (x d -y d + α d ) 2 /τ 2 D = D d=1 sinh(2πτ D ) cosh(2πτ D ) -cos(2π(x d -y d )) . (4.4c)
The diagonal contribution enjoys the following asymptotics

K per M (x, x) = D d=1 α d ∈Z 1 1 + α 2 d /τ 2 D = sinh(2πτ D ) cosh(2πτ D ) -1 D = coth(πτ D /2) D . Hence we choose τ D = 2 coth -1 (1+1/D) π = 2 π (ln(2 + 1/D) -ln(1/(D + 1)
)) to ensure that K per M (x, x) ≤ e uniformly for any dimension. Finally, we write

K per M (x, y) = χ per M (x -y) = d χ per M (x d -y d ) (
with some abuse of notation), and we plot the function χ per M in Figure 4.5 for several dimensions. Periodic truncated kernel. The truncated kernel is translation-invariant and Lipschitz continuous only, and reads

K T (x, y) = sup(1 -|x -y|, 0) = χ(x -y), ϕ(ξ) = sin 2 (πξ/2) (πξ/2) 2 . (4.5a)
We emphasize that its tensor product

D d=1 sup(1 -2|x d -y d |, 0) = D d=1 χ(x d -y d ), ϕ(ξ) = D d=1 sin 2 (πξ d /2) (πξ d /2) 2 (4.5b)
K per T (x, y) := α∈Z D D d=1 τ D sup(1 -τ D |x d -y d + α d |, 0) = α * ∈Z D D d=1 τ D sin 2 (πα d /τ D ) (πα * d /τ D ) 2 e 2iπ(x d -y d )α * d , (4.5c)
where 1 ≤ τ D ≤ 2 is a parameter, chosen so that the sum is finite, thus

K per T (x, y) = D d=1 α d ∈{-1,0,1} τ D sup(1 -τ D |x d -y d + α d |, 0). (4.5d)
The diagonal term is K per T (x, x) = τ D D , and we thus choose τ D = 1+1/D in order to ensure that K per T (x, x) ≤ e. Writing K per T (x, y) = χ per T (x -y) = d χ per T (x d -y d ), we plot χ per T in Figure 4.6 for several dimensions.

Transported kernels of interest

Objective. An admissible kernel K defined on an open set Ω is said to be a compactly supported if it extends to a continuous function on the closure Ω and this extension vanishes on the boundary ∂Ω. For our numerical experiments, we design four compactly supported kernels defined on the unit cube [0, 1] D , determined from the four examples in Table 1.1. Figure 4.7 displays K(x, 1/2) for the two-dimensional case.

Transported tensorial exponential kernel. In view of the standard expression of the exponential kernel, we introduce the following transported version defined on the unit cube

[0, 1] D K tran E (x, y) = 1 β D exp -τ D erf -1 (2x -1) -erf -1 (2y -1) 1 , (4.6a)
where erf -1 (x) = erf -1 (x 1 ), . . . , erf -1 (x D ) and |x| 1 = d |x d | and we the erf function reads erf(x) =

(2/ √ π)

x 0 e -y 2 dy. By Proposition 3.1, this kernel corresponds to the transport of the localized kernel

K loc E (x, y) = 1 β D exp -τ D |x -y| 1 exp -|x| 2 -|y| 2 , (4.6b)
obtained from the map S(x) = erf -1 (2x -1). Here, we choose Transported multiquadric kernel. In view of the expression of multiquadric kernel, we consider the following localized version on the unit cube

τ D = √ π D , β D = e τ 2 D /4 (1 -erf(τ D /2) D with, independently of the dimension D, K tran E (x, x) = 1 β D (1 + τ D √ π ) D ≤ e, [0,1] D ×[0,1] D K tran E (x, y) dxdy = 1.
[0, 1] D K tran M (x, y) = β D (1 + τ 2 D |erf -1 (2x -1) -erf -1 (2y -1))| 2 ) (D+1)/2 , (4.7a)
corresponding to the following transported kernel

K loc M (x, y) = β D exp(-|x| 2 -|y| 2 ) (1 + τ 2 D |x -y| 2 ) D , x, y ∈ R D , (4.7b)
where we choose

τ D = 2 D , β D = 1 √ π exp(1/τ 2 D )(1 -erf(1/τ D ))(1/τ D ) D 1 + τ 2 D 2 D ,
as τ D → 0, determined so that

K tran M (x, x) = β D ≤ e, [0,1] D [0,1] D K tran (x, y)dxdx β D √ π exp(τ 2 D /2)(1 -erf(1/τ D )(1/τ D ) D = 1.
Transported Gaussian kernel. In view the expression of the Gaussian kernel, we introduce the following localized version on the unit cube

[0, 1] D K G (x, y) = φ(x -y) = β D exp -τ 2 D D d=1 erf -1 (2x d -1) -erf -1 (2y d -1) 2 , x, y ∈ [0, 1] D , (4.8a)
where we choose

τ D = 2 D and β D = (1 + τ 2 D ) D/2 , determined so that K tran G (x, x) = β D ≤ e, [0,1] D [0,1] D K tran G (x, y)dxdx β D (1 + τ 2 D ) D/2 = 1.
Transported truncated kernel. In view of the expression of the truncated kernel, we consider the following localized version on the unit cube [0, 1]

D K tran T (x, y) = β D (1 + τ D |erf -1 (2x d -1) -erf -1 (2y d -1))| 2 ) D , x, y ∈ [0, 1] D , (4.9a)
corresponding to the following transported kernel

K loc T (x, y) = β D exp(-|x| 2 -|y| 2 ) (1 + τ D |x -y|) D , x, y ∈ R D . (4.9b)

Further constructions

A compactly supported, tensor product kernel. The example given now is not translation-invariant and is supported on the unit interval Ω = [0, 1]:

K seed (x, y) := y(1 -x)1 y≤x + x(1 -y)1 x≤y , x, y, ∈ R. (4.10a)
Observe that K seed (0, y) = K seed (1, y) = 0, and that the Fourier technique above does not apply. For future reference we compute

∂ y K seed (x, y) = (1 -x)1 y<x -x1 y>x , ∂ 2 y K seed (x, y) = -δ x [0,1] K seed (x, y)dy = 1 2 x(1 -x), [0,1] 2 K seed (x, y)dxdy = 1 12 , (4.10b)
together with the following integration by part formula:

[0,1] ∂ 2 K seed (•, x)∂ 2 K seed (•, y)dz = [0,1] ∂ 2 (K seed (•, x)∂ 2 K seed (•, y))dz + K seed (x, y) = K seed (x, y).
(4.10c) Hence, the space H K seed (0, 1) is the homogeneous Sobolev space Ḣ1 0 ([0, 1]) of functions that vanish at the boundary, defined by density from the set of smooth functions supported on (0, 1) in the 2 (0, 1)-norm. The eigenvalues and eigenfunctions of the spectral decomposition is given by solving

[0,1] K seed (•, x)ζ i (x)dx = λ i ζ i or, equivalently, ∂ 2 x ζ i = -1 λi ζ i . We find ζ i (x) = sin(iπx), λ i = 1 (iπ) 2 for i = 1, 2, . . . A multi-dimensional version.
Using the kernel K seed in (4.10a), in general dimensions we consider Ω = [0, 1] D and the following kernel

K(x, y) = 1≤d≤D K seed (x d , y d ), = 1≤d≤D y d (1 -x d )1 y d ≤x d + x d (1 -y d )1 x d ≤y d . (4.11a)
In view of

Ω (∂ 1 . . . ∂ D )K(•, x)(∂ 1 . . . ∂ D )K(•, y)dz = 1≤d≤D [0,1] ∂ d K seed (x d , •)∂ d K seed (y d , •)dz = K(x, y), (4.11b)
the space generated by this kernel corresponds to the norm (0,1)

D |ϕ(x)| 2 + |∂ 1 . . . ∂ D ϕ(x)| 2 dx.
Similarly to the one-dimensional case, we compute

∂ d K(x, y) = (∂K seed )(x d , y d ) e =d K seed (x e , y e ), 1≤d≤D ∂ 2 d K(x, y) = (-1) D δ x [0,1] D K(x, y)dy = 1≤d≤D 1 2 x d (1 -x d ), [0,1] D ×[0,1] D K(x, y)dxdy = 1 12 D . (4.11c)
The eigenvalues and eigenfunctions of the the spectral decomposition are now

ζ α (x) = 1≤d≤D sin(α d πx d ), λ α = 1 1≤d≤D (α d π) 2 , α = (α 1 , . . . , α D ) ∈ N D . (4.11d)
5 Three formulation of the discrepancy function

Formulation in physical variables

We now factor out the integration error into (1) a factor depending upon the regularity of the function ϕ under consideration, measured in the H K -norm, which is is independent of the choice of the interpolation points, and (2) a factor depending solely upon the kernel K and the mesh points, which is independent of the choice of the function. Three equivalent formulations of the second term are now derived in the physical, spectral, or discrete Fourier variables. Although these formulations are in principle equivalent, they shed a very different light on the problem of interest. We begin with an expression in the physical variables, which we will not use directly for the derivation of actual estimates, as it appears to be difficult to work with.

Ω ϕ(x)dx - 1 N 1≤n≤N ϕ(y n ) ≤ E K (Y, N, D) ϕ H K (Ω) (5.1)
for any set Y = (y 1 , . . . , y N ) of points in Ω, with the error function

E K (Y, N, D) 2 := 1 N 2 N n,m=1 Ω×Ω K(x, y) + K(y n , y m ) -K(y n , y) -K(x, y m ) dxdy (5.2a)
or, equivalently,

E K (Y, N, D) = e K (Y ) H K (Ω) , e K (Y ) := Ω K(•, x) - 1 N 1≤n≤N K(•, y n ) dx.
(5.2b) Remark 5.2. We observe that the integrand in (5.2a) is non-negative and can be also written in terms of the pseudo-distance D (see (2.2)), that is,

E K (Y, N, D) 2 = - 1 2N 2 N n,m=1 Ω×Ω D(x, y) + D(y n , y m ) -D(y n , y) -D(x, y m ) dxdy. (5.3)
Recall that (2.12) shows that D(x, y) measures the "distance" between K(•, x) and K(•, y), and (5.3) provides an equivalent representation of the error.

Formulation in spectral variables

E s,p K (N, D) := inf y 1 ,...,y N ∈Ω i≥1 λ sp /2 i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) p (5.4) (for s ≥ 0 and p begin the dual exponent to p ∈ [1, = ∞)) satisfies E s,p K (N, D) ≤ i>N λ sp /2 i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) p , (5.5a)
where y 1 , . . . , y N solve the system of equations

Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) = 0. (5.5b)
Whenever the eigenvalues and functions are known explicitly, the spectral decomposition associated with the kernel can be used in combination with the following formula for the error function. Observe that we can now deal with the spaces H s,p K (Ω) with general integrability and differentiability exponents, and

recall that E K (Y, N, D) = E 1,2 K (Y, N, D).
Proposition 5.4 (Factorization in spectral variables). Consider an admissible kernel K = K(x, y) defined on an open set Ω ⊂ R D . Then, for any function ϕ in the corresponding Banach space H s,p K (Ω) associated with exponents s > 0 and p ∈ [1, +∞). Then, for any set Y = (y 1 , . . . , y N ) of points in Ω and any function ϕ ∈ H s,p K (Ω), one has

Ω ϕ(x)dx - 1 N 1≤n≤N ϕ(y n ) ≤ E s,p K (Y, N, D) ϕ H s,p K (Ω) , (5.6)
where the error function (with 1/p + 1/p = 1)

E s,p K (Y, N, D) := i≥1 λ sp /2 i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) p 1/p = e K (Y ) H s,p K (Ω) (5.7)
is expressed in terms of the eigenfunctions and eigenvectors ζ i , λ i of the Mercer representation of the kernel K.

Formulation in discrete Fourier variables

Assuming now more structure on the kernel and, specifically, assuming that it is based on a discrete lattice, we express the error function in a third form based on the Poisson formula.

Proposition 5.5 (Factorization in discrete Fourier variables). Consider a kernel based on a discrete lattice L with elementary cell C, say K per (x, y) := 1 |C| e 2iπ<x-y,•> , ρ 2 (L * ) , determined from a generating function ρ defined on the dual lattice L * . Then, for any set Y = (y 1 , . . . , y N ) of points in C and any L-periodic function ϕ ∈ H s,p K per (Ω) (for some s > 0 and p ∈ [1, +∞) with 1/p + 1/p = 1), one has (5.6) with

E s,p K per (Y, N, D) = α * ∈L * ρ s C e 2iπ<x,α * > dx - 1 N n e 2iπ<y n ,α * > p 1/p . (5.8)
6 Controlling the discrepancy function

Physical variables

We seek for a set of points minimizing the error function and, first in the physical variables, we obtain the following result whose proof is given in the appendix.

E K (N, D) := inf Y =(y 1 ,...,y N ) E K (Y, N, D) > 0.
Moreover, if Y be a minimizer, then the gradient of the functional vanishes at Y , namely

Ω ∇K(x, y n )dx = 1 N N m=1
(∇K)(y n , y m ), n = 1, 2, . . . , N. (6.1)

Furthermore, for all functions H Y K (Ω) (as defined in (2.3)), the following integration formula holds (without error term):

Ω ϕ(x)dx = 1 N 1≤n≤N ϕ(y n ), ϕ ∈ H Y K (Ω). (6.2)
Here, ∇K(x, y) stands here for any of ∇ x K(x, y) or ∇ y K(x, y), since the kernel K is symmetric in its two arguments. The following comments are in order:

• The functional is be totally symmetric with respect to its arguments so, clearly, a minimizer is never unique. Yet, suppose that our kernel K(x, y) is concave, then the following semi-discrete algorithm

d dt y n = Ω ∇K(x, y n )dx - 1 N N m,n=1 (∇K)(y n , y m ), n = 1, 2, . . . , N (6.3)
converges toward a minimum of the functional E K . Observe also that the dynamical system under consideration involves two terms: ∇K(y n , y m ) tends to push points away from each other, while the integral term tends to attract the points toward the mass-center of Ω. These two competitive effects leads to a non-trivial distribution of the points.

• With the partition of unity defined in (2.10b), any minimizer must satisfy Ω δ n (x)dx = 1 N , for n = 1, 2, . . . , N.

Best discrepancy sequences in one dimension. We are now in a position to analyze an example in dimension D = 1, where we have fully explicit expressions. Proposition 6.2 (The case of dimension D = 1). For the kernel (4.10), the equally-spaced sequence

y n = 2n -1 2N , n = 1, 2, . . . , N (6.4)
is the unique solution (up to re-ordering) that minimizes the error function Y → E K (Y ) introduced in (5.1) and, with this sequence, the following estimate holds:

[0,1] ϕ(x)dx - 1 N 1≤n≤N ϕ(y n ) ≤ 1 √ 6N ϕ H 1,2 K ([0,1]) , ϕ ∈ H 1,2 K ([0, 1]).

Spectral variables

For clarity in the presentation, let us repeat here our previous observation. Since the following results will be of a direct application in the numerical section, we include the remaining proofs in the main text.

one has

E s,p K (N, D) ≤ i>N λ sp /2 i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) p 1/p , (6.6a)
where (y 1 , . . . , y N ) ∈ R N ×D is given by the following system of equations

Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) = 0. (6.6b)
Proof. The expression K N (x, y) = 1≤n≤N λ n ζ n (x)ζ n (y) defines an admissible kernel, for which we have an orthogonal decomposition of the form

H K (Ω) = H K N (Ω)⊕H ⊥ K N (Ω).
It suffices to consider the minimization problem

E K N (N, D) = inf y 1 ,...,y N ∈Ω 1≤i≤N λ i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) 2 1/2 and 1≤i≤N λ i Ω ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) (∇ζ i )(y n ) = 0, n = 1, . . . , N.
Best discrepancy sequences in one dimension. We can now revisit the example in dimension D = 1 treated in Proposition 6.2. Proposition 6.4 (The case of dimension D = 1.). For the kernel (4.10) and the equally-spaced sequence

y n = 2n -1 2N , n = 1, 2, . . . , N, (6.7)
the error function Y → E s,p K (Y, N, 1) introduced in (5.7) can be bounded by, for s > 1 -1/p and some constant C s,p > 0:

[0,1] ϕ(x)dx - 1 N 1≤n≤N ϕ(y n ) ≤ C s,p N s ϕ H s,p K ([0,1]) , ϕ ∈ H s,p K ([0, 1]).
Proof. The general case follows using the spectral formulation (5.7), where ζ i (x) = sin(iπx) and λ i = 1 (iπ) 2 . We introduce

c i := [0,1] ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) = [0,1] sin(iπx)dx - 1 N 1≤n≤N sin(iπ 2n -1 2N ) = 1 -cos(iπ) iπ - sin 2 ( iπ 2 ) N sin( iπ 2N )
.

Hence, c 2i = 0, c 2i+1 = 2 (2i+1)π - 1 N sin( (2i+1)π 2N )
and, in view of (5.7),

E s,p K (Y, N, 1) p = i≥0 1 ((2i + 1)π) sp 2 (2i + 1)π - 1 N sin( (2i+1)π 2N ) p = i≥0 λ sp /2 2i+1 2λ 1/2 2i+1 - 1 N sin( 1 2N λ 1/2 2i+1 ) p where 1/p + 1/(p ) = 1. In particular, E 1,2 K (Y, N, 1) = 1 √ 6N
, as computed in Proposition 6.2. For general s, p, we get the expansion 2λ

1/2 2i+1 - 1 N sin( 1 2N λ 1/2 2i+1 ) λ -1/2 2i+1 N 2 as λ -1/2 2i+1
2N < 1, that is iπ < N , and is bounded elsewhere.)

The error term in (5.7) is bounded as follows:

E s,p K (Y, N, 1) p ≤ πi≤N C N 2p λ (s-1)p /2 2i-1 + C iπ≥N λ sp /2 2i+1 ,
where the first term in the right hand-side is uniformly bounded by C N 2p provided (s -1)p > 1, that is, s > 1/p. The second term is bounded by the first term C N sp provided sp > 1. Thus we can control E s,p K (Y, N, 1) by C/N s in the range sp > 1 at least. Observe finally that we are led to the same estimate, since

(E s,p K (N, 1)) p ≤ i>N λ sp /2 i [0,1] D ζ i (x)dx - 1 N 1≤n≤N ζ i (y n ) p 1 N sp .

Discrete Fourier variables: toward a sharp estimate

Our third formulation provides us with the most practical setup for the applications, when the spectral decomposition may not be explicitly available. is an admissible kernel with period C and in the Hilbert space H K per (C) one has

Ω ϕ(x)dx - 1 N N n=1 ϕ(y n ) ≤ E K per (Y, N, D) ϕ H K , (6.9)
where

E K per (Y, N, D) 2 = 1 N 2 N n,m=1 K per (y n , y m ) -1, (6.10)
where Y is any sequence of N points in C. Moreover, setting

Y := arg inf Y ∈C N E K (Y, N, D), one has E K (N, D) 2 ≤ n>N ρ(α * n ) N 2 N n=1 e 2iπ<y n ,α * n > 2 , (6.11)
where the ordering n → α n is defined so that n → ρ(α * n ) is a decreasing sequence along the points α ∈ L * while y 1 , . . . , y N are defined by solving the set of equations (when

α * n = 0) N m=1 e 2iπ<y m ,α * n > = 0, 1 ≤ n ≤ N. (6.12)
We conjecture that the system (6.12) does admit a solution for any set {α * n ∈ L * } 1≤n≤N . Numerically, for this system for each of the periodic kernels of interest we computed (see below) a numerical approximation for a broad range of dimensions D and integers N . Moreover, the existence of such points can be established rigorously for several examples of kernels. Importantly, we tested numerically (see below) that the following estimate is sharp:

E K (N, D) 1 N n>N ρ(α * n ). (6.13)
In our examples, we will be able to compute the above sum, both numerically and analytically, and it will be proven to provide us with a sharp estimate of the discrepancy error, altough we cannot establish rigourously the validity of this estimate. In the context of Proposition 5.5, we will compute numerically

n>N ρ(α * n ) = n≥0 ρ(α * n ) - n≤N ρ(α * n ) = K per (0, 0) - n≤N ρ(α * n ). (6.14)
Another theoretical standpoint is obtained by the following "level-set argument", namely by relying on the approximation

n>N ρ(α * n ) ξ/ ρ(ξ)<m -1 (N ) ρ(ξ)dξ, (6.15)
where m -1 is the inverse of the function m( ) = |{ρ(ξ) > }|, where ρ(ξ) is a smooth extension of ρ on the whole space R D satisfying R D ρ(ξ) dξ = 1 with ρ(ξ) ≥ 0 and ρ(0) = 1. Unfortunately, this approximation is very inaccurate in the examples we have considered.

In particular, if a transport map of the function is known, that is a one-to-one map y : R D → Ω such that dy = ρ(ξ)dξ, then, denoting its inverse as ξ(y), the last integral reduces to measuring the following set

n>N ρ(α * n ) y:ρ(ξ(y))<m -1 (N ) dy = {y : ρ(ξ(y)) < m -1 (N )} . (6.16)
Let us make some further remarks:

• If one want to study a general kernel K(x, y) = χ(x -y) with x, y ∈ Ω, then one can compute an upper-bound using 6.11 with the Fourier coefficients on the doubled lattice L 2 , generated by the cell Ω 2 := {x -y / x, y ∈ Ω}. This is done using the Poisson formula, i.e.

K(x, y) = χ(x -y) = e 2iπ<x-y,•> , ρ l 2 (L * 2 ) |Ω 2 | , ρ(ξ) = 1 (2π) D/2 Ω2 χ(ζ)e -2iπ<ξ,ζ> dζ.
(6.17)

However, in practical terms, computing Fourier coefficients can be a quite difficult task.

• We can specify directly the Fourier coefficients χ(α * ) in the light of Proposition 5.5. This obviously simplifies the problem of finding suitable parameters since they are automatically set by the two conditions χ(0) = 1 and χ ∈ 1 (L * ).

• It is clear from this result that the price to pay, when the dimension of the problem increases, is increasing also the regularity of the kernel, and the space H K contains functions that are more regular as the dimension increases.

• The lattice L certainly plays an important role in finding the best constant E K (N, D). If χ is symmetric, the present discussion connects with the sphere packing problem. For instance, for D = 2, the best lattice is the hexagonal one, for D = 3 the dodecahedral one (Kepler conjecture), while D = 8, 24 was treated in [START_REF] Viazovska | The sphere packing problem in dimension 8[END_REF]. Other cases remain opened.

Proof of Proposition 5.5. To derive (6.9), we recall from (5.2a) that

E K (Y, N, D) 2 = C×C K(x, y)dxdy + 1 N 2 N n,m=1 K(y n , y m ) - 2 N n Ω K(x, y n )dx.
We use the fact that K(x, y) = χ(x -y) is L-periodic and deduce that the right-hand side is a constant, and does not depend on y n , therefore

E K (Y, N, D) 2 = α * ∈L * χ(α * ) C e 2iπ<x,α * > dx - 1 N N n=1 e 2iπ<y n ,α * > 2 .
We have C e 2iπ<x,α * n > dx = 1 if α * n = 0, while it is 0 otherwise.

When α * = 0 we consider

E K (Y, N, D) 2 = α * ∈L * ,α * =0 χ(α * ) N 2 N n=1 e 2iπ<y n ,α * > 2 . (6.18)
Using the ordering n → α n associated with the function χ(α * n we arrive at the system of equations (6.12).

The existence of solutions to (6.12) can for instance be established for the canonical lattice C = [0, 1] D and L = L * = Z D , and any set of points {α * :

L d ≤ α * d < R d }, with size N = D d=1 (R d -L d ).
Without loss of generality we assume that L d = 0 and we set N e = e-1 d=1 R d for e = 1, . . . , D (with the convention that N 1 = 1). We define a one-to-one map from the set {1 ≤ n ≤ N } to {0 ≤ n d < R d } d=1,...,D by using the map in

n(n 1 , . . . , n D ) = d≤D N d n d , α * n = n 1 , . . . , n D , y n = n d (n) R d d=1,...,D .
We compute .

N n=1 e 2iπ<y n ,α * m > = N n=1 e 2iπ d n d (n) R d n d (m) = R1-1 n1=L1 e 2iπ n 1 R 1 n m 1 . . .
It is regular and admits the derivatives

∂ y k I(Y ) = N n=1 4iπα * n N m=1 e 2iπ<y m -y k ,α * n >
. Therefore, we can use a gradient-based method in order to numerically compute the solutions, and we expect that any local minimum of the functional I(Y ) will also be a global minimum.

A comparative study for a selection of kernels 7.1 Our strategy for the numerical study

The aim of this section is to numerically compute the optimal sequences Y and the discrepancy functions sdefined by (with Ω = [0, 1] D )

Y = arg inf Y ∈Ω N E K (Y, N, D), E K (N, D) = E K (Y , N, D). (7.1)
We treat two classes of interest:

• Transported kernel: we recall (see (5.2a)) that the following quantity provides the worst Monte-Carlotype integration error for any admissible kernel K = K loc localized to the set Ω:

E K (Y, N, D) 2 = 1 N 2 N n,m=1 Ω×Ω
K loc (x, y)+K loc (y n , y m )-K loc (y n , y)-K loc (x, y m ) dxdy. (7.2)

• Periodic kernels: in the context of Proposition 5.5, the expression of the discrepancy function takes a simpler form and in the limit N, D → ∞ we can compare our results to the asymptotic expression

E K (Y, N, D) 2 = 1 N 2 N n,m=1 K per (y n , y m ) -1, E K (N, D) ≤ n>N ρ(α * n ). (7.3)
We recall that α n are chosen on the dual lattice so that n → χ(α * n ) is decreasing. We choose here the canonical lattice L = Z D for simplicity.

We start from the four examples in Table 1.1, that are defined in the whole space R D , and design the eight kernels listed in Table 1.2. A normalization parameter, depending on the dimension was determined explicitly for each case. In each case, we compute E K (Y, N, D) for the following two sequences:

• Y is a randomly chosen sequence.

• Y is a numerical sequence that approximates the optimal solution Y .

The numerical results for E K (Y, N, D) are displayed as tables for N = 16, 32, . . . , 512 and D = 1, 2, 4, . . . , 128. and we can compare the error obtained for various values N, D and various kernels. The kernels in the second line of Table 1.2 are periodic, so that we can compute the upper bounds given by the formula (6.11) and we can compare it to the numerically computed results. As a general rule, the computation time grows like C DN 3 for some rather large constant C > 0. For instance, on a machine Intel I7 7700K, each result in the tables presented below took between 10 secondes and up to two minutes for the tests in dimension D = 128 and with N = 512. In order to ease the comparison between the various tests and kernels we have considered, we are collecting our results in Tables 7.1, 7.2, and 7.3 when N = 512 and in dimensions 1, 16, 128. Remarkably, our asymptotic analysis provides an accurate prediction of the numerical error (except for the truncated kernel in the largest dimension D = 512).

Periodic kernels

Strategy for periodic kernels. For each of the four periodic kernels in Table 1.2, we present three tables:

(A) E K (Y, N, D) for a randomly chosen sequence Y .

(B) E K (Y, N, D) for a numerical sequence that approximates the optimal solution Y .

(C)

1 N n>N ρ(α * n ) which is our theoretical (but yet heuristic) asymptotic rate of convergence (6.13).

From a computational point of view, some remarks are in order:

• In item (A) above, the random sequences Y are computed using the standard Mersenne-Twister mt19937. The error bound E K (Y, N, D) can be computed using the left-hand formula in (7.3).

• In item (B) above, in order to compute the optimal Y we have two strategies:

-We can either minimizing directly the left-hand side of (7.3). This optimization problem can be solved easily by a gradient descent and is computationally tractable. However, depending upon the choice of the kernel this method can show poor performances (i presence of almost vanishing gradients.

-Or else, we can compute the first values α * n for n = 0, . . . , N -1 and then solve the equation (6.12). We recall here that N m=1 e 2iπ<y m ,α * n > = 0 for 1 ≤ n ≤ N . We used this method and it always gave ni principle good numerical results. However, to solve the relevant system we must use an algorithm whose complexity is of order N 3 , that leads to extremely long execution times as N becomes large.

-Once the sequence Y is computed, the error E K (Y, N, D) is obtained using the left-hand side of (7.3).

• Finally, in item (C) above, we compute first the set of values χ(α * n ) for n = 1, . . . , N by using a graph-search type algorithm. Observe that this step can be challenging in the case of a rather oscillatory function χ such as the truncated kernel.

For each periodic kernel we provide below the asymptotic formula derived from our level-set arguments, although it is not very accurate within the range of N, D under consideration. We plot in this section the distribution of sequences of points approximating the best discrepancy sequences Y . These figures corresponds to N = 256 and D = 2. Observe that none of these distributions is radially-symmetric, as the corresponding periodic kernels are not.

Numerical results for the periodic tensorial exponential kernel. For the periodic kernel (4.1a) the level-set method (6.15) applies and lead us to the bound E K (N, D) log(N ) D-1 N , which is probably far from being optimal but coincides with the Koksma-Hlawka inequality for low-discrepancy sequences. Here, the space H 1,1 K ([0, 1] D ) (see (4.1d)) is similar to the space of functions with bounded variation BV ([0, 1] D ). The three relevant tables are given here. Observe that the discrepancy error for the numerical points is always much smaller than the discrepancy error for the random points, as was expected.

Note that the asymptotic convergence rate is quite close to the one obtained with the numerical sequence in Table 7.5. There exist cases for which this convergence rate is smaller than the theoretical bound. This might seem to be a contradiction, since the latter bound is supposed to be the minimum over all possible sequences. However, several sources of (yet small) numerical error arise due to the search algorithm of the decreasing sequence α * n . In addition, a second (small) error arises from the estimate 6.14, which is probably sharp but yet only an approximation of the theoretical bound. Numerical results for the periodic multiquadric kernel. For the periodic multiquadric kernel and using the level-set argument in (6.15) we can derive the bound E

K (N, D) 2 D τ D D 1 -exp -N 1/D D 2τ D
. The three tables are now presented and the same observations as above for the exponential kernel can be made. Note again the existence of a small numerical error for the two-dimensional case, as the error discrepancy for the optimized sequence is not decreasing, namely for D = 2 and N = 256, 512. Observe that the error vanishes for the one-dimensional case provides N ≥ 64. This was expected since we are dealing with a kernel generating a space consisting of periodic function, whose Fourier coefficient decreases at exponential rate, as predicted by our theoretical convergence rate above. From a numerical point of view, it is expected to be a finite dimensional space, having few basis functions in low dimensions. Numerical results for the periodic Gaussian kernel. For the periodic kernel (4.3) we can use our level-set arguments in (6.15) and we c an arribve at the bound E K (N, D) ≤ N 1-1/D exp(-N 2/D ). We now present the three tables of interest and again the results are quite similar to the exponential and multiquadric cases. Numerical results for the periodic truncated kernel. For the periodic kernel (4.3) we can use our level-set method arguments and we arrive at the bound E K (N, D) ≤ log(N ) D-1 N 2

. We now present the three relevant tables as above and the observations we made concerning these results are quite similar.

Transported kernels

The numerical results we performed on transported kernel kerne are very similar to the one we just described for periodic kernels. We content with pointing out here significant differences between periodic and transported kernels:

• We did not numerically compute the theoretical convergence rates for the four transported kernels. Indeed, we would need to compute the Fourier coefficients appearing in Proposition 5.5. This appears to be too costly and computationally out of reach with the existing techniques.

• On the other hand, we did compute the error function thanks to its expression (7.2), evaluated via a direct Monte-Carlo method. This is still computationally expensive, and adds some additional white noise to the final results.

Thus, in this section, we restrict ourselves with presenting the plot in Figure 7.2 displaying the distribution of points that approximates the best discrepancy sequences Y . This corresponds to the choice N = 256 and the dimension D = 2. Some observations are in order:

• Among our four transportation-based kernels, three are radially-symmetric, that is, Gaussian, multiquadric, and truncated. Their best discrepancy sequences are expected to be radially-symmetric. Indeed, our numerical tests confirm this property for the multiquadric and truncated kernels, which enjoy similar properties.

• However the Gaussian kernel did not lead us to a radially-symmetric result. This is due to a numerical challenge we can refer to as the "vanishing gradient problem" (a terminology used in the artificial intelligence community), due to the fact that the functional we want to minimize is almost flat. This problem did arise also with the multiquadric kernel, but we solved it by using the same trick discussed already in Section 7.2.

• The exponential kernel is not radially-symmetric so that its best discrepancy sequence is not expected to be so, and this is fully consistent with our results in Figure 7.2.

1 +

 1 |x| 2 ) -(D+1)/2 exp(-|x| 2 /2) sup(1 -|x|, 0) D χ (1 + |ξ| 2 ) -(D+1)/2 exp(-|ξ|) exp(-|ξ| 2 /2)

3. 1

 1 Translation-invariant and radially-symmetric kernels on R D

  .1b) Positivity property. Under the assumption (3.1b), let us consider a collection Y = (y 1 , . . . , y N ) in R D and the associated bilinear form K

Figure 4 . 1 :

 41 Figure 4.1: Plots of periodic kernels in two dimensions: exponential, multiquadric, Gaussian, and truncated, respectively.

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: Level sets of periodic kernels in two dimensions: exponential, multiquadric, Gaussian, and truncated, respectively.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Periodic Gaussian kernel in dimension D = 1, 4, 64

Figure 4 . 6 :

 46 Figure 4.6: Periodic truncated kernels in dimension D = 1, 4, 64

Figure 4 . 7 :

 47 Figure 4.7: Transported kernels in two dimensions: exponential, multiquadric, Gaussian, and truncated, respectively.

Proposition 5 . 1 (

 51 Factorization in physical variables). Consider an admissible kernel K = K(x, y) defined on an open set Ω ⊂ R D . Then, for any function ϕ in the Hilbert Space H K (Ω), one has

Proposition 5 . 3 (

 53 Minimizing the error function in spectral variables). Consider an admissible kernel K : Ω × Ω → R together with its Mercer representation (λ i , ζ i ) i≥1 . Then the error function3

Proposition 6 . 1 (

 61 Minimizing the error function in physical variables). Let K : Ω×Ω → R be an admissible kernel defined on some open set Ω, and consider the error function (5.1). Provided the kernel is convex with respect to each variable, the functional Y = (y 1 , . . . , y N ) → (E K (Y, N, D)) 2 is a positive (non-strictly) convex with positive infimum

Proposition 6 . 3 (

 63 Minimizing the error function in spectral variables). Consider an admissible kernel K : Ω × Ω → R and denote by (λ i , ζ i ) i≥1 its Mercer representation. Then for the error function (5.7) after minimization E s,p K (N, D) := inf y 1 ,...,y N ∈Ω E s,p K (Y, N, D)

Proposition 6 . 5 .

 65 Let L a lattice with elementary cell C normalized such that |C| = 1, and consider L * its dual lattice. Let α * ∈ L * → ρ(α * ) be any discrete function satisfying ρ ∈ 1 (L * ), ρ ≥ 0, and ρ(0) = 1. Then the kernel K per (x, y) = α * ∈L * ρ(α * )e 2iπ<x-y,α * > (6.8)

C

  K(x, y) dx = C 0 = C×C K(x, y) dxdy, y ∈ C, since |C| = 1. Consider (6.8) and the error function expressed in Fourier variable (5.8) with s = 1 and p = 2, that is,

R D - 1 n 1 -e 2iπn m 1 1 -e 2iπ n m 1

 1111 computing the solutions to (6.12) for an arbitrary set {α * n ∈ Z D } with size N , we can consider the functional I(Y ) := N n=1 N m=1 e 2iπ<y m ,α * n > 2

Table 1 . 2

 12

: Transported (Tra) and periodic (Per) kernels on [-1, 1] D (sum over α ∈ Z D )

  .1b)

				Matern				multiquadric
	4					4
	2					2
	0	0.05	0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95		0	0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
	-4 -2					-4 -2
		-4	-2	0	2	4	-4	-2	0	2	4
			Gaussian
			0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95

Table 7 .

 7 1: E K (Y, N, D) for each choice of periodic tensorial kernels when N = 512 in dimension D = 1

		random numerical asymptotics
	exponential	0.035	0.002	0.002
	multiquadric	0.036	0.000	0.000
	Gaussian	0.036	0.000	0.000
	truncated	0.036	0.002	0.002

Table 7 .

 7 2: E K (Y, N, D) for each choice of periodic tensorial kernels when N = 512 in dimension D = 16

		random numerical asymptotics
	exponential	0.054	0.034	0.034
	multiquadric	0.055	0.021	0.021
	Gaussian	0.055	0.018	0.018
	truncated	0.056	0.036	0.038
	Table 7.3: E K (Y, N, D) for each choice of periodic tensorial kernels when N = 512 in dimension D = 128
		random numerical asymptotics
	exponential	0.059	0.049	0.043
	multiquadric	0.059	0.037	0.037
	Gaussian	0.059	0.037	0.037

Table 7 .

 7 4: E K (Y, N, D) for periodic tensorial exponential with random points Y

		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0.228 0.365 0.355 0.312	0.304	0.308	0.326	0.319
	N= 32	0.307 0.222 0.216 0.199	0.210	0.228	0.231	0.234
	N= 64	0.154 0.142 0.138 0.139	0.150	0.168	0.163	0.165
	N= 128 0.117 0.088 0.087 0.097	0.111	0.116	0.114	0.115
	N= 256 0.060 0.053 0.061 0.077	0.075	0.081	0.081	0.082
	N= 512 0.035 0.038 0.050 0.049	0.054	0.057	0.057	0.059
	Table 7.5: E K (Y, N, D) for periodic tensorial exponential with numerical points Y
		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0.062 0.126 0.172 0.195	0.211	0.217	0.221	0.223
	N= 32	0.031 0.075 0.114 0.131	0.143	0.149	0.151	0.153
	N= 64	0.016 0.049 0.076 0.090	0.099	0.103	0.106	0.107
	N= 128 0.008 0.030 0.051 0.063	0.069	0.073	0.074	0.077
	N= 256 0.004 0.020 0.034 0.043	0.048	0.051	0.054	0.061
	N= 512 0.002 0.012 0.022 0.030	0.034	0.037	0.042	0.049

Table 7 .

 7

		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0.069 0.143 0.202 0.245	0.288	0.308	0.318	0.323
	N= 32	0.034 0.082 0.129 0.157	0.179	0.207	0.220	0.226
	N= 64	0.017 0.046 0.078 0.102	0.116	0.129	0.147	0.156
	N= 128 0.009 0.026 0.048 0.067	0.077	0.084	0.092	0.105
	N= 256 0.004 0.014 0.029 0.042	0.052	0.056	0.060	0.066
	N= 512 0.002 0.008 0.018 0.027	0.034	0.038	0.040	0.043

6: E K (N, D) for periodic tensorial exponential with the asymptotic formula

Table 7 .

 7 7: E K (Y, N, D) for periodic tensorial multiquadric with random points Y

		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0.249 0.410 0.392 0.315	0.302	0.301	0.325	0.311
	N= 32	0.349 0.245 0.230 0.196	0.201	0.227	0.226	0.234
	N= 64	0.171 0.150 0.140 0.136	0.143	0.169	0.163	0.167
	N= 128 0.130 0.090 0.086 0.093	0.108	0.117	0.115	0.114
	N= 256 0.066 0.050 0.059 0.075	0.075	0.081	0.081	0.082
	N= 512 0.036 0.036 0.049 0.048	0.055	0.057	0.057	0.059

Table 7 .

 7 8: E K (Y, N, D) for periodic tensorial multiquadric with numerical points Y D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0.002 0.078 0.172 0.204	0.261	0.277	0.313	0.306
	N= 32	0.000 0.030 0.095 0.128	0.149	0.198	0.210	0.227
	N= 64	0.000 0.005 0.045 0.081	0.103	0.106	0.140	0.157
	N= 128 0.000 0.001 0.018 0.044	0.067	0.074	0.075	0.098
	N= 256 0.000 0.002 0.007 0.024	0.042	0.051	0.052	0.053
	N= 512 0.000 0.007 0.003 0.014	0.021	0.034	0.037	0.037

Table 7 .

 7 9: E K (N, D) for periodic tensorial multiquadric with the asymptotic formula D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0.004 0.081 0.171 0.207	0.272	0.301	0.314	0.321
	N= 32	0.000 0.027 0.092 0.134	0.148	0.194	0.213	0.223
	N= 64	0.000 0.005 0.044 0.085	0.100	0.105	0.137	0.151
	N= 128 0.000 0.001 0.017 0.043	0.067	0.073	0.075	0.097
	N= 256 0.000 0.000 0.008 0.025	0.043	0.050	0.052	0.053
	N= 512 0.000 0.000 0.003 0.014	0.021	0.034	0.036	0.037

Table 7 .

 7 10: E K (Y, N, D) for periodic Gaussian with random points Y

		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0.279 0.445 0.400 0.316	0.301	0.301	0.325	0.311
	N= 32	0.415 0.267 0.235 0.196	0.200	0.227	0.226	0.234
	N= 64	0.205 0.159 0.138 0.136	0.142	0.168	0.163	0.167
	N= 128 0.152 0.090 0.085 0.092	0.107	0.117	0.115	0.114
	N= 256 0.075 0.047 0.059 0.074	0.075	0.081	0.081	0.082
	N= 512 0.036 0.034 0.046 0.048	0.055	0.057	0.057	0.059

Table 7 .

 7 11: E K (Y, N, D) for periodic Gaussian with numerical points Y

		D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128
	N= 16	0 0.008 0.164 0.191	0.258	0.276	0.313	0.306
	N= 32	0 0.000 0.051 0.123	0.145	0.197	0.209	0.227
	N= 64	0 0.000 0.013 0.075	0.099	0.105	0.140	0.157
	N= 128	0 0.000 0.003 0.033	0.064	0.072	0.075	0.098
	N= 256	0 0.002 0.000 0.019	0.041	0.050	0.052	0.053
	N= 512	0 0.008 0.000 0.008	0.018	0.033	0.036	0.037

Table 7 .

 7 12: E K (N, D) for periodic Gaussian with the asymptotic formula D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0 0.018 0.145 0.198	0.270	0.300	0.314	0.321

Table 7 .

 7 13: E K (Y, N, D) for periodic tensorial truncated with random points Y D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0.279 0.398 0.354 0.327	0.315	0.317	0.324	0.328
	N= 32	0.382 0.240 0.216 0.210	0.225	0.225	0.232	0.231
	N= 64	0.189 0.148 0.142 0.149	0.159	0.163	0.163	0.163
	N= 128 0.142 0.092 0.090 0.105	0.113	0.116	0.115	0.115
	N= 256 0.071 0.051 0.063 0.081	0.079	0.082	0.082	0.082
	N= 512 0.038 0.040 0.055 0.052	0.056	0.057	0.057	0.058

Table 7 .

 7 14: E K (Y, N, D) for periodic tensorial truncated with numerical points Y D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0.062 0.100 0.176 0.209	0.233	0.276	0.303	0.315
	N= 32	0.031 0.058 0.116 0.140	0.156	0.173	0.196	0.214
	N= 64	0.016 0.035 0.079 0.096	0.107	0.119	0.125	0.139
	N= 128 0.007 0.021 0.049 0.067	0.075	0.081	0.085	0.090
	N= 256 0.004 0.013 0.030 0.046	0.052	0.056	0.058	0.061
	N= 512 0.002 0.010 0.021 0.032	0.036	0.039	0.040	0.042

Table 7 .

 7 15: E K (N, D) for periodic tensorial truncated with the asymptotic formula D= 1 D= 2 D= 4 D= 8 D= 16 D= 32 D= 64 D= 128

	N= 16	0.062 0.127 0.217 0.289	0.314	0.322	0.325	0.327
	N= 32	0.031 0.077 0.133 0.188	0.218	0.227	0.230	0.231
	N= 64	0.016 0.042 0.086 0.114	0.148	0.159	0.162	0.163
	N= 128 0.007 0.023 0.054 0.073	0.096	0.110	0.114	0.115
	N= 256 0.004 0.013 0.034 0.050	0.059	0.075	0.080	0.081
	N= 512 0.002 0.007 0.022 0.034	0.038	0.049	0.055	0.057

which is sometimes also refered to as the Matérn kernel

δ ij being the Kronecker symbol

This expression is equivalent to the one given in the introduction.

Acknowledgments. The first author (PLF) gratefully acknowledges support from the Simons Center for Geometry and Physics, Stony Brook University, as well as from the Innovative Training Network (ITN) grant 642768 (ModCompShocks). This paper was written during the Academic year 2018-2019 when he was visiting the Courant Institute for Mathematical Sciences, New York University.

which is positive. Hence, being the product of two admissible kernels, we see that the kernel K loc is admissible. The transported kernel K tran is also admissible since for all relevant ψ Ω×Ω K tran (y, y )ψ(y)ψ(y )dydy =

Proof of Proposition 3.3. Since ρ 1/2 ∈ 2 (L * ) it is clear that the expression Proof of Proposition 5.1. In view of (2.4) and the identity ϕ(x) = ϕ, K(•, x) H K (Ω) , we find

In agreement with the statement of the proposition, we write

and, by the Cauchy-Schwarz inequality,

. Hence, the integration error splits into two contributions, as expected. By expanding the integrand using (2.4) we find the equivalent expression

Moreover, in view of K(x, y) = 1 2 K(x, x) + K(y, y) -D(x, y) , the derivation of (5.2a) is completed.

Proof of Proposition 5.3. The expression KN (x, x ) = 1≤n≤N λnζn(x)ζn(x ) defines an admissible kernel, for which we have the following orthogonal decomposition: HK (Ω) = HK N (Ω) ⊕ H ⊥ K N (Ω). Consider the minimization problem

Proof of Proposition 5.4. We consider the Mercer representation and use the relations TK ζi = λiζi and with eK = i eK , ζi 2 (Ω) ζi. We have

and, by Hölder inequality, we obtain

Proof of Proposition 5.5 . Provided supp(ϕ) ⊂ C, we have

where the error function e is defined as the Fourier transform of e :

Proof of Proposition 6.1. Let us set e := 1

Thus EK (Y, N, D) 2 ≥ e -eY 2 H K > 0, and, as we are going to check, a minimizer makes the term eY 2

In view of the definition (5.7) and the symmetry property ∂yn K(y n , y m ) = ∂yn K(y m , y n ), we obtain

. By computing second-order derivatives, we see that Hessian ∂ y n d ∂yn e (EK (Y, N, D)) 2 reads as follows (d, e ∈ 1, . . . , D ):

Moreover, recalling (2.10b), we compute

Hence, we arrive at

In particular, the minimum vanishes when Ω δ n (x)dx = 1 N , and this establishes (6.2). Hence, suppose that this relation is true, then from (2.10b) we deduce

) m is a stochastic matrix, which is (6.2).

Proof of Proposition 6.2. We begin with the case s = and p = 2, and consider an arbitrary (ordered) sequence of points, say 0 < y 1 < . . . < y N < 1. We rewrite our functional (5.2a) as

In view of Proposition 6.1, any minimizer satisfies the system of equations (n = 1, 2, . . . , N )

y n + 1 N 2 (2n -N -1) = 0, which we put in the form 1 N 1≤n≤N y n -1 2 + 1 N (n -1/2) = y n . This is a linear algebraic system, which is readily solved explicitly; this leads us to (6.4). We check immediately that 1 N 1≤n≤N y n = 1 2 and 1≤n≤N (y n ) 2 = N 3 -1 12N , and, finally, we evaluate the functional to be (EK (Y, D, N ))