Kernel quadrature with DPPs

Introduction

Numerical integration is an important tool for Bayesian methods and model-based machine learning . Formally, numerical integration consists in approximating

In this paper, we propose a new quadrature rule for functions in a given RKHS. Our nearest scientific neighbour is , but instead of sampling nodes independently, we leverage dependence and use a repulsive distribution called a projection determinantal point process (DPP), while the weights are obtained through a simple quadratic optimization problem. DPPs were originally introduced by as probabilistic models for beams of fermions in quantum optics. Since then, DPPs have been thoroughly studied in random matrix theory , and have more recently been adopted in machine learning and Monte Carlo methods .

Related work on kernel-based quadrature

When the integrand $f$ belongs to the RKHS $\mathcal{F}$ of kernel $k$ , the quadrature error reads

Bayesian Quadrature initially considered a fixed set of nodes and put a Gaussian process prior on the integrand $f$. Then, the weights were chosen to minimize the posterior variance of the integral of $f$. If the kernel of the Gaussian process is chosen to be $k$, this amounts to minimizing the RHS of (2). The case of the Gaussian reference measure was later investigated in detail , while parametric integrands were considered in . Rates of convergence were provided in for specific kernels on compact spaces, under a fill-in condition that encapsulates that the nodes must progressively fill up the (compact) space.

2 Leverage-score quadrature

In , the author proposed to sample the nodes $(x_{j})$ i.i.d. from some proposal distribution $q$, and then pick weights $\hat{\bm{w}}$ in (1) that solve the optimization problem

for some regularization parameter $\lambda>0$. Proposition 1 gives a bound on the resulting approximation error of the mean element for a specific choice of proposal pdf, namely the leverage scores

Let $\delta\in$, and $d_{\lambda}=\operatorname{Tr}\bm{\Sigma}(\bm{\Sigma}+\lambda\bm{I})^{-1}$. Assume that $N\geq 5d_{\lambda}\log(16d_{\lambda}/\delta)$, then

Projection determinantal point processes

not to be mistaken for the RKHS kernel $k$. One can show [25, Lemma 21] that

The probability of co-occurrence is thus always smaller than that of a Poisson process with the same intensity. In this sense, a projection DPP with symmetric kernel is a repulsive distribution, and $\operatorname*{\mathfrak{K}}$ encodes its repulsiveness.

One advantage of DPPs is that they can be sampled exactly. Because of the orthonormality of $(\psi_{n})$, one can write the chain rule for (9); see . Sampling each conditional in turn, using e.g. rejection sampling , then yields an exact sampling algorithm. Rejection sampling aside, the cost of this algorithm is cubic in $N$ without further assumptions on the kernel. Simplifying assumptions can take many forms. In particular, when $d=1$, and $\omega$ is a Gaussian, gamma , or beta pdf, and $(\psi_{n})$ are the orthonormal polynomials with respect to $\omega$, the corresponding DPP can be sampled by tridiagonalizing a matrix with independent entries, which takes the cost to $\mathcal{O}(N^{2})$ and bypasses the need for rejection sampling. For further information on DPPs see .

Kernel quadrature with projection DPPs

where we recall that $(e_{n})$ are the normalized eigenfunctions of the integral operator $\bm{\Sigma}$. The weights $\bm{w}$ are obtained by solving the optimization problem

is the reconstruction operatorThe reconstruction operator $\bm{\Phi}$ depends on the nodes $x_{j}$, although our notation doesn’t reflect it for simplicity.. In Section 4.1 we prove that (11) almost surely has a unique solution $\hat{\bm{w}}$ and state our main result, an upper bound on the expected approximation error $\|\mu_{g}-\bm{\Phi}\hat{\bm{w}}\|^{2}_{\mathcal{F}}$ under the proposed Projection DPP. Section 4.2 gives a sketch of the proof of this bound.

Assuming that nodes $(x_{j})_{j\in[N]}$ are known, we first need to solve the optimization problem (11) that relates to problem (5) without regularization ($\lambda=0$). Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$, then

where $\bm{K}(\bm{x})=(k(x_{i},x_{j}))_{i,j\in[N]}$. The right-hand side of (13) is quadratic in $\bm{w}$, so that the optimization problem (11) admits a unique solution $\hat{\bm{w}}$ if and only if $\bm{K}(\bm{x})$ is invertible. In this case, the solution is given by $\hat{\bm{w}}=\bm{K}(\bm{x})^{-1}\mu_{g}(x_{j})_{j\in[N]}$. A sufficient condition for the invertibility of $\bm{K}(\bm{x})$ is given in the following proposition.

Assume that the matrix $\bm{E}(\bm{x})=(e_{i}(x_{j}))_{i,j\in[N]}$ is invertible, then $\bm{K}(\bm{x})$ is invertible.

The proof of Proposition 2 is given in Appendix D.1. Since the pdf (9) of the projection DPP with kernel (10) is proportional to $\operatorname{Det}^{2}\bm{E}(\bm{x})$, the following corollary immediately follows.

We now give our main result that uses nodes $(x_{j})_{j\in[N]}$ drawn from a well-chosen projection DPP.

This compares favourably with herding, for instance, which comes with a rate in $\mathcal{O}(\frac{1}{N})$ for the quadrature based on herding with uniform weights .

2 Bounding the approximation error under the DPP

In this section, we give the skeleton of the proof of Theorem 1, referring to the appendices for technical details. The proof is in two steps. First, we give an upper bound for the approximation error $\|\mu_{g}-\bm{\Phi}\hat{\bm{w}}\|^{2}_{\mathcal{F}}$ that involves the maximal principal angle between the functional subspaces of $\mathcal{F}$

DPPs allow closed form expressions for the expectation of trigonometric functions of such angles; see and Appendix E.1 for the geometric intuition behind the proof. The second step thus consists in developing the expectation of the bound under the DPP.

Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$ be such that $\operatorname{Det}\bm{E}(\bm{x})\neq 0$. By Proposition 2, $\bm{K}(\bm{x})$ is non singular and $\dim\mathcal{T}(\bm{x})=N$. The optimal approximation error writes

In other words, (16) equates the approximation error to $\|\bm{\Pi}_{\mathcal{T}(\bm{x})^{\perp}}\mu_{g}\|^{2}_{\mathcal{F}}$, where $\bm{\Pi}_{\mathcal{T}(\bm{x})^{\perp}}$ is the orthogonal projection onto $\mathcal{T}(\bm{x})^{\perp}$. Now we have the following lemma.

Similarly, we can define the $N$ principal angles $\theta_{n}(\mathcal{T}(\bm{x}),\mathcal{E}^{\mathcal{F}}_{N})\in\left[0,\frac{\pi}{2}\right]$ for $n\in[N]$ between the subspaces $\mathcal{E}^{\mathcal{F}}_{N}$ and $\mathcal{T}(\bm{x})$. These angles quantify the relative position of the two subspaces. See Appendix C.3 for more details about principal angles. Now, we have the following lemma.

Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$ such that $\operatorname{Det}\bm{E}(\bm{x})\neq 0$. Then

To sum up, we have so far bounded the approximation error by the geometric quantity in the right-hand side of (19). Where projection DPPs shine is in taking expectations of such geometric quantities.

2.2 Taking the expectation under the DPP

The analysis in Section 4.2.1 is valid whenever $\operatorname{Det}\bm{E}(\bm{x})\neq 0$. As seen in Corollary 1, this condition is satisfied almost surely when $\bm{x}$ is drawn from the projection DPP of Theorem 1. Furthermore, the expectation of the right-hand side of (19) can be written in terms of the eigenvalues of the kernel $k$.

3 Discussion

In comparison with , we emphasize that the dependence of our bound on the eigenvalues of the kernel $k$, via $r_{N}$, is explicit. This is in contrast with Proposition 1 that depends on the eigenvalues of $\bm{\Sigma}$ through the degree of freedom $d_{\lambda}$ so that the necessary number of samples $N$ diverges when $\lambda\rightarrow 0$. On the contrary, our quadrature requires a finite number of points for $\lambda=0$. It would be interesting to extend the analysis of our quadrature in the regime $\lambda>0$.

Numerical simulations

so that $\mathcal{F}=\mathcal{F}_{s}$ is the Sobolev space of order $s$ on $$. Note that$k_{s}$can be expressed in closed form using Bernoulli polynomials . We take$g\equiv 1$in (1), so that the mean element$\mu_{g}\equiv 1$. We compare the following algorithms:$(i)$the quadrature rule DPPKQ we propose in Theorem 1,$(ii)$the quadrature rule DPPUQ based on the same projection DPP but with uniform weights, implicitly studied in ,$(iii)$the kernel quadrature rule (5) of , which we denote LVSQ for leverage score quadrature, with regularization parameter$\lambda\in\{0,0.1,0.2\}$(note that the optimal proposal is$q_{\lambda}^{*}\equiv 1$),$(iv)$herding with uniform weights ,$(v)$sequential Bayesian quadrature (SBQ) with regularization to avoid numerical instability, and$(vi)$Bayesian quadrature on the uniform grid (UGBQ). We take$N\in$. Figures 1(a) and 1(b) show log-log plots of the worst case quadrature error w.r.t.$N$, averaged over 50 samples for each point, for$s\in\{1,3\}$.

We observe that the approximation errors of all first four quadratures converge to with different rates. Both UGBQ and DPPKQ converge to with a rate of $\mathcal{O}(N^{-2s})$, which indicates that our $\mathcal{O}(N^{2-2s})$ bound in Theorem 1 is not tight in the Sobolev case. Meanwhile, the rate of DPPUQ is $\mathcal{O}(N^{-2})$ across the three values of $s$: it does not adapt to the regularity of the integrands. This corresponds to the CLT proven in . LVSQ without regularization converges to slightly slower than $\mathcal{O}(N^{-2s})$. Augmenting $\lambda$ further slows down convergence. Herding converges at an empirical rate of $\mathcal{O}(N^{-2})$, which is faster than the rate $\mathcal{O}(N^{-1})$ predicted by the theoretical analysis in . SBQ is the only one that seems to plateau for $s=3$, although it consistently has the best performance for low $N$. Overall, in the Sobolev case, DPPKQ and UGBQ have the best convergence rate. UGBQ – known to be optimal in this case – has a better constant.

Now, for a multidimensional example, consider the “Korobov" kernel $k_{s}$ defined on $^{d}$ by

We still take $g\equiv 1$ in (1) so that $\mu_{g}\equiv 1$. We compare $(i)$ our DPPKQ, $(ii)$ LVSQ without regularization ($\lambda=0$), $(iii)$ the kernel quadrature based on the uniform grid UGBQ, $(iv)$ the kernel quadrature SGBQ based on the sparse grid from , $(v)$ the kernel quadrature based on the Halton sequence HaltonBQ . We take $N\in$ and $s=1$. The results are shown in Figure 1(c). This time, UGBQ suffers from the dimension with a rate in $\mathcal{O}(N^{-2s/d})$, while DPPKQ, HaltonBQ and LVSQ $(\lambda=0)$ all perform similarly well. They scale as $\mathcal{O}((\log N)^{2s(d-1)}N^{-2s})$, which is a tight upper bound on $\sigma_{N+1}$, see and Appendix B. SGBQ seems to lag slightly behind with a rate $\mathcal{O}((\log N)^{2(s+1)(d-1)}N^{-2s})$ .

2 The Gaussian kernel

Conclusion

In this article, we proposed a quadrature rule for functions living in a RKHS. The nodes are drawn from a DPP tailored to the RKHS kernel, while the weights are the solution to a tractable, non-regularized optimization problem. We proved that the expected value of the squared worst case error is bounded by a quantity that depends on the eigenvalues of the integral operator associated to the RKHS kernel, thus preserving the natural feel and the generality of the bounds for kernel quadrature . Key intermediate quantities further have clear geometric interpretations in the ambient RKHS. Experimental comparisons suggest that DPP quadrature favourably compares with existing kernel-based quadratures. In specific cases where an optimal quadrature is known, such as the uniform grid for 1D periodic Sobolev spaces, DPPKQ seems to have the optimal convergence rate. However, our generic error bound does not reflect this optimality in the Sobolev case, and must thus be sharpened.

We have discussed room for improvement in our proofs. Further work should also address exact sampling algorithms, which do not exist yet when the spectral decomposition of the integral operator is not known. Approximate algorithms would also suffice, as long as the error bound is preserved.

We acknowledge support from ANR grant BoB (ANR-16-CE23-0003) and région Hauts-de-France. We also thank Adrien Hardy and the reviewers for their detailed and insightful comments.

References

Appendix A Implementation details

In this section, we give details on the repulsion kernels in each example of the main paper, and explain how we sampled from the corresponding DPPs. In short, we relied on matrix models for univariate cases, and vanilla DPP sampling for multivariate settings.

A.2 The one-dimensional Gaussian kernel

For notational convenience, we further let

Now, the Mercer decomposition of $k_{\gamma}$ reads

and $H_{m}$ is the $m$-th Hermite polynomial (i.e., orthonormal polynomials for the pdf of a unit Gaussian). Now, denote

A.3 The case of a tensor product of RKHSs

We consider the case where $\mathcal{F}$ writes as a tensor product of RKHSs, with the associated kernel

A.3.2 Fixing an order on multi-indices

The definition of the projection DPP and its kernel $\mathfrak{K}$ now require that we fix an order on multi-indices. We choose an order $\prec$ that keeps eigenvalues decreasing, as in the univariate case where $\sigma_{1}\geq\sigma_{2}\geq\dots$. Whenever the univariate eigenvalues take the form $\sigma_{i}=\frac{1}{(1+i)^{\eta}}$ with $\eta>0$, such as in the Korobov case, it holds

Now, if the eigenvalues takes the form $\displaystyle\sigma_{i}={\eta^{-i}}$, with $\eta>1$, as in the Gaussian case,

In the multivariate Korobov and the Gaussian cases, we thus define in this work $\bm{u}\prec\bm{v}$ as (44) or (47), respectively.

We sampled from the corresponding DPP using the generic sampling algorithm in , using the uniform and Gaussian distributions as proposal in the successive rejection sampling steps for the Korobov and Gaussian cases, respectively.

Appendix B Supplementary simulations

For the Gaussian RKHS in dimension $d$, it holds

We consider the case of Korobov spaces with $d\in\{2,3\}$ and $s\in\{1,2\}$ and compare the quadrature error of the same algorithms as in 5.1. The results are compiled in Figure 3. The numerical simulations confirm the dependencies of the theoretical bounds of the different algorithms to the dimension $d$ and the regularity $s$. In particular, UGBQ have better performance for high values of $s$ and low values of $N$ while its asymptotic behaviour is still the same $\mathcal{O}(N^{-2s/d})$. Moreover, the empirical rate of SGBQ is similar to its theoretical rate $\mathcal{O}((\log N)^{2(s+1)(d-1)}N^{-2s})$ . Finally, the rate $\mathcal{O}((\log N)^{2s(d-1)}N^{-2s})$ is confirmed also for the algorithms DPPKQ, LVSQ $(\lambda=0)$ and HaltonBQ.

B.2 The multi Gaussian ensemble

We consider the case of Gaussian spaces with $d\in\{2,3\}$. The kernel $k_{\gamma,d}$ and the reference measure are the tensor product of respectively the same kernel and the same measure used in Section 5.2. We compare DPPKQ and Bayesian quadrature based on the tensor product of Gauss-Hermite nodes noted GHBQ. Note that a variant of this algorithm was proposed in : the quadrature nodes are the tensor product of the Gauss-Hermite nodes however the weights were calculated differently. The authors proved under an assumption on the stability of the weights (that was verified empirically) that the rate of convergence is $\mathcal{O}(dr^{d}\beta^{\prime d}e^{-\delta^{\prime}dN^{1/d}})$, where $r$ is a constant that quantify the stability of the weights, and $\beta^{\prime},\delta^{\prime}$ are constants that depend simultaneously on the the stability of the weights and length scales of the kernel and the measure. The results are compiled in Figure 4.

The numerical simulations shows that the empirical rate of DPPKQ is $\mathcal{O}(e^{-\delta dN^{1/d}})$ that is slightly better than its theoretical rate $\mathcal{O}(e^{-\delta d!^{1/d}N^{1/d}})$. Moreover, we observe that the empirical rate of DPPKQ is better than the empirical rate of HGBQ.

Appendix C Mercer’s theorem, leverage scores, and principal angles

For the sake of completeness, this section gathers some known results, which will be used to prove our own. We will need a general version of Mercer’s theorem, as usual for kernel methods, see Section C.1. On a more technical ground, we will also need formulas for leverage score changes under rank 1 updates, see Section C.2. Finally, Section C.3 covers principal angles between subspaces of a Hilbert space, which bridge the gap between pairs of Hilbert subspaces and determinants, and facilitate taking expectations in Theorem 1.

where the convergence is absolute and uniform.

is compact. A sufficient condition for this assumption is the integrability of the diagonal (Lemma 2.3, ):

Note that this condition is not necessary (Example 2.9, ). Now, under the compact embedding assumption, the pointwise convergence of the Mercer decomposition to the kernel $k$ is equivalent to the injectivity of the embedding $I_{\mathcal{F}}$ (Theorem 3.1, ).

C.2 Leverage score changes under rank 1 updates

In this section we prove a lemma inspired from Lemma 5 in . This lemma concerns the changes of leverage scores under rank 1 updates.

while the cross-leverage score between the $i$-th column and the $j$-th column is defined by

The proof of this lemma is similar to Lemma 5 in . We recall the proof for completeness.

(Adapted from ) The Sherman-Morrison formula applied to $\bm{A}\bm{W}\bm{W}^{\operatorname{\intercal}}\bm{A}^{\operatorname{\intercal}}$ and the vector $\sqrt{\rho}\bm{a}_{i}$ yields

By definition of $\tau_{i}(\bm{A}\bm{W})$

Now let $j\in[M]-\{i\}$. By definition of $\tau_{j}(\bm{A}\bm{W})$

C.3 Principal angles between subspaces in Hilbert spaces

We recall in this section the definition of principal angles between subspaces in Hilbert spaces and connect them to the determinant of the Gramian matrix of their orthonormal bases.

Let $\mathcal{H}$ be a Hilbert space. Let $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ be two finite-dimensional subspaces of $\mathcal{H}$ with $N=\dim\mathcal{P}_{1}=\dim\mathcal{P}_{2}$. Denote $\bm{\Pi}_{\mathcal{P}_{1}}$ and $\bm{\Pi}_{\mathcal{P}_{2}}$ the orthogonal projections of $\mathcal{H}$ onto these two subspaces. There exist two orthonormal bases for $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ denoted $(\bm{v}_{i}^{1})_{i\in[N]}$ and $(\bm{v}_{i}^{2})_{i\in[N]}$, and a set of angles $\theta_{i}(\mathcal{P}_{1},\mathcal{P}_{2})\in[0,\frac{\pi}{2}]$ such that

We refer to for the proof in the finite-dimensional case and for the general case. The following result shows that the principal angles are somewhat independent of the choice of orthonormal bases. It can be found in for the finite dimensional case. We give here the proof for the general case, for the sake of completeness.

Let $(\bm{w}^{1}_{i})_{i\in[N]}$ be any orthonormal basis of $\mathcal{P}_{1}$ and $(\bm{w}^{2}_{i})_{i\in[N]}$ be any orthonormal basis of $\mathcal{P}_{2}$, and let $\bm{W}=(\langle\bm{w}^{1}_{i},\bm{w}^{2}_{j}\rangle_{\mathcal{H}})_{1\leq i,j\leq N}$ and $\bm{G}=\bm{W}\bm{W}^{\operatorname{\intercal}}$. Then the eigenvalues of $\bm{G}$ are the $\cos^{2}\theta_{i}(\mathcal{P}_{1},\mathcal{P}_{2})$. In particular, $\operatorname{Det}^{2}\bm{W}=\operatorname{Det}\bm{G}=\prod\limits_{i\in[N]}\cos^{2}\theta_{i}(\mathcal{P}_{1},\mathcal{P}_{2})$.

where $\displaystyle\bm{V}=(\langle\bm{v}^{1}_{i},\bm{v}^{2}_{j}\rangle_{\mathcal{H}})_{1\leq i,j\leq N}$. Then

Thus the eigenvalues of $\bm{G}$ are the eigenvalues of $\bm{V}\bm{V}^{\operatorname{\intercal}}$. By Proposition 5, the diagonal elements of $\bm{V}$ are

We finish the proof by showing that the anti-diagonal elements satisfy

Finally, $\bm{V}$ is a diagonal matrix and the eigenvalues of $\bm{G}$ are the $\cos^{2}\theta_{i}(\mathcal{P}_{1},\mathcal{P}_{2})$. ∎

Appendix D Proofs of our results

The rest of Section D deals with Theorem 1, our upper bound on the approximation error of DPP-based kernel quadrature. The proof is rather long, but can be decomposed in four steps, which we now introduce for ease of reading.

First, we prove Lemma 1, which separates the search for an upper bound into examining the contribution of the three terms in (17); this is Section D.2. The first two terms of (17) only depend on the function $g$ in (1), and we leave them be. The third term is more geometric, and relates to the approximation error of the space spanned by $(e^{\cal F}_{n})_{n\in[N]}$ by the (random) subspace ${\cal T}(\bm{x})$.

Second, in Section D.3, we bound this geometric term for a fixed DPP realization $\bm{x}$. We pay attention to obtain a bound that will later yield a tractable expectation under that DPP. This is done in Proposition 4, which in turn requires two intermediate results, Lemma 4 and Proposition 6.

Third, we take the expectation of the bound in Proposition 4 under the proposed DPP. This is done in Proposition 3, which is proven thanks to Proposition 2, Lemmas 2, 5 & 6. This is Section D.4.

Fourth, Theorem 1 is obtained in Section D.5, using the results of the previous steps, and an argument to reduce the proof to RKHSs with flat initial spectrum.

Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$ such that $\operatorname{Det}\bm{E}(\bm{x})\neq 0$, and define

Thus to prove that $\operatorname{Det}\bm{K}(\bm{x})>0$, it is enough to prove that the $\operatorname{Det}\bm{K}_{M}(\bm{x})$ is larger than a positive real number for $M$ large enough. We write

with $\bm{F}_{M}(\bm{x})=(e_{i}(x_{j}))_{(i,j)\in[M]\times[N]}$ and $\Sigma_{M}$ is a diagonal matrix containing the first $M$ eigenvalues $(\sigma_{m})$. The Cauchy-Binet identity yields

so that $\bm{K}(\bm{x})$ is a.s. invertible. ∎

D.2 Proof of Lemma 1

Note that $\bm{\Sigma}^{1/2}=\bm{\Sigma}_{N}^{1/2}+\bm{\Sigma}_{N}^{\perp 1/2}$ and

Now, recall that the $(e_{n}^{\mathcal{F}})_{n\in[N]}$ is orthonormal. Moreover for $n\in[N]$, $e_{n}^{\mathcal{F}}$ is an eigenfunction of $\bm{\Sigma}_{N}^{1/2}$ and the corresponding eigenvalue is $\sqrt{\sigma}_{n}$. Thus

D.3 Proof of Proposition 4

Proposition 4 gives an upper bound to the term $\max\limits_{n\in[N]}\sigma_{n}\|\bm{\Pi}_{\mathcal{T}(\bm{x})^{\perp}}e_{n}^{\mathcal{F}}\|_{\mathcal{F}}^{2}$ that appears in Lemma 1. We first prove a technical result, Lemma 4, and then combine it with Proposition 6 to finish the proof. We conclude with the proof of Proposition 6.

Indeed, $\displaystyle\|\bm{\Pi}_{\mathcal{T}(\bm{x})^{\perp}}e_{n}^{\mathcal{F}}\|^{2}_{\mathcal{F}}=1-\|\bm{\Pi}_{\mathcal{T}(\bm{x})}e_{n}^{\mathcal{F}}\|^{2}_{\mathcal{F}}$ since $\|e_{n}^{\mathcal{F}}\|^{2}_{\mathcal{F}}=1$. Thus it is sufficient to prove that $\|\bm{\Pi}_{\mathcal{T}(\bm{x})}e_{n}^{\mathcal{F}}\|^{2}_{\mathcal{F}}=\Delta_{n}^{\cal F}(\bm{x})$. This boils down to showing that $\bm{K}(\bm{x})^{-1}$ is the matrix of the inner product $\langle\cdot,\cdot\rangle_{\cal F}$ restricted to ${\mathcal{T}(\bm{x})}$.

We give the proof of (108); the proof of (109) follows the same lines.

where the $c_{i}$ are the elements of the vector $\bm{c}=\bm{K}(\bm{x})^{-1}e_{n}^{\mathcal{F}}(\bm{x})$. Then

Combining (112) and (113) along with the definition of the vector $\bm{c}=\bm{K}(\bm{x})^{-1}e_{n}^{\mathcal{F}}(\bm{x})$ yields

D.3.2 End of the proof of Proposition 4

By Lemma 4, the inequality (22) in Proposition 4 is equivalent to

As an intermediate remark, note that in the special case $n=1$, by construction

where $\prec$ is the Loewner order, the partial order defined by the convex cone of positive semi-definite matrices. Thus

For $n\neq 1$, the proof is much more subtle. Indeed, a naive application of the inequality (117) would lead to the following inequality

Now we have the following useful proposition.

For $n\in[N]\smallsetminus\{1\}$, we have

For ease of reading, we first show that inequality (115) and therefore Proposition 4 is easily deduced from this Proposition 6 and then give its proof.

Then we use (125) that is connected to the rank-one update from the kernel $k^{(n-1)}$ to $k^{(n)}$ so that

Then we apply (126) to the r.h.s. again $N-n-1$ times to finally get:

D.3.3 Proof of Proposition 6

which prepares the use of Lemma 3 in Section C.2. By definition of the $n$-th leverage score of the matrix $\bm{A}$, see (54) in Section C.2,

where $\displaystyle\rho_{n}=\frac{\sigma_{1}}{\sigma_{n}}-1$. Thus

As above, we conclude the proof by considering the limit $M\to\infty$

This proves inequality (126) and concludes the proof of Proposition 6. ∎

D.4 Proof of Proposition 3

We first prove two lemmas that are necessary to prove Proposition 3.

Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$ such that $\operatorname{Det}^{2}\bm{E}(\bm{x})\neq 0$. Then,

The condition $\operatorname{Det}^{2}\bm{E}(\bm{x})\neq 0$ yields by Proposition 2 that $\bm{K}(\bm{x})$ is non singular. Thus $\dim\mathcal{T}(\bm{x})=N$. Let $(t_{i})_{i\in[N]}$ an orthonormal basis of $\mathcal{T}(\bm{x})$ with respect to $\langle.,.\rangle_{\mathcal{F}}$. Using Corollary 2, and the fact that $(e_{n}^{\mathcal{F}})_{n\in[N]}$ is an orthonormal basis of $\mathcal{E}^{\mathcal{F}}_{N}$ according to $\langle.,.\rangle_{\mathcal{F}}$,

Now, let $\bm{c}_{i}$ the columns of the matrix $\bm{C}(\bm{x})$. $(t_{i})_{i\in[N]}$ is an orthonormal basis of $\mathcal{T}(\bm{x})$ with respect to $\langle.,.\rangle_{\mathcal{F}}$, then by (147)

Combining (146), (152) and (155) concludes the proof of Lemma 5:

Let $\bm{x}=(x_{1},\dots,x_{N})\in\mathcal{X}^{N}$. From (79)

Then by monotone convergence theorem, $\displaystyle\bm{x}\mapsto\frac{1}{N!}\operatorname{Det}\bm{K}(\bm{x})$ is mesurable and

D.4.2 End of the proof of Proposition 3

Then by Lemma 5 and the fact that $\operatorname{Det}^{2}\bm{E}^{\mathcal{F}}(\bm{x})=\prod\limits_{n\in[N]}\sigma_{n}\operatorname{Det}^{2}\bm{E}(\bm{x})$

Then, taking the expectation with respect to $\bm{x}$ resulting from a DPP of kernel $\operatorname*{\mathfrak{K}}(x,y)$,

D.5 Proof of Theorem 1

As a consequence, by writing $r_{N}=\sum\limits_{m\geq N+1}\sigma_{m}$,

which can be plugged in Lemma 1 to conclude the proof. ∎

Appendix E The intuitions behind the algorithm

The algorithm presented in this article is based on several intuitions. In this section, we summarize these intuitions.

which has a tractable expectation under the projection DPP that we consider in this paper. As an illustration of (180), Figure 6 compares the quality of approximation of a mean element $\mu_{g}$ using kernel interpolation based on two configurations of nodes: the first configuration (top) is well spread and the second configuration (bottom) is not. Observe that the largest principal angle $\theta_{N}$ for the first configuration is around $\pi/4$, so that $\tan^{2}\theta_{N}\approx 1$; while it is around $\pi/2$ for the second configuration so that $\tan^{2}\theta_{N}\gg 1$. Now observe that the first design of nodes gives the best reconstruction. This observation is consistent with (180).

E.2 The inclusion probability of DPPs and the Christoffel functions

The optimal distribution $q_{\lambda}$, see section 2.2, can be linked to the so-called Christoffel functions . These functions are rooted in the literature on orthogonal polynomials . To make it simpler, we introduce them in dimension $d=1$. They are defined by

Christoffel functions have a more explicit form that can be used for pointwise evaluation

The authors derived an asymptotic equivalent of the function $C_{\lambda,w,k}$ in the regime $\lambda\rightarrow 0$ under some assumptions on the kernel. Furthermore, they proved that $C_{\lambda,w,k}$ is tied to $q_{\lambda}$ by the following relationship (Lemma 5, ):

In other words, the inclusion probability of the corresponding projection DPP is related to the inverse of the Christoffel function as defined in (183). Figure 7 illustrates the evaluations of the inclusion probability of the projection DPP in the case of RKHS defined by the Gaussian kernel along with the Gaussian measure in the real line. Recall that in this case the eigenfunctions are given by