On the Equivalence between Kernel Quadrature Rules and Random Feature Expansions

Introduction

The numerical computation of high-dimensional integrals is one of the core computational tasks in many areas of machine learning, signal processing and more generally applied mathematics, in particular in the context of Bayesian inference (Gelman, 2004), or the study of complex systems (Robert and Casella, 2005). In this paper, we focus on quadrature rules, that aim at approximating the integral of a certain function from only the (potentially noisy) knowledge of the function values at as few as possible well-chosen points. Key situations that remain active areas of research are problems where the measurable space where the function is defined on is either high-dimensional or structured (e.g., presence of discrete structures, or graphs). For these problems, techniques based on positive definite kernels have emerged as having the potential to efficiently deal with these situations, and to improve over plain Monte-Carlo integration (O’Hagan, 1991; Rasmussen and Ghahramani, 2003; Huszár and Duvenaud, 2012; Oates and Girolami, 2015). In particular, the quadrature problem may be cast as the one of approximating a fixed element, the mean element (Smola et al., 2007), of a Hilbert space as a linear combination of well chosen elements, the goal being to minimize the number of these factors as it corresponds to the required number of function evaluations.

A seemingly unrelated problem on positive definite kernels have recently emerged, namely the representation of the corresponding infinite-dimensional feature space from random sets of features. If a certain positive definite kernel between two points may be represented as the expectation of the product of two random one-dimensional (typically non-linear) features computed on these two points, the full kernel (and hence its feature space) may be approximated by sufficiently many random samples, replacing the expectation by a sample average (Neal, 1995; Rahimi and Recht, 2007; Huang et al., 2006). When using these random features, the complexity of a regular kernel method such as the support vector machine or ridge regression goes from quadratic in the number of observations to linear in the number of observations, with a constant proportional to the number of random features, which thus drives the running time complexity of these methods.

In this paper, we make the following contributions:

After describing the functional analysis framework our analysis is based on and presenting many examples in Section 2, we show in Section 3 that these two problems are strongly related; more precisely, optimizing weights in kernel-based quadrature rules can be seen as decomposing a certain function in a special class of random features for a particular decomposition that always exists for all positive definite kernels on a measurable space.

We provide in Section 4 a theoretical analysis of the number of required samples for a given approximation error, leading to both upper and lower bounds that are based solely on the eigenvalues of the associated integral operator and match up to logarithmic terms. In particular, we show that the upper bound may be obtained as independent and identically distributed samples from a specific non-uniform distribution, while the lower bound if valid for any set of points.

Applying our results to kernel quadrature, while our results are fairly general, we recover known upper and lower bounds for the special cases of Sobolev spaces (Section 4.4). Moreover, our results extend to the more general problem of full function approximations (beyond simply computing an integral), with results in $L_{2}$- and $L_{\infty}$-norm that match known results for special cases (Section 5).

Applying our results to random feature expansions, we show in Section 4.5 an improvement of the number of random features needed for preserving the generalization guarantees for learning with Lipshitz-continuous losses.

Random Feature Expansions of Positive Definite Kernels

Throughout this paper, we consider a topological space ${\mathcal{X}}$ equipped with a Borel probability measure $d\rho$, which we assume to have full support. This naturally defines the space of square-integrable functionsFor simplicity and following most of the literature on kernel methods, we identify functions and their equivalence classes for the equivalence relationship of being equal except for a zero-measure (for $d\rho$) subset of ${\mathcal{X}}$..

Membership of kernel evaluations: for any $x\in{\mathcal{X}}$, the function $k(\cdot,x):y\mapsto k(y,x)$ is an element of ${\mathcal{F}}$.

Reproducing property: for all $f\in{\mathcal{F}}$ and $x\in{\mathcal{X}}$, $f(x)=\langle f,k(\cdot,x)\rangle_{\mathcal{F}}$. In other words, function evaluations are equal to dot-products with a specific element of ${\mathcal{F}}$.

Moreover, throughout the paper, we assume that the function $x\mapsto k(x,x)$ is integrable with respect to $d\rho$ (which is weaker than $\sup_{x\in\mathcal{X}}k(x,x)<\infty$). This implies that ${\mathcal{F}}$ is a subset of $L_{2}(d\rho)$; that is, functions in the RKHS ${\mathcal{F}}$ are all square-integrable for $d\rho$. In general, ${\mathcal{F}}$ is strictly included in $L_{2}(d\rho)$, but, in this paper, we will always assume that it is dense in $L_{2}(d\rho)$, that is, any function in $L_{2}(d\rho)$ may be approximated arbitrarily closely by a function in ${\mathcal{F}}$. Finally, for simplicity of our notation (to make sure that the sequence of eigenvalues of integral operators is infinite) we will always assume that $L_{2}(d\rho)$ is infinite-dimensional, which excludes finite sets for ${\mathcal{X}}$. Note that the last two assumptions (denseness and infinite dimensionality) can easily be relaxed.

Reproducing kernel Hilbert spaces are often studied through a specific integral operator which leads to an isometry with $L_{2}(d\rho)$ (Smale and Cucker, 2001). Let $\Sigma:L_{2}(d\rho)\to L_{2}(d\rho)$ be defined as

Since $\int_{\mathcal{X}}k(x,x)d\rho(x)$ is finite, $\Sigma$ is self-adjoint, positive semi-definite and trace-class (Simon, 1979). Given that $\Sigma f$ is a linear combination of kernel functions $k(\cdot,y)$, it belongs to ${\mathcal{F}}$. More precisely, since we have assumed that ${\mathcal{F}}$ is dense in $L_{2}(d\rho)$, $\Sigma^{1/2}$, which is the unique positive self-adjoint square root of $\Sigma$, is a bijection from $L_{2}(d\rho)$ to our RKHS ${\mathcal{F}}$; that is, for any $f\in{\mathcal{F}}$, there exists a unique $g\in L_{2}(d\rho)$ such that $f=\Sigma^{1/2}g$ and $\|f\|_{\mathcal{H}}=\|g\|_{L_{2}(d\rho)}$ (Smale and Cucker, 2001). This justifies the notation $\Sigma^{-1/2}f$ for $f\in{\mathcal{F}}$ and means that $\Sigma^{1/2}$ is an isometry from $L_{2}(d\rho)$ to ${\mathcal{F}}$; in other words, for any functions $f$ and $g$ in ${\mathcal{F}}$, we have:

This justifies the view of ${\mathcal{F}}$ as the subspace of functions $f\in L_{2}(d\rho)$ such that $\|\Sigma^{-1/2}f\|_{L_{2}(d\rho)}^{2}$. This relationship is even more transparent when considering a spectral decomposition of $\Sigma$.

Mercer decomposition.

From extensions of Mercer’s theorem (König, 1986), there exists an orthonormal basis $(e_{m})_{m\geqslant 1}$ of $L_{2}(d\rho)$ and a summable non-increasing sequence of strictly positive eigenvalues $(\mu_{m})_{m\geqslant 1}$ such that $\Sigma e_{m}=\mu_{m}e_{m}$. Note that since we have assumed that ${\mathcal{F}}$ is dense in $L_{2}(d\rho)$, there are no zero eigenvalues.

Since $\Sigma^{1/2}$ is an isometry from $L_{2}(d\rho)$ to ${\mathcal{F}}$, $(\mu_{m}^{1/2}e_{m})_{m\geqslant 1}$ is an orthonormal basis of ${\mathcal{F}}$. Moreover, we can use the eigendecomposition to characterize elements of ${\mathcal{F}}$ as the functions in $L_{2}(d\rho)$ such that

is finite. In other words, once projected in the orthonormal basis $(e_{m})_{m\geqslant 1}$, elements $f$ of ${\mathcal{F}}$ correspond to a certain decay of its decomposition coefficients $(\langle f,e_{m}\rangle_{L_{2}(d\rho)})_{m\geqslant 1}$.

Finally, by decomposing the function $k(\cdot,y):x\mapsto k(x,y)$, we obtain the Mercer decomposition:

Properties of the spectrum.

The sequence of eigenvalues $(\mu_{m})_{m\geqslant 1}$ is an important quantity that appears in the analysis of kernel methods (Hastie and Tibshirani, 1990; Caponnetto and De Vito, 2007; Harchaoui et al., 2008; Bach, 2013; El Alaoui and Mahoney, 2014). It depends both on the kernel $k$ and the chosen distribution $d\rho$.

Some modifications of the kernel $k$ or the distribution $d\rho$ lead to simple behaviors for the spectrum. For example, if we have a second distribution so that $\frac{d\rho^{\prime}}{d\rho}$ is upper-bounded by a constant $c$, then, as a consequence of the Courant-Fischer minimax theorem (Horn and Johnson, 2012), the eigenvalues for $d\rho^{\prime}$ are less than than $c$ times that the ones for $d\rho$. Similarly, if the kernel $k^{\prime}$ is such that $ck-k^{\prime}$ is a positive definite kernel, then we have a similar bound between eigenvalues.

In this paper, for any strictly positive $\lambda$, we will also consider the quantity $m^{\ast}(\lambda)$ equal to the number of eigenvalues $\mu_{m}$ that are greater than or equal to $\lambda$. Since we have assume that the sequence $m$ is non-increasing, we have $m^{\ast}(\lambda)=\max\{m\geqslant 1,\ \mu_{m}\geqslant\lambda\}$. This is a left-continuous non-increasing function, that tends to $+\infty$ when $\lambda$ tends to zero (since we have assumed that there are infinitely many strictly positive eigenvalues), and characterizes the sequence $(\mu_{m})_{m\geqslant 1}$, as we can recover $\mu_{m}$ as $\mu_{m}=\sup\{\lambda\geqslant 0,\ m^{\ast}(\lambda)\geqslant m\}$.

Potential confusion with covariance operator.

Note that the operator $\Sigma$ is a self-adjoint operator on $L_{2}(d\rho)$. It should not be confused with the (non-centered) covariance operator $C$ (Baker, 1973), which is a self-adjoint operator on a different space, namely the RKHS ${\mathcal{F}}$, defined by $\langle g,Cf\rangle_{\mathcal{F}}=\int_{{\mathcal{X}}}f(x)g(x)d\rho(x)$. Given that $\Sigma^{1/2}$ is an isometry from $L_{2}(d\rho)$ to ${\mathcal{F}}$, the operator $C$ may also be used to define an operator on $L_{2}(d\rho)$, which happens to be exactly $\Sigma$. Thus, the two operators have the same eigenvalues. Moreover, we have, for any $y\in{\mathcal{X}}$:

that is, $C$ is equal to the restriction of $\Sigma$ on ${\mathcal{F}}$.

2 Kernels as expectations

Such additional structure allows to give an explicit characterization of the RKHS ${\mathcal{F}}$ in terms of the features $\varphi$. In terms of operators, the function $\varphi$ leads to a specific square-root of the integral operator $\Sigma$ defined in Section 2.1 (which is not the positive self-adjoint square-root $\Sigma^{1/2}$).

We consider the bounded linear operator $T:L_{2}(d\tau)\to L_{2}(d\rho)$ defined as

Given $T:L_{2}(d\tau)\to L_{2}(d\rho)$, the adjoint operator $T^{\ast}:L_{2}(d\rho)\to L_{2}(d\tau)$ is the unique operator such that $\langle g,T^{\ast}f\rangle_{L_{2}(d\tau)}=\langle Tg,f\rangle_{L_{2}(d\rho)}$ for all $f,g$. Given the definition of $T$ in Eq. (2), we simply inverse the role of ${\mathcal{V}}$ and ${\mathcal{X}}$ and have:

that is we have an expression of the integral operator $\Sigma$ as $\Sigma=TT^{\ast}$. Thus, the decomposition of the kernel $k$ as an expectation corresponds to a particular square root $T$ of the integral operator—there are many possible choices for such square roots, and thus many possible expansions like Eq. (1). It turns out that the positive self-adjoint square root $\Sigma^{1/2}$ will correspond to the equivalence with quadrature rules (see Section 3.2).

Decomposition of functions in ℱℱ{\mathcal{F}}.

Since $\Sigma=TT^{\ast}$ and $\Sigma^{1/2}$ is an isometry between $L_{2}(d\rho)$ and ${\mathcal{F}}$, we can naturally expressed any elements of ${\mathcal{F}}$ through the operator $T$ and thus the features $\varphi$.

As a linear operator, $T$ defines a bijection from the orthogonal of its null space $({\rm Ker}\ T)^{\perp}\subset L_{2}(d\tau)$ to its image ${\rm Im}(T)\subset L_{2}(d\rho)$, and this allows to define uniquely $T^{-1}f\in({\rm Ker}\ T)^{\perp}$ for any $f\in{\rm Im}(T)$, and a dot-product on ${\rm Im}(T)$ as

As shown by Bach (2014, App. A), ${\rm Im}(T)$ turns out to be equal to our RKHSThe proof goes as follows: (a) for any $y\in{\mathcal{X}}$, $k(\cdot,y)$ can be expressed as $\int_{{\mathcal{V}}}\varphi(v,y)\varphi(v,\cdot)d\tau(v)=T\varphi(\cdot,y)$ and thus belongs to ${\rm Im}(T)$; (b) for any $f\in{\rm Im}(T)$, and $y\in{\mathcal{X}}$, we have $\langle f,k(\cdot,y)\rangle_{{\rm Im}(T)}=\langle T^{-1}f,\varphi(\cdot,y)\rangle_{L_{2}(d\tau)}=(TT^{-1}f)(y)=f(y)$, that is, the reproducing property is satisfied. These two properties are characteristic of ${\mathcal{F}}$.. Thus, the norm $\|f\|_{\mathcal{F}}^{2}$ for $f\in{\mathcal{F}}$ is equal to the squared $L_{2}$-norm of $T^{-1}f\in({\rm Ker}\ T)^{\perp}$, which is itself equal to the minimum of $\|g\|^{2}_{L_{2}(d\tau)}$ over all $g$ such that $Tg=f$. The resulting $g$ may also be defined through pseudo-inverses.

In other words, a function $f\in L_{2}(d\rho)$ is in ${\mathcal{F}}$ if and only if it may be written as

Singular value decomposition.

The operator $T$ is an Hilbert-Schmidt operator, to which the singular value decopomposition can be applied (Kato, 1995). That is, there exists an orthonormal basis $(f_{m})_{m\geqslant 1}$ of $({\rm Ker}\ T)^{\perp}\subset L_{2}(d\tau)$, together with the orthonormal basis $(e_{m})_{m\geqslant 1}$ of $L_{2}(d\rho)$ which we have from the eigenvalue decomposition of $\Sigma=TT^{\ast}$, such that $Tf_{m}=\mu_{m}^{1/2}e_{m}$. Moreover, we have:

with a convergence in $L_{2}(d\tau\otimes d\rho)$. This extends the Mercer decomposition of the kernel $k(x,y)$.

Integral operator as an expectation.

Given the expansion of the kernel $k$ in Eq. (1), we may express the integral operator $\Sigma$ as follows, explicitly as an expectation:

where $a\otimes_{L_{2}(d\rho)}b$ is the operator $L_{2}(d\rho)\to L_{2}(d\rho)$ so that $(a\otimes_{L_{2}(d\rho)}b)f=\langle b,f\rangle_{L_{2}(d\rho)}a$. This will be useful to define empirical versions, where the integral over $d\tau$ will be replaced by a finite average.

3 Examples

In this section, we provide examples of kernels and usual decompositions. We first start by decompositions that always exist, then focus on specific kernels based on Fourier components.

The Mercer decomposition provides an expansion for all kernels, as follows:

Periodic kernels on [0,1]01[0,1].

We consider ${\mathcal{X}}=$ and translation-invariant kernels $k(x,y)$ of the form $k(x,y)=t(x-y)$, where $t$ is a square-integrable 1-periodic function. These kernels are positive definite if and only if the Fourier series of $t$ is non-negative (Wahba, 1990). An orthonormal basis of $L_{2}()$ is composed of the constant function $c_{0}:x\mapsto 1$ and the functions $c_{m}:x\mapsto\sqrt{2}\cos 2\pi mx$ and $s_{m}:x\mapsto\sqrt{2}\sin 2\pi mx$. A kernel may thus be expressed as

A particularly interesting example is obtained through derivatives of $f$. If $f$ is differentiable and has a derivative $f^{\prime}\in L_{2}()$, then, on the Fourier series coefficients of $f$, taking the derivative corresponds to multiplying the two $m$-th coefficients by $2\pi m$ and swapping them. Sobolev spaces for periodic functions on $$(i.e., such that$f(0)\!=\!f(1)$) are defined through integrability of derivatives (Adams and Fournier, 2003). In the Hilbert space set-up, a function$f$belongs to the Sobolev space of order$s$if one can define a$s$-th order square-integrable derivative in$L_{2}$(for the Lebesgue measure, which happens to be equal to$d\rho$), that is,$f^{(s)}\in L_{2}()$. The Sobolev squared norm is then defined as any positive linear combination of the quadratic forms$\int_{0}^{1}f^{(t)}(x)^{2}dx$,$t\in\{0,\dots,s\}$, with non-zero coefficients for$t=0$and$t=s$(all of these norms are then equivalent). If only using$t=0$and$t=s$with non-zero coefficients, we need$\nu_{0}^{-1}=1$and$\nu_{m}^{-1}=1+m^{2s}$. An equivalent (i.e, with upper and lower bounded ratios) sequence is obtained by replacing$\nu_{m}=(1+m^{2s})^{-1}$by$\nu_{m}=m^{-2s}$, leading to a closed-form formula:

where $\{x-y\}$ denotes the fractional part of $x-y$, and $B_{2s}$ is the $2s$-th Bernoulli polynomial (Wahba, 1990). The RKHS ${\mathcal{F}}$ is then the Sobolev space of order $s$ on $$, with a norm equivalent to any of the family of Sobolev norms; it will be used as a running example throughout this paper.

In order to extend to $d>1$, we may consider several extensions as described by Oates and Girolami (2015), and compute the resulting eigenvalues of the integral operators. For simplicity, we consider the Sobolev space on $$, with$\nu_{0}=1$and$\nu_{m}^{-1}=m^{2s}$for$m>0$. The first possibility to extend to$^{d}$is to take a kernel which is simply the pointwise product of individual kernels on$$. That is, if$k(x,y)$is the kernel on$$, define$K(X,Y)=\prod_{j=1}^{d}k(x_{j},y_{j})$between$X$and$Y$in$^{d}$. As shown in Appendix A, this leads to eigenvalue decays bounded by$(\log m)^{2s(d-1)}m^{-2s}$, and thus up to logarithmic terms at the same speed$m^{-2s}$as$d=1$. While this sounds attractive in terms of generalization performance, it corresponds to a space a function which is not a Sobolev space in$d$dimensions. That is the associated squared norm on$f$would be equivalent to a linear combination of squared$L_{2}$-norm of partial derivatives

for all $t_{1},\dots,t_{d}$ in $\{0,\dots,s\}$. This corresponds to functions which have square-integrable partial derivatives with all individual orders less than $s$. All values of $s\geqslant 1$ are allowed and lead to an RKHS.

This is thus to be contrasted with the usual multi-dimensional Sobolev space which is composed of functions which have square-integrable partial derivatives with all orders $(t_{1},\dots,t_{d})$ with sum $t_{1}+\cdots+t_{d}$ less than $s$. Only $s>d/2$ is then allowed to get an RKHS. The Sobolev norm is then of the form

In the expansion on the $d$-th order tensor product of the Fourier basis, the norm above is equivalent to putting a weight on the element $(m_{1},\dots,m_{d})$ asymptotically equivalent to \big{(}\sum_{j=1}^{d}m_{j}\big{)}^{2s}, which thus represent the inverse of the eigenvalues of the corresponding kernel for the uniform distribution $d\rho$ (this is simply an explicit Mercer decomposition). Thus, the number of eigenvalues which are greater than $\lambda$ grows as the number of $(m_{1},\dots,m_{d})$ such that their sum is less than $\lambda^{-1/(2s)}$, which itself is less than a constant times $\lambda^{-d/(2s)}$ (see a proof in Appendix A). This leads to an eigenvalue decay of $m^{-2s/d}$, which is much worse because of the term in $1/d$ in the exponent.

For these kernels, the decay of eigenvalues has been well-studied by Widom (1963), who relates the decay of eigenvalues to the tails of the distribution $d\rho$ and the decay of the Fourier transform of $t$. For example, for the Gaussian kernel where $k(x,y)=\exp(-\alpha\|x-y\|_{2}^{2})$, on sub-Gaussian distributions, the decay of eigenvalues is geometric, and for kernels leading to Sobolev spaces of order $s$, such as the Matern kernel (Furrer and Nychka, 2007), the decay is of the form $m^{-2s/d}$. See also examples by Birman and Solomyak (1977); Harchaoui et al. (2008).

Finally, note that in terms of computation, there are extensions to avoid linear complexity in $d$ (Le et al., 2013).

Kernels on hyperspheres.

As shown by Smola et al. (2001); Bach (2014), the equivalent of Fourier series (which corresponds to $d=1$) is then the basis of spherical harmonics, which is organized by integer frequencies $k\geqslant 1$; instead of having 2 basis vectors (sine and cosine) per frequency, there are $O(k^{d-1})$ of them. As shown by Bach (2014, page 44), we have an explicit expansion of $k(x,y)$ in terms of spherical harmonics, leading to a sequence of eigenvalues equal to $k^{-d-2s-1}$ on the entire subspace associated with frequency $k$. Thus, by taking multiplicity into account, after $\sum_{j=1}^{k}j^{d-1}\approx k^{d}$ (up to constants) eigenvalues, we have an eigenvalue of $k^{-d-2s-1}$; this leads to an eigenvalue decay (where all eigenvalues are ordered in decreasing order and we consider the $m$-th one) as $(m^{1/d})^{-d-2s-1}=m^{-1-1/d-2s/d}$.

4 Approximation from randomly sampled features

Given the formulation of $k$ as an expectation in Eq. (1), it is natural to consider sampling $n$ elements $v_{1},\dots,v_{n}\in{\mathcal{V}}$ from the distribution $d\tau$ and define the kernel approximation

which defines a finite-dimensional RKHS $\hat{{\mathcal{F}}}$.

In this paper, we rather consider approximations of functions in ${\mathcal{F}}$ by functions in $\hat{{\mathcal{F}}}$, the RKHS associated with $\hat{k}$. A key difficulty is that in general $\hat{{\mathcal{F}}}$ is not even included in $\mathcal{F}$, and therefore, we cannot use the norm in ${\mathcal{F}}$ to characterize approximations. In this paper, we choose the $L_{2}$-norm associated with the probability measure $d\rho$ on ${\mathcal{X}}$ to characterize the approximation. Given $f\in{\mathcal{F}}$ with norm $\|f\|_{\mathcal{F}}$ less than one, we look for a function $\hat{f}\in\hat{{\mathcal{F}}}$ of the smallest possible norm and so that $\|f-\hat{f}\|_{L_{2}(d\rho)}$ is as small as possible.

Note that the measure $d\tau$ is associated to the kernel $k$ and the random features $\varphi$, while the measure $d\rho$ is associated to the way we want to measure errors (and leads to a specific defintion of the integral operator $\Sigma$).

Goals.

We thus aim at sampling $n$ points $v_{1},\dots,v_{n}\in{\mathcal{V}}$ from a distribution with density $q$ with respect to $d\tau$. Then the kernel approximation using importance weights is equal to

(so that the law of large numbers leads to an approximation converging to $k$), and we thus aim at minimizing \Big{\|}\sum_{i=1}^{n}\frac{\beta_{i}}{q(v_{i})^{1/2}}\varphi(v_{i},\cdot)-\int_{\mathcal{V}}g(v)\varphi(v,\cdot)d\tau(v)\Big{\|}_{L_{2}(d\rho)}, with $n\|\beta\|_{2}^{2}$ (which represents the norm of the approximation in $\hat{{\mathcal{F}}}$ because of our importance weights are taken into account) as small as possible.

Quadrature in RKHSs

Following Smola et al. (2007), the error may be expressed using the reproducing property as:

and by Cauchy-Schwarz inequality its supremum over $\|h\|_{\mathcal{F}}\leqslant 1$ is equal to

In this paper, we explore sampling points $x_{i}$ from a probability distribution on ${\mathcal{X}}$ with density $q$ with respect to $d\rho$. Note that when $g$ is a constant function, it is sometimes required that the coefficients $\alpha$ are non-negative and sum to a fixed constant (so that constant functions are exactly integrated). We will not pursue this here as our theoretical results do not accommodate such constraints (see, e.g., Chen et al., 2010; Bach et al., 2012, and references therein).

In practice, independent (but not necessarily identically distributed) noise $\varepsilon_{i}$ may be present with variance $\sigma^{2}(x_{i})$. Then, the worst (with respect to $\|h\|_{\mathcal{F}}\leqslant 1$) expected (with respect to the noise) squared error is

2 Reformulation as random features

For any $x\in{\mathcal{X}}$, the function $k(\cdot,x)$ is in ${\mathcal{F}}$, and since we have assumed that $\Sigma^{1/2}$ is an isometry from $L_{2}(d\rho)$ to ${\mathcal{F}}$, there exists a unique element, which we denote $\psi(\cdot,x)$, of $L_{2}(d\rho)$ such that $\Sigma^{1/2}\psi(\cdot,x)=k(\cdot,x)$. Given the Mercer decomposition $k(\cdot,x)=\sum_{m\geqslant 1}\mu_{m}e_{m}(x)e_{m}$, we have the expansion $\psi(\cdot,x)=\sum_{m\geqslant 1}\mu_{m}^{1/2}e_{m}(x)e_{m}$ (with convergence in the $L_{2}$-norm for the measure $d\rho\otimes d\rho$; note that we do not assume that $\mu_{m}^{1/2}$ is summable), and thus we may consider $\psi$ as a symmetric function. Note that $\psi$ may not be easy to compute in many practical cases (except for some periodic kernels on $$).

We thus have for $(x,y)\in{\mathcal{X}}\times{\mathcal{X}}$:

That is, the kernel $k$ may always be written as an expectation. Moreover, we have the quadrature error in Eq. (7) equal to (again using the isometry $\Sigma^{1/2}$ from $L_{2}(d\rho)$ to ${\mathcal{F}}$):

which is exactly an instance of the approximation result in Eq. (6) with ${\mathcal{V}}={\mathcal{X}}$ and $\varphi=\psi$, that is the random feature is indexed by the same set ${\mathcal{X}}$ as the kernel. Thus, the quadrature problem, that is finding points $x_{i}$ and weights $(\alpha_{i})$ to get the best possible error over all functions of the unit ball of ${\mathcal{F}}$, is a subcase of the random feature problem for a specific expansion. Note that this random decomposition in terms of $\psi$ is always possible (although not in closed form in general).

As shown in Section 2.2, random feature expansions correspond to square-roots of the integral operator $\Sigma:L_{2}(d\rho)\to L_{2}(d\rho)$ as $\Sigma=TT^{\ast}$. Among the many possible square roots, the quadrature case corresponds exactly to the positive self-adjoint square root $T=\Sigma^{1/2}$. In this situation, the basis $(f_{m})_{m\geqslant 1}$ of the singular value decomposition of $T=\Sigma^{1/2}$ is equal to $(e_{m})_{m\geqslant 1}$, recovering the expansion $\psi(x,y)=\sum_{m\geqslant 1}\mu_{m}^{1/2}e_{m}(x)e_{m}(y)$ which we have seen above.

In this important situation, we have two different expansions: the one based on Fourier features, where the random variable indexing the one-dimensional feature is a frequency, while for the one based on the square root $\psi$, the random variable is a spatial variable in ${\mathcal{X}}$. As we show in Section 4, our results are independent of the chosen expansions and thus apply to both. However, (a) when the goal is to do quadrature, we need to use $\psi$, and (b) in general, the decomposition based on Fourier features can be easily computed once samples are obtained, while for most kernels, $\psi(x,y)$ does not have any closed-form simple expression. In Section 6, we provide a simple example with ${\mathcal{X}}=$ where the two decompositions are considered.

Goals.

In order to be able to make the parallel with random feature approximations, we consider importance-weighted coefficients $\beta_{i}=\alpha_{i}q(x_{i})^{1/2}$, and we thus aim at minimizing the approximation error

3 Relationship with column sampling

The problem of quadrature is related to the problem of column sampling. Given $n$ observations $x_{1},\dots,x_{n}\in{\mathcal{X}}$, the goal of column-sampling methods is to approximate the $n\times n$ matrix of pairwise kernel evalulations, the so-called kernel matrix, from a subset of its columns. It has appeared under many names: Nyström method (Williams and Seeger, 2001), sparse greedy approximations (Smola and Schölkopf, 2000), incomplete Cholesky decomposition (Fine and Scheinberg, 2001), Gram-Schmidt orthonormalization (Shawe-Taylor and Cristianini, 2004) or CUR matrix decompositions (Mahoney and Drineas, 2009).

While column sampling has typically been analyzed for a fixed kernel matrix, it has a natural extension which is related to quadrature problems: selecting $n$ points $x_{1},\dots,x_{n}$ from ${\mathcal{X}}$ such that the projection of any element of the RKHS ${\mathcal{F}}$ onto the subspace spanned by $k(\cdot,x_{i})$, $i=1,\dots,n$ is as small as possible. Natural functions from ${\mathcal{F}}$ are $k(\cdot,x)$, $x\in{\mathcal{X}}$, and thus the goal is to minimize, for such $x\in{\mathcal{X}}$,

In the usual sampling approach, several points are considered for testing the projection error, and it is thus natural to consider the criterion averaged through the measure $d\rho$, that is:

In fact, when $d\rho$ is supported on a finite set, this formulation is equivalent to minimizing the nuclear norm between the kernel matrix and its low-rank approximation. There are thus several differences and similarities between recent work on column sampling (Bach, 2013; El Alaoui and Mahoney, 2014) and the present paper on quadrature rules and random features:

Different error measures: The column sampling approach corresponds to a function $g$ in Eq. (7) which is a Dirac function at the point $x$, and is thus not in $L_{2}(d\rho)$. Thus the two frameworks are not equivalent.

Approximation vs. prediction: The works by Bach (2013); El Alaoui and Mahoney (2014) aim at understanding when column sampling leads to no loss in predictive performance within a supervised learning framework, while the present paper looks at approximation properties, mostly regardless of any supervised learning problem, except in Section 4.5 for random features (but not for quadrature).

Lower bounds: In Section 4.3, we provide explicit lower bounds of approximations, which are not available for column sampling.

Similar sampling issues: In the two frameworks, points $x_{1},\dots,x_{n}\in{\mathcal{X}}$ are sampled i.i.d. with a certain distribution $q$, and the best choice depends on the appropriate notion of leverage scores (Mahoney, 2011), while the standard uniform distribution leads to an inferior approximation result. Moreover, the proof techniques are similar and based on concentration inequalities for operators, here in Hilbert spaces rather in finite dimensions.

4 Related work on quadrature

Many methods have been designed for the computation of integrals of a function given evaluations at certain well-chosen points, in most cases when $g$ is constant equal to one. We review some of these below.

When the underlying set ${\mathcal{X}}$ is a compact interval of the real line, several methods exists, such as the trapezoidal or Simpson’s rules, which are based on interpolation between the sample points, and for which the error decays as $O(1/n^{2})$ and $O(1/n^{4})$ for functions with uniformly bounded second or fourth derivatives (Cruz-Uribe and Neugebauer, 2002).

Gaussian quadrature is another class of methods for one-dimensional integrals: it is based on a basis of orthogonal polynomials for $L_{2}(d\rho)$ where $d\rho$ is a probability measure supported in an interval, and their zeros (Hildebrand, 1987, Chap. 8). This leads to quadrature rules which are exact for polynomials of degree $2n-1$ but error bounds for non-polynomials rely on high-order derivatives, although the empirical performance on functions of a Sobolev space in our experiments is as good as optimal quadrature schemes (see Section 6); depending on the orthogonal polynomials, we get various quadrature rules, such as Gauss-Legendre quadrature for the Lebesgue measure on $$.

Quasi Monte-carlo methods employ a sequence of points with low discrepancy with uniform weights (Morokoff and Caflisch, 1994), leading to approximation errors of $O(1/n)$ for univariate functions with bounded variation, but typically with no adaptation to smoother functions.

Higher-dimensional integrals.

All of the methods above may be generalized for products of intervals $^{d}$, typically with $d$ small. For larger problems, Bayes-Hermite quadrature (O’Hagan, 1991) is essentially equivalent to the quadrature rules we study in this paper.

Some of the quadrature rules are constrained to have positive weights with unit sum (so that the positivity properties of integrals are preserved and constants are exactly inegrated). The quadrature rules we present do not satisfy these constraints. If these constraints are required, kernel herding (Chen et al., 2010; Bach et al., 2012) provides a novel way to select a sequence of points based on the conditional gradient algorithm, but with currently no convergence guarantees improving over $O(1/\sqrt{n})$ for infinite-dimensional spaces.

Theoretical results.

The best possible error for a quadrature rule with $n$ points has been well-studied in several settings; see Novak (1988) for a comprehensive review. For example, for ${\mathcal{X}}=$ and the space of Sobolev functions, which are RKHSs with eigenvalues of their integral operator decreasing as $m^{-2s}$, Novak (1988, Prop. 2 and 3, page 38) shows that the best possible quadrature rule for the uniform distribution and $g=1$ leads to an error rate of $n^{-s}$, as well as for any squared-integrable function $g$. The proof of these results (both upper and lower bounds) relies on detailed properties of Sobolev spaces. In this paper, we recover these results using only the decay of eigenvalues of the associated integral operator $\Sigma$, thus allowing straightforward extensions to many situations, like Sobolev spaces on manifolds such as hyperspheres (Hesse, 2006), where we also recover existing results (up to logarithmic terms).

Moreover, Novak (1988, page 17) shows that adaptive quadrature rules where points are selected sequentially with the knowledge of the function values at previous points cannot improve the worst-case guarantees. Our results do not recover this lower bound result for adaptivity.

Finally, Langberg and Schulman (2010) consider multiplicative errors in computing integrals and mainly focuses on different function spaces, such as ones used in clustering functionals. Although sampling quadrature points from a well-chosen density is common in the two approaches, the analysis tools are different. It would be interesting to see if some of these tools can be transferred to our RKHS setting.

From quadrature to function approximation and optimization.

The problem of quadrature, uniformly over all functions $g\in L_{2}(d\rho)$ that define the integral, is in fact equivalent to the full approximation of a function $h$ given values at $n$ points, where the approximation error is characterized in $L_{2}$-norm. Indeed, given the observations $h(x_{i})$, $i=1,\dots,n$, we build $\sum_{i=1}^{n}\alpha_{i}h(x_{i})$ as an approximation of $\int_{{\mathcal{X}}}g(x)h(x)d\rho(x)$. It turns out that the coefficients $\alpha_{i}$ are linear in $g$, that is, there exists $a_{i}\in L_{2}(d\rho)$ such that $\alpha_{i}=\langle a_{i},g\rangle_{L_{2}(d\rho)}$. This implies that $\sum_{i=1}^{n}h(x_{i})\langle a_{i},g\rangle_{L_{2}(d\rho)}$ is an approximation of $\langle h,g\rangle_{L_{2}(d\rho)}$. Thus, the worst case error with respect to $g$ in the unit ball of $L_{2}(d\rho)$ is \big{\|}\sum_{i=1}^{n}h(x_{i})a_{i}-h\big{\|}_{L_{2}(d\rho)}, that is, we have an approximation result of $h$ through observations of its values at certain points.

Novak (1988) considers the approximation problem in $L_{\infty}$-norm and shows that for Sobolev spaces, going from $L_{2}$- to $L_{\infty}$-norms incurs a loss of performance of $\sqrt{n}$. We recover partially these results in Section 5 from a more general perspective. When optimizing the points at which the function is evaluated (adaptively or not), the approximation problem is often referred to as experimental design (Cochran and Cox, 1957; Chaloner and Verdinelli, 1995).

Finally, a third problem is of interest (and outside of the scope of this paper), namely the problem of finding the minimum of a function given (potentially noisy) function evaluations. For noiseless problems, Novak (1988, page 26) shows that the approximation and optimization problems have the same worst-case guarantees (with no influence of adaptivity); this optimization problem has also been studied in the bandit setting (Srinivas et al., 2012) and in the framework of “Bayesian optimization” (see, e.g. Bull, 2011).

Theoretical Analysis

In this section, we provide approximation bounds for the random feature problem outlined in Section 2.4 (and thus the quadrature problem in Section 3). In Section 4.1, we provide generic upper bounds, which depend on the eigenvalues of the integral operator $\Sigma$ and present matching lower bounds (up to logarithmic terms) in Section 4.3. The upper-bound depends on specific distributions of samples that we discuss in Section 4.2. We then consider consequences of these results on quadrature (Section 4.4) and random feature expansions (Section 4.5).

The following proposition (see proof in Appendix B.1) determines the minimal number of samples required for a given approximation accuracy:

For $\lambda>0$ and a distribution with positive density $q$ with respect to $d\tau$, we consider

Let $v_{1},\dots,v_{n}$ be sampled i.i.d. from the density $q$, then for any $\delta\in(0,1)$, if

with probability greater than $1-\delta$, we have $\frac{1}{n}\sum_{i=1}^{n}q(v_{i})^{-1}\|\varphi(v_{i},\cdot)\|_{L_{2}(d\rho)}^{2}\leqslant\frac{2\mathop{\rm tr}\Sigma}{\delta}$ and

We can interpret the proposition above as follows: given any squared error $4\lambda>0$ and a distribution with density $q$, the number $n$ of samples from $q$ needed so that the unit ball of ${\mathcal{F}}$ is approximated by the ball of radius $2$ of $\hat{{\mathcal{F}}}$ is, up to logarithmic terms, at most a constant times $d_{\max}(q,\lambda)$, defined in Eq. (9). The result above is a statement for a fixed $q$ and $\lambda$ and this number of samples $n$ depends on these.

We could also invert the relationship between $\lambda$ and $n$, that is, answer the following question: given a fixed number $n$ of samples, what is the approximation error $\lambda$? This requires inverting the function $\lambda\mapsto d_{\max}(q,\lambda)$. This will be done in Section 4.2 for a specific distribution $q$ where the expression simplifies, together with specific examples from Section 2.3.

Finally, note that we also have a bound on $\frac{1}{n}\sum_{i=1}^{n}q(v_{i})^{-1}\|\varphi(v_{i},\cdot)\|_{L_{2}(d\rho)}^{2}$, which shows that our random functions are not too large on average (this constraint will be needed in the lower bound as well in Section 4.3).

It turns out that the final bound on the squared error is exactly proportional to the regularization parameter $\lambda$. As shown in Appendix B.1, this leads to an approximation $\hat{f}$ which is a linear function of $f$, as $\hat{f}=(\hat{\Sigma}+\lambda I)^{-1}\hat{\Sigma}f$, where $\hat{\Sigma}$ is a properly defined empirical integral operator and $\lambda>0$ is the regularization parameter. Then, Bernstein concentration inequalities for operators (Minsker, 2011) can be used in a way similar to the work of Bach (2013); El Alaoui and Mahoney (2014) on column sampling, to provide a bound on all desired quantities.

Result in expectation.

In Section 4.5, we will need a result in expectation. As shown at the end of Appendix B.1, as soons as, $\lambda\leqslant(\mathop{\rm tr}\Sigma)/4$ and $\displaystyle n\geqslant 5d_{\max}(\lambda)\log\frac{2(\mathop{\rm tr}\Sigma)d_{\max}(\lambda)}{\lambda},$ then

2 Optimized distribution

We may now consider a specific distribution that depends on the kernel and on $\lambda$, namely

for which $d_{\max}(q^{\ast}_{\lambda},\lambda)=d(\lambda)=\mathop{\rm tr}\Sigma(\Sigma+\lambda I)^{-1}$. With this distribution, we thus need to have $n\geqslant 5d(\lambda)\log\frac{16d(\lambda)}{\delta}$ with $d(\lambda)=\mathop{\rm tr}\Sigma(\Sigma+\lambda I)^{-1}$ is the degrees of freedom, a traditional quantity in the analysis of least-squares regression (Hastie and Tibshirani, 1990; Caponnetto and De Vito, 2007), which is always smaller than $d_{\max}(1,\lambda)$ and can be upper-bounded explicitly for many examples, as we now explain. The computation of $d_{\max}(1,\lambda)$ in the operator setting (for which we may use $q=1$), a quantity often referred to as the maximal leverage score (Mahoney, 2011), remains an open problem.

The quantity $d(\lambda)$ only depends on the integral operator $\Sigma$, that is, for all possible choices of square roots, i.e., all possible choices of feature expansions, the number of samples that our results guarantee is the same. This being said, some expansions may be more computationally practical than others, and when using the distribution with $q(v)=1$, the bounds will be different.

Given the singular value decomposition of $\varphi$ in Eq. (3), we have, for any $v\in{\mathcal{V}}$, $\varphi(v,\cdot)=\sum_{m\geqslant 1}\mu_{m}^{1/2}f_{m}(v)e_{m}$ and thus

which provides an explicit expression for the density $q_{\lambda}^{\ast}$.

For a given squared error value $\lambda$, the optimized distribution $q^{\ast}_{\lambda}$, while leading to the degrees of freedom that will happen to be optimal in terms of approximation, has two main drawbacks:

Dependence on $\lambda$: this implies that if we want a reduced error (i.e., a smaller $\lambda$), then the samples obtained from a higher $\lambda$, may not be reused to provably obtain the desired bound; in other words, the sampling is not anytime. For specific examples, e.g., quadrature with periodic kernels on $$with the uniform distribution, then$q=1$happens to be optimal for all$\lambda$, and thus, we may reuse samples for different values of the error.

Hard to compute in practice: the optimal distribution depends on a leverage score $\langle\varphi(v,\cdot),(\Sigma+\lambda I)^{-1}\varphi(v,\cdot)\rangle_{L_{2}(d\rho)}$, which may be hard to use for several reasons; first, it requires access to the infinite-dimensional operator $\Sigma$, which may be difficult; moreover, even if it possible to invert $\Sigma+\lambda I$, the set ${\mathcal{V}}$ might be particularly large and impractical to sample from. At the end of Section 4.1, we propose a simple algorithm based on sampling.

Eigenvalues and degrees of freedom.

In order to relate more directly to the eigenvalues of $\Sigma$, we notice that we may lower bound the degrees of freedom by a constant times the number $m^{\ast}(\lambda)$ of eigenvalues greater than $\lambda$:

We now make the assumption that there exists a $\gamma>0$ independent of $j$ such that

This assumption essentially states that the eigenvalues decay sufficiently homogeneously and is satisfied by $\mu_{m}\propto m^{-2\alpha}$ with $\gamma=(2\alpha-1)^{-1}$, $\mu_{m}\propto r^{m}$ with $\gamma=(1-r)^{-1}$ and similar bounds also hold for all examples in Section 2.3. It allows us to relate the degrees of freedom directly to eigenvalue decays.

Indeed, this implies that $\frac{1}{\lambda}\sum_{\mu_{m}<\lambda}\mu_{m}\leqslant\gamma{\rm max}(\{m,\ \mu_{m}\geqslant\lambda\})=m^{\ast}(\lambda)$ for all $\lambda\leqslant\mu_{1}$ (the largest eigenvalue) and thus

We can now restate the approximation result of Prop. 1 from Section 4.1 with the optimized distribution (see proof in Appendix B.2):

For $\lambda>0$ and the distribution with density $q^{\ast}_{\lambda}$ defined in Eq. (10) with respect to $d\tau$, with degrees of freedom $d(\lambda)$. Let $v_{1},\dots,v_{n}$ be sampled i.i.d. from the density $q$, defining the kernel (and its associated RKHS $\hat{{\mathcal{F}}}$) $\hat{k}(x,y)=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{q(v_{i})}\varphi({v_{i}},x)\varphi({v_{i}},y)$. Then, for any $\delta\in(0,1)$, with probability $1-\delta$, we have:

if \displaystyle n\geqslant 5\ d(\lambda)\log\big{[}16d(\lambda)/\delta\big{]},

if Eq. (11) is satisfied, and, by choosing $m\leqslant\frac{n}{5(1+\gamma)\log\frac{16n}{5\delta}}$, and $\lambda=\mu_{m}$.

The statement (a) above, is a simple corollary of Prop. 1, and goes from level of error $\lambda$ to minimum number $n$ of samples. The statement (b) goes in the other direction, that is, from the number of samples $n$ to the achieved approximation error. It depends on the eigenvalues $\mu_{m}$ of the integral operator taken at $m=O(n/\log(n))$. For example, for polynomial decays of eigenvalues of the form $\mu_{m}=O(m^{-2s})$, we get (non squared) errors proportional to $(\log n)^{s}n^{-s}$ for $n$ samples, while for geometric decays, we get geometric errors as a function of the number $n$ of samples.

Note however that for the statement (b) to hold, we need to sample the points $v_{1},\dots,v_{n}$ from the distribution $q_{\mu_{m}}^{\ast}$, that is, for different numbers of samples $n$, the distribution is unfortunately different (except in special cases). It would be interesting to study the properties of independent but not identically distributed samples $v_{1},\dots,v_{n}$ and the possibility of achieving the same rate adaptively.

Corollary for Sobolev spaces.

The same approximation results holds for translation-invariant kernels on $^{d}$; but when $d\rho$ is the uniform distribution, as shown in Section 4.4, the optimized distribution for the quadrature case is still the uniform distribution, for all values of $\lambda$, and can thus be computed.

Algorithm to estimate the optimized distribution.

This implies that the density of the optimized distribution with respect to the uniform measure on $\{v_{1},\dots,v_{N}\}$ is proportional to \big{(}T^{\ast}T(T^{\ast}T+\lambda I)^{-1}\big{)}_{ii}. We can then sample any number $n$ of points from resampling from $\{v_{1},\dots,v_{N}\}$ from the density above. The computational complexity is $O(N^{3})$. A detailed analysis of the approximation properties of this algorithm is outside the scope of this paper.

We have $(T^{\ast}T)_{ij}={\eta_{i}^{1/2}\eta_{j}^{1/2}}\int_{\mathcal{X}}\varphi(v_{i},x)\varphi(v_{j},x)d\rho(x)$. In some cases, it can be computed in closed form—such as for quadrature where this is equal to ${\eta_{i}^{1/2}\eta_{j}^{1/2}}k(v_{i},v_{j})$. In some others, it requires i.i.d. samples $x_{1},\dots,x_{M}$ from $d\rho$, and the estimate: ${{\eta_{i}^{1/2}\eta_{j}^{1/2}}}{M^{-1}}\sum_{k=1}^{M}\varphi(v_{i},x_{k})\varphi(v_{j},x_{k})$.

3 Lower bound

In this section, we aim at providing lower-bounds on the number of samples required for a given accuracy. We have the following result (see proof in Appendix B.3):

For $\delta\in(0,1)$, if we have a family $\psi_{1},\dots,\psi_{n}\in L_{2}(d\rho)$ such that

then $\displaystyle n\geqslant\frac{{\rm max}\{m,\ \mu_{m}\geqslant 144\lambda\}}{4\log\frac{10\mathop{\rm tr}\Sigma}{\lambda\delta}}$.

The proof technique not surprisingly borrows tools from minimax estimation over ellipsoids, namely the Varshamov-Gilbert’s lemma.

We obtain matching upper and lower bounds up to logarithmic terms, using only the decay of eigenvalues $(\mu_{m})_{m\geqslant 1}$ of the integral operator $\Sigma$ (of course, if sampling from the optimized distribution $q^{\ast}_{\lambda}$ is possible). Indeed in that case, as shown in Prop. 2, we have shown that we need at most 10\ d(\lambda)\log\big{[}2d(\lambda)\big{]}, where $d(\lambda)$ is the degrees of freedom, which is upper and lower bounded by a constant times $m^{\ast}(\lambda)={\rm max}\{m,\ \mu_{m}\geqslant\lambda\}$.

In order to obtain such a bound, we need to constrain both $\|\beta\|_{2}$ and the norms of the vectors $\psi_{i}$, which correspond to bounded features for the random feature interpretation and tolerance to noise for the quadrature interpretation. We choose our scaling to match the constraints we have in Prop. 1, for which the parameter $\delta$ ends up entering the lower bound logarithmically.

4 Quadrature

We may specialize the results above to the quadrature case, namely give a formulation where the features $\varphi$ do not appear (or equivalently using $\psi$ defined in Section 3.2). This is a special case where ${\mathcal{V}}={\mathcal{X}}$ and $\varphi=\psi$. In terms of operators $T$ in Section 2.2, this corresponds to $T=\Sigma^{1/2}$.

Following Section 4.1, we have an expression for the optimized distribution, both in terms of operators, as follows,

and in terms of eigenvalues and eigenvectors of $k$, that is,

While this is uniform in some special cases (uniform distribution on $$ and Sobolev kernels, as shown below), this is typically hard to compute and sample from. An algorithm for approximating it was presented at the end of Section 4.1.

A weakness of our result is that in general our optimized distribution $q^{\ast}_{\lambda}(x)$ depends on $\lambda$ and thus on the number of samples. In some cases with symmetries (i.e., uniform distribution on $$or the hypersphere),$q^{\ast}_{\lambda}$happens to be constant for all$\lambda$. Note also that we have observed empirically that in some cases,$q^{\ast}_{\lambda}$converges to a certain distribution when$\lambda$ tends to zero (see an example in Section 6).

Sobolev spaces.

For the special case of Sobolev spaces on $^{d}$ with $d\rho$ the uniform distribution, the optimized distribution in Eq. (12) is also the uniform distribution. Indeed, the eigenfunctions of the integral operator $\Sigma$ are $d$-th order tensor products of the uni-dimensional Fourier basis (the constant and all pairs of sine/cosine at a given frequency), with the same eigenvalue for the $2^{d}$ possibilities of sines/cosines for a given multi-dimensional frequency $(m_{1},\dots,m_{d})$. Therefore, when summing all squared values of the eigenfunctions corresponding to $(m_{1},\dots,m_{d})$, we end up with the sum $\sum_{a\in\{0,1\}^{d}}\prod_{i=1}^{d}\cos^{2a_{i}}(2\pi m_{i}x_{i})\sin^{2(1-a_{i})}(2\pi m_{i}x_{i})$, which ends up being constant equal to one (and thus independent of $x$) because $\cos^{2a_{i}}(2\pi m_{i}x_{i})+\sin^{2a_{i}}(2\pi m_{i}x_{i})=1$.

Quadrature rule.

We assume that points $x_{1},\dots,x_{n}$ are sampled from the distribution with density $q$ with respect to $d\rho$. The quadrature rule for a function $h\in{\mathcal{F}}$ is $\sum_{i=1}^{n}\frac{\beta_{i}h(x_{i})}{q(x_{i})^{1/2}}$. To compute $\beta$, we need to minimize with respect to $\beta$ the error: