Monte Carlo integration of non-differentiable functions
on $[0, 1]^{l}$, $l = 1, ..., d$, using a single determinantal point
pattern defined on $[0, 1]^{d}$

Jean-François Coeurjolly, Adrien Mazoyer, Pierre-Olivier Amblard

► To cite this version:

Jean-François Coeurjolly, Adrien Mazoyer, Pierre-Olivier Amblard. Monte Carlo integration of non-differentiable functions on $[0, 1]^{l}$, $l = 1, ..., d$, using a single determinantal point pattern defined on $[0, 1]^{d}$. Electronic Journal of Statistics, 2021, 15 (2), pp.6228-6280. 10.1214/21-EJS1929. hal-02508401

HAL Id: hal-02508401

https://hal.science/hal-02508401

Submitted on 14 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Monte Carlo integration of non-differentiable
functions on $[0,1]^{\ell}$, $\ell = 1, \dots, d$, using a single
determinantal point pattern defined on $[0,1]^d$

Jean-François Coeurjolly, Adrien Mazoyer and Pierre-Olivier Amblard

February 28, 2020

Abstract

This paper concerns the use of a particular class of determinantal point processes (DPP), a class of repulsive spatial point processes, for Monte Carlo integration. Let $d \ge 1$, $I \subseteq \bar{d} = \{1, \dots, d\}$ with $\ell = |I|$. Using a single set of $N$ quadrature points $\{u_1, \dots, u_N\}$ defined, once for all, in dimension $d$ from the realization of the DPP model, we investigate “minimal” assumptions on the integrand in order to obtain unbiased Monte Carlo estimates of $\mu(f_I) = \int_{[0,1]^\ell} f_I(u)du$ for any known $\ell$-dimensional integrable function on $[0, 1]^\ell$. In particular, we show that the resulting estimator has variance with order $N^{-1-(2s)^2/d}$ when the integrand belongs to some Sobolev space with regularity $s > 0$. When $s > 1/2$ (which includes a large class of non-differentiable functions), the variance is asymptotically explicit and the estimator is shown to satisfy a Central Limit Theorem.

Introduction

In the context of computer experiments (see for example [24, Chapter 5]), complex
phenomena are simulated using a mathematical model to replace the process which
generates the data. Usually, the model depends on a large number of parameters
(inputs). An objective of the experiments is to quantify the influence of the vari-
ability of the inputs on the variable of interest. An experiment consists in running
simulations, where each simulation represents a possible combination of the inputs.
It is impossible in practice to consider all possible configurations, the number of
simulations being limited. Therefore, the design of experiments, i.e. the choice of
combinations of inputs, is of great importance. Under a lack of information on how
inputs are linked to outputs, a strategy is to spread chosen inputs to cover as much
as possible all the input space. This technique is called *space-filling design*. It can
be summarized by generating $N$ points in a given space which regularly cover this
space. Latin hypercubes [16, 20], low discrepancy sequences (see e.g. [10, 26]) are
standard methods to generate designs. The goal of computer experiments is not
only to examine the influence of all the inputs on an output of interest, but also the
influence of a subset of these inputs, or also the influence of particular combination
of subsets of these inputs. Since computer experiments may be very expensive in
terms of computation load and/or storage capacity, the regularity of the coverage
of the designs should be conserved when the initial configuration is projected onto
lower dimensional spaces. This would allow to use the initial configuration to study
the influence of subsets of inputs for example with the same efficiency.

We investigate in this paper a similar problem in the context of Monte Carlo
integration. Let $d \ge 1$, $I \subseteq \bar{d} = \{1, \dots, d\}$ with $\iota = |I|$, we aim to estimate,
under “minimal” assumptions, the integral $\mu(f_I) = \int_{[0,1]^\iota} f_I(u)du$ for any known
$\iota$-dimensional integrable function on $[0, 1]^\iota$ using a single set of $N$ quadrature points
$\{u_1, \dots, u_N\}$ defined once for all in dimension $d$. We insist on the fact that the same
set of nodes is used to estimate different integrals. Hence, this set does not exploit
the form of $f_I$, the locations of its possible singularities, etc.

Monte Carlo integration has a long history and it is not the aim of this paper
to make a detailed bibliography. We refer the interested reader by an extensive
treatment and bibliography to the electronic book [21]. Let us however cite a few
methods keeping in mind what we mean by “minimal” assumptions in the situation
$\iota = d$. Ordinary Monte Carlo methods and importance sampling methods (see e.g.
[23]) consist in using i.i.d. nodes $\{u_1, \dots, u_N\}$ with a so-called proposal density (the
uniform density in the usual situation). Under some $L^2([0, 1]^d)$ type assumption on
the integrand, the resulting estimator denoted by $\hat{\mu}_N(f_{\bar{d}})$ has a variance proportional
to $N^{-1}$ and satisfies a Central Limit Theorem. When $d$ is large, Monte Carlo Markov
Chains (MCMC) methods, where the set of quadrature points is the realization
of a particular Markov chain, are usually preferred. When $f \in L^1([0, 1]^d)$, the
variance of $\hat{\mu}_N(f_{\bar{d}})$ is still of order $N^{-1}$ and satisfies a CLT (see e.g. [8]). To
improve the rate of convergence, the price to pay is to require some regularity
assumptions on $f_{\bar{d}}$. Many methods exist in the literature: grid-based stratified
methods [9], possibly combined with antithetic sampling (see [21, Chapter 10]),
Quasi Monte Carlo and randomized versions, scrambled nets [6, 15, 22, 23], etc.
For example, a version of scrambled nets with antithetic sampling can lead to an
estimator with variance $O(N^{-3-2/d}\log(N)^{d-1})$ if, to simplify Owen's assumption
[19], $f_{\bar{d}}$ is $d$ times continuously differentiable on $[0, 1]^d$. Additional assumptions on
the scrambled net are required to obtain a CLT. Grid-based stratified methods which
are maybe the first simple alternative to ordinary Monte Carlo methods require that
$f_{\bar{d}}$ is continuously differentiable on $[0, 1]^d$ and yield an estimator satisfying a CLT
with variance asymptotically proportional to $N^{-1-2/d}$. Let us also mention that
[18] showed that a version of scrambled net has a variance $o(N^{-1})$ under the sole
assumption that $f_{\bar{d}} \in L^2([0, 1]^d)$ but the rate is not explicit until strong regularity assumptions are made on $f_{\bar{d}}$.

In a recent work, another alternative by defining the nodes as the realization of a repulsive point pattern, and in particular a determinantal point pattern, has been proposed in [3]. The class of determinantal point processes (DPPs for short) has received a growing attention in the last decades (see e.g. [27, 25, 12, 14, 6]), thanks to its very appealing properties in particular in terms of tractability and exact simulation. [3] have defined an Orthogonal Polynomial Ensemble, which is a particular inhomogeneous DPP whose kernel is defined through orthonormal polynomials. Under the assumption that $f_{\bar{d}}$ is continuously differentiable and compactly supported in $B' \subset [0, 1]^d$ the authors obtained an estimator with variance equivalent to an explicit constant times $N^{-1-1/d}$.

In this paper, we investigate a different DPP model. To be more explicit, we consider the most natural kernel, called the Dirichlet kernel in this paper, which is based on the Fourier decomposition of a rectangular subset of $N$ indices of $\mathbb{Z}^d$. It is a projection DPP which produces almost surely $N$ points. It has the advantage to lead to a homogeneous DPP pattern, an interesting characteristic as we want the pattern to be used to estimate any integral. A second advantage is that the marginals are fully characterized and explicit which means that marginals can efficiently be used to estimate $\hat{\mu}_N(f_I)$ for any $I \subset \bar{d}$. Last but not least, our main result Theorem 3.1, shows that the resulting estimator $\hat{\mu}_N(f_I)$ has asymptotic variance proportional to $N^{-1-(2s\wedge 1)/d}$ for any $f_I \in \mathcal{H}^s([0, 1]^l)$ where $\mathcal{H}^s([0, 1]^l)$ is some Sobolev space with regularity $s \ge 0$, see (3.4) for more details. We remind that for periodic functions $L^2([0, 1]^l) = \mathcal{H}^0([0, 1]^l)$ and if $f_I$ is periodic and continuously differentiable then $f_I \in \mathcal{H}^l([0, 1]^l)$. In particular, our result states that when $s > 1/2$ (thus potentially for non-differentiable functions), the variance is asymptotically equivalent to an explicit constant times $N^{-1-1/d}$. In this case, we also obtain a central limit theorem for the estimator (assuming in addition that the integrand is bounded). As a summary, the estimator proposed has the characteristic to exhibit a variance that decreases faster than the ordinary Monte Carlo as soon as $s > 0$. The decay is slower than methods such as grid-based methods, scrambled nets, etc, but require much less regularity assumptions and can be applied to any $l$-dimensional function, $l = 1, \dots, d$.

The paper is organized as follows. Section 1 contains a background on spatial point processes, generalities on the projection of spatial point processes and DPPs. We also outline the interest of repulsive point processes and in particular DPPs for Monte Carlo integration. Section 2 introduces the Dirichlet DPP and exposes some of its properties. Our main result, Theorem 3.1, is presented in Section 3. It details convergence results for Monte Carlo integration based on the realization of a Dirichlet DPP. A multivariate version of this CLT is also proposed. Section 4 discusses computational aspects for the simulation of Dirichlet DPPs and contains
a simulation study which illustrates our results. Finally, proofs of the results are postponed to Appendix B.

# 1 Background and notation

## 1.1 Spatial point processes

A spatial point process $\mathbf{X}$ defined on a Borel set $B \subseteq \mathbb{R}^d$ is a locally finite measure on $B$, see for example [17] and references therein for measure theoretical details, whose realization is of the form $\{x^{(1)}, \dots, x^{(k)}\} \in B^k$ where $k$ is the realization of a random variable and the $x^{(i)}$'s represent the events. We assume that $\mathbf{X}$ is simple meaning that two events cannot occur at the same location. Thus, $\mathbf{X}$ is viewed as a locally finite random set.

In most cases, the distribution of a point process $\mathbf{X}$ can be described by its intensity functions $\rho_{\mathbf{X}}^{(k)}: B^k \to \mathbb{R}^+$, $k \in \mathbb{N} \setminus \{0\}$. By Campbell theorem, see e.g. [17], $\rho_{\mathbf{X}}^{(k)}$ is characterized by the following integral representation: for any non-negative measurable function $h: B^k \to \mathbb{R}^+$

$$ \mathbb{E} \left[ \sum_{x^{(1)}, \dots, x^{(k)} \in \mathbf{X}} h(x^{(1)}, \dots, x^{(k)}) \right] \tag{1.1} \\ = \int_{B^k} \rho_{\mathbf{X}}^{(k)}(x^{(1)}, \dots, x^{(k)}) h(x^{(1)}, \dots, x^{(k)}) dx^{(1)} \dots dx^{(k)} $$

where $\neq$ over the summation means that $x^{(1)}, \dots, x^{(k)}$ are pairwise distinct points. Intuitively, for any pairwise distinct points $x^{(1)}, \dots, x^{(k)} \in B$, $\rho_{\mathbf{X}}^{(k)}(x^{(1)}, \dots, x^{(k)}) dx^{(1)} \dots dx^{(k)}$ is the probability that $\mathbf{X}$ has a point in each of the $k$ infinitesimally small sets around $x^{(1)}, \dots, x^{(k)}$ with volumes $dx^{(1)}, \dots, dx^{(k)}$, respectively. When $k=1$, this yields the intensity function simply denoted by $\rho_{\mathbf{X}} = \rho_{\mathbf{X}}^{(1)}$. The second order intensity $\rho_{\mathbf{X}}^{(2)}$ is used to define the pair correlation function

$$ g_{\mathbf{X}}(x^{(1)}, x^{(2)}) = \frac{\rho_{\mathbf{X}}^{(2)}(x^{(1)}, x^{(2)})}{\rho_{\mathbf{X}}(x^{(1)}) \rho_{\mathbf{X}}(x^{(2)})} \tag{1.2} $$

for pairwise distinct $x^{(1)}, x^{(2)} \in B$ and where $g_{\mathbf{X}}(x^{(1)}, x^{(2)})$ is set to 0 if $\rho_{\mathbf{X}}(x^{(1)})$ or $\rho_{\mathbf{X}}(x^{(2)})$ is zero. By convention, $\rho_{\mathbf{X}}^{(k)}(x^{(1)}, \dots, x^{(k)})$ is set to 0 if $x^{(i)} = x^{(j)}$ for some $i \neq j$. Therefore $g_{\mathbf{X}}(x, x)$ is also set to 0 for all $x \in B$ by convention. The pair correlation function (pcf for short) can be used to determine the local interaction between points of $\mathbf{X}$ located at $x$ and $y$: $g_{\mathbf{X}}(x, y) > 1$ characterizes positive correlation between the points; $g_{\mathbf{X}}(x, y) = 1$ means there is no interaction (typically a Poisson point process); $g_{\mathbf{X}}(x, y) < 1$ characterizes negative correlations. A point pattern is often referred to as a repulsive point process, if $g_{\mathbf{X}}(x, y) < 1$ for
any $x, y \in B$ (see e.g. [13, Section 6.5]). Finally, a point process $\mathbf{X}$ with constant intensity function on $B$ is said to be homogeneous.

## 1.2 Projection of a spatial point process

In this work, we sometimes consider projection of spatial point processes. By projection, we mean that we keep a given number of coordinates from the original spatial point process. Such a framework requires that the original point process $\mathbf{X}$ must be defined on a compact set $B \subset \mathbb{R}^d$: otherwise, the configuration of points of the projected point processes may not form a locally finite configuration, as also noticed in the two-dimensional case in [1, p. 17].

This section presents a few notation in this context. Let $I \subseteq \bar{d} := \{1, \dots, d\}$ with cardinality $|I| = \iota$. Let $B_1, \dots, B_d$ compact sets of $\mathbb{R}$ and $B = B_1 \times \dots \times B_d$, denote by $B_I$ its orthogonal projection onto $\mathbb{R}^\iota$. In particular $B_I = \prod_{i \in I} B_i$ with $B = B_{\bar{d}}$. We denote by $P_I$ the orthogonal projection of $\mathbb{R}^d$ onto $\mathbb{R}^\iota$. To ease the reading, we let $B_\ell = B_0$ for $\ell = 1, \dots, d$ and even to fix ideas let $B_0 = [0,1]$. Thus, $B_I = B_0^\iota = [0,1]^\iota$. For any $x \in B$, we let $x_I = P_I x$ and for a point process $\mathbf{X}$ defined on $B$, the projected point process $\mathbf{X}_I = P_I \mathbf{X}$ is then defined on $B_I$. Intensity functions and Laplace functionals for $P_I \mathbf{X}$ can be derived from the corresponding functions and functionals from $\mathbf{X}$, see [15] for more details and Section 2.2 for the particular model considered in this paper.

## 1.3 Determinantal point processes

In this section, the class of continuous DPPs is introduced. Again, we restrict our attention to DPPs defined on a compact set $B \subset \mathbb{R}^d$. A point process $\mathbf{X}$ on $B$ is said to be a DPP on $B$ with kernel $K: B \times B \to \mathbb{C}$ if for any $k \ge 1$ its $k$th order intensity function is given by

$$ \rho_{\mathbf{X}}^{(k)}(x^{(1)}, \dots, x^{(k)}) = \det \left[ K(x^{(i)}, x^{(j)}) \right]_{i,j=1}^{k} \quad (1.3) $$

and we simply denote by $\mathbf{X} \sim \text{DPP}_B(K)$. Note that $K$ needs to be non-negative definite to ensure $\rho_{\mathbf{X}}^{(k)} \ge 0$. Our results rely on the spectral decomposition of $K$, see (1.5). Therefore, we assume that $K$ is a continuous covariance function. The intensity of $\mathbf{X}$ is given by $\rho_{\mathbf{X}}(x) = K(x, x)$ and its pcf by

$$ g_{\mathbf{X}}(x, y) = 1 - \frac{|K(x, y)|^2}{K(x, x)K(y, y)}. \quad (1.4) $$

The popularity of DPPs relies mainly upon (1.3)-(1.4): all moments of $\mathbf{X}$ are explicit and since $K$ is Hermitian, $g_{\mathbf{X}}(x, y) < 1$ for any $x, y \in B$. The kernel $K$ defines an
integral operator $\mathcal{K}$ (see e.g. [5]) defined for any $f \in L^2(B)$ by

$$ \mathcal{K}(f)(x) = \int_B K(x, y)f(y)dy, \quad x \in B. $$

From Mercer's Theorem [22, Sec. 98], $\mathcal{K}$ admits the following spectral decomposition for any $x, y \in B$

$$ K(x, y) = \sum_{j \in \mathbb{N}} \lambda_j \phi_j(x) \overline{\phi_j(y)} \quad (1.5) $$

where $\{\phi_j\}_{j \in \mathbb{N}}$ are eigenfunctions associated to $\mathcal{K}$ and form an orthonormal basis of $L^2(B)$, and where $\{\lambda_j\}_{j \in \mathbb{N}}$ are the eigenvalues of $\mathcal{K}$ satisfying $\lambda_j \ge 0$ for any $j \in \mathbb{N}$. We abuse notation in the sequel and refer $\lambda_j$'s to as the eigenvalues of $\mathcal{K}$.

The existence of a DPP on $B$ with kernel $\mathcal{K}$ is ensured if its eigenvalues satisfy $\lambda_j \le 1$ for any $j \in \mathbb{N}$, see e.g. [12, Theorem 4.5.5.]. Eigenvalues and eigenfunctions are indexed here by $\mathbb{N}$ in (1.5), but other countable sets could be considered. In particular, the $d$-dimensional Fourier basis is indexed by $\mathbb{Z}^d$. A DPP $\mathbf{X} \sim \text{DPP}_B(\mathcal{K})$ is said to be homogeneous if $\mathcal{K}$ is the restriction on $B \times B$ of a kernel $C$ defined on $\mathbb{R}^d \times \mathbb{R}^d$ which satisfies $C(x, y) = C(o, x - y)$ for any $x, y \in \mathbb{R}^d$ where $o$ is the origin in $\mathbb{R}^d$. In that case, we will refer to $\mathcal{K}$ as a stationary kernel and will use the abusive notation $K(x, y) \equiv K(x - y)$.

A kernel $\mathcal{K}$ such that $\lambda_j \in \{0, 1\}$ for $j \in \mathbb{N}$ is called a “projection kernel” and the corresponding DPP a “projection DPP”. The number of points in $B$ of such a model is almost surely constant and equal to the number of non-zero eigenvalues of $\mathcal{K}$ (see e.g. [14]).

## 1.4 Why are DPPs interesting for Monte Carlo integration?

The repulsive nature of DPPs can be exploited to generate quadrature points that explore nicely the input space. To see this, let $B \subset \mathbb{R}^d$ be a bounded set, $f \in L^2(B)$ and $\mathbf{Y}$ a homogeneous point process on $B$ with intensity parameter $\rho_{\mathbf{Y}}$ and pair correlation function $g_{\mathbf{Y}}$ (a similar result would hold in the inhomogeneous case). Campbell's Theorem (1.1) ensures that the estimator

$$ \hat{\mu}(f) = \frac{1}{\rho_{\mathbf{Y}}} \sum_{u \in \mathbf{Y}} f(u) \quad (1.6) $$

is an unbiased estimator of $\mu(f) = \int_B f(u)du$ with variance

$$ \mathrm{Var}[\hat{\mu}(f)] = \frac{1}{\rho_{\mathbf{Y}}} \int_B f(u)^2 du + \int_{B^2} (g_{\mathbf{Y}}(u, v) - 1)f(u)f(v)du dv. \quad (1.7) $$

If $f$ is non-negative (or non-positive), Equation (1.7) suggests that using a point processes satisfying $g_{\mathbf{Y}} < 1$ makes the variance smaller than the first term which
turns out to be the variance under the Poisson case. It is worth noting that the use of a DPP for this task does not require any sign assumption for $f$. Indeed, given the fact that $1 - g_Y(u, v) = |K(u, v)|^2/\rho_Y^2$ and from Mercer's decomposition (1.5), we obtain

$$ \mathrm{Var}[\hat{\mu}(f)] = \frac{1}{\rho_Y} \int_B f(u)^2 du - \frac{1}{\rho_Y^2} \sum_{j,k \in \mathbb{N}} \lambda_j \lambda_k \left| \int_B f(u) \phi_j(u) \bar{\phi}_k(u) du \right|^2 \quad (1.8) $$

and so the second term is always negative.

The use of a general DPP does not seem to be of great interest. First, the number of points is random and thinking $\rho_Y$ as a number of points, we claim that the rate of convergence remains the same as in the independent case.

Therefore it seems natural to focus on the subclass of projections DPPs, i.e. a class for which the number of points is almost surely constant. An approach using an ad-hoc Orthogonal Polynomial Ensemble with $N$ points has been proposed in [3]. This class is a particular inhomogeneous projection DPP on $B$. As already outlined in the introduction, under the assumption that $f \in C^1(B)$ and that $f$ is compactly supported on some bounded $B' \subset B$, it is proved that a central limit theorem holds for the integral estimator with variance decreasing as $N^{-1-1/d}$. We propose a similar approach, based on a realization of an $(N, d)$-Dirichlet DPP $X$. The $(N, d)$-Dirichlet DPP, detailed in the next section, is a projection DPP based on the Fourier basis. Unlike, the model proposed by [3], this DPP has the advantage to be homogeneous and its projections $X_I$ for $I \subseteq \{1, \dots, d\}$ are fully characterized (see Section 2.2). Finally, an advantage of our approach is that we do not require that $f$ is continuously differentiable. We only assume that $f$ belongs to some Sobolev space with low regularity parameter, see Section 3 and in particular Theorem 3.1 for more details.

# 2 The $(N, d)$-Dirichlet DPP and its projections

## 2.1 The $(N, d)$-Dirichlet DPP

Let us consider the Fourier basis in $B = [0,1]^d$ defined for any $j \in \mathbb{Z}^d$, $x \in B$ by $\phi_j^{(d)}(x) = e^{2i\pi j^\top x}$. Given a vector of $d$ positive integers $n = (n_i)_{i=1...d}$, we construct the following kernel

$$ K(x, y) = \sum_{j \in E_N} \phi_j^{(d)}(x) \overline{\phi_j^{(d)}(y)} = \sum_{j \in E_N} e^{2i\pi j^\top (x-y)} \quad (2.1) $$

where

$$ E_N = E_1 \times \dots \times E_d \quad \text{and for } i = 1 \dots d, \quad E_i = \{0, 1, \dots, n_i - 1\}. \quad (2.2) $$
Thus, $E_N$ is the rectangular subset of $\mathbb{Z}^d$ with cardinality $N$ which identifies eigenvalues that are all equal to 1. Due to the invariance by translation of the Fourier basis, the kernel (2.1) is a homogeneous kernel. This construction implies that for any $x, y \in [0,1]^d$, $K(x - y) = \prod_{i=1}^d K_i(x_i - y_i)$ where the $K_i$'s are one-dimensional stationary kernels defined for any $x_i, y_i \in [0,1]$ by $K_i(x_i - y_i) = \sum_{j=0}^{n_i-1} e^{2i\pi j(x_i - y_i)}$. We point out that defining $E_i$ as a block of successive $n_i$ frequencies (e.g. frequencies centered around 0), would lead to the same DPP [12, Remark 4, p.48]. In particular, if $n_i$ is odd, we could consider $E_i = \{-\lfloor n_i/2 \rfloor, \dots, \lfloor n_i/2 \rfloor\}$, which leads to the standard Dirichlet kernel, see e.g. [30]. This justifies the name Dirichlet DPP for this model.

Such a DPP, which produces almost surely $N = \prod_{i=1}^d n_i$ points in $B$, will be referred to as an $(N, d)$-Dirichlet DPP. We could wonder why we impose $E_N$ to be a rectangular subset of $\mathbb{Z}^d$ instead of, for instance, the graded lexicographic order used in [3]. As seen in the next section, the rectangular nature of $E_N$ allows us to characterize the distribution of $X_I$ for any $I \subseteq \bar{d}$.

Let us add that, when $d = 1$, the kernel $K$ corresponds to the Fourier approximation [14] of the one-dimensional sine-kernel $\sin(\pi n_i(x_i - y_i))/(\pi(x_i - y_i))$, which takes its origins in the joint distribution of the eigenvalues (called the Weyl measure) of a unitary matrix. Asymptotic results involving one-dimensional linear functionals from the sine-kernel appear in several papers, see e.g. [28]. The present paper provides therefore an extension to the $d$-dimensional case and a more thorough treatment of the statistical application to Monte Carlo integration.

## 2.2 Projections of an $(N, d)$-Dirichlet kernel

An $(N, d)$-Dirichlet kernel (2.1) can be written as the product of $d$ one-dimensional kernels. More precisely, for any $I \subset \bar{d}$, by denoting $N_I = \prod_{i \in I} n_i$ and $N_{I^c} = N/N_I$, the $(N, d)$-Dirichlet kernel can always be written as

$$K(x - y) = K_I(x_I - y_I)K_{I^c}(x_{I^c} - y_{I^c}) \quad (2.3)$$

where $K_I$ (resp. $K_{I^c}$) is the $(N_I, \iota)$-Dirichlet kernel (resp. $(N_{I^c}, d - \iota)$-Dirichlet kernel). Projected point processes $X_I$ from models with kernels satisfying (2.3) have been studied and characterized in [15]. In particular, for an $(N, d)$-Dirichlet DPP we have the following result.

**Proposition 2.1.** Let $X$ be an $(N, d)$-Dirichlet DPP on $B$, let $I \subseteq \{1, \dots, d\}$, then $X_I$ is an $(-1/N_{I^c})$-DPP on $B_I$ with kernel $N_{I^c}K_I$, i.e. $X_I \sim (-1/N_{I^c})$-DPP$_{B_I}(N_{I^c}K_I)$. In particular, $X_I$ has (obviously) $N$ points and $k$-th order intensity

$$\rho_{X_I}(x^{(1)}, \dots, x^{(k)}) = \det_{-1/N_{I^c}} [N_{I^c} K_I(x^{(i)}, x^{(j)})]_{i,j=1}^k$$
for any pairwise distinct $x^{(1)}, \dots, x^{(k)} \in B_I$. Its pcf is therefore given, for any pairwise distinct $x, y \in B_I$, by

$$g_{\mathbf{X}_I}(x, y) = 1 - \frac{|K_I(x, y)|^2}{NN_I}. \qquad (2.4)$$

In the above result, the notation $\det_\alpha$ stands for an $\alpha$ determinant, see e.g. [25] for details on such quantities and for general properties of $(\alpha)$-DPPs. Proposition 2.1 therefore proposes a full characterization of the distribution of $\mathbf{X}_I$. Its main consequence (and interest for the present paper), as seen from (2.4), is that the pcf of $\mathbf{X}_I$ is bounded by 1 and therefore, for any $I$, $\mathbf{X}_I$ remains in the class of repulsive point patterns.

# 3 Numerical integration with Dirichlet kernel

## 3.1 Objective

In this section, we study the use of specific DPPs for Monte Carlo integration. To this end, we use notation introduced in Section 1.2. Our objective is to estimate any $\iota$-dimensional integral, for $1 \le \iota \le d$, using a Monte Carlo approach and using the same quadrature points. More precisely, let $d \ge 1$, $I \subset \bar{d} = \{1, \dots, d\}$ with cardinality $\iota = |I|$ and let $f_I: B_I \to \mathbb{R}$ be a measurable function on $B_I = [0,1]^\iota$, such that $f_I \in L^2(B_I)$. More assumptions on $f_I$ will be given later. We intend to estimate

$$\mu(f_I) = \int_{B_I} f_I(u) du$$

using the projection onto $B_I$ of $\mathbf{X}$ an $(N, d)$-Dirichlet kernel on $B$. In particular, we estimate $\mu(f_I)$ by

$$\hat{\mu}_N(f_I) = \frac{1}{N} \sum_{u \in \mathbf{X}_I} f_I(u) = \frac{1}{N} \sum_{j=1}^{N} f_I((u_j)_I) \qquad (3.1)$$

where $\mathbf{X} = \{u_1, \dots, u_N\}$ is an $(N, d)$-Dirichlet DPP on $B$ and where we remind the notation $(u_j)_I = P_I u_j$ for any $u_j \in B$. In the following, we study asymptotic properties for $\hat{\mu}_N(f_I)$.

In this paper, we have chosen to focus on integrals on $[0,1]^\iota$ for simplicity. It can straightforwardly be extended to rectangles. Indeed, let $a = (a_i)_{i=1,\dots,d}$, $b = (b_i)_{i=1,\dots,d} \in \mathbb{R}^d$ such that $a_i < b_i$, $i = 1, \dots, d$, and $R = [a_1, b_1] \times \dots \times [a_d, b_d]$. An estimate of $\int_{R_I} f_I(u) du$ where now $f_I \in L^1(R_I)$ is simply given by

$$\hat{\mu}_N(f_I) = \prod_{i \in I} (b_i - a_i) \left\{ \frac{1}{N} \sum_{i=1}^{N} f_I(a_I + (b_I - a_I)(u_j)_I) \right\}. \qquad (3.2)$$
We introduce two additional fundamental pieces of notation induced by the choice of the Fourier basis and used in our results. Let $f_I \in L^2(B_I)$. The notation $\hat{f}_I(j)$ for $j \in \mathbb{Z}^l$ stands for the $j$th Fourier coefficient, i.e.

$$ \hat{f}_I(j) = \int_{[0,1]^l} f_I(u)e^{-2i\pi j^\top u} du. \qquad (3.3) $$

Finally, we define the space $\mathcal{H}^s(B_I)$ as the isotropic (with respect to the sup norm) Sobolev space with index $s \ge 0$ of squared integrable periodic functions by

$$ \mathcal{H}^s(B_I) = \left\{ f_I \in L_{\text{per}}^2(B_I) : \sum_{j \in \mathbb{Z}^l} (1 + \|j\|_{\infty})^{2s} |\hat{f}(j)|^2 < \infty \right\} \quad (3.4) $$

where $L_{\text{per}}^2(B_I)$ is the set of squared integrable periodic functions. By $f_I$ periodic, we mean that all one-dimensional marginals of $f_I$ say $\hat{f}_{I,k}$, $k = 1, \dots, l$ are such that $\tilde{f}_{I,k}(0) = \tilde{f}_{I,k}(1)$. To remind some known facts about $\mathcal{H}^s(B_I)$, let us define the following spaces

$$ \begin{aligned} C_{\text{per}}^{1}(B_{I}) &= \{f_{I} \text{ periodic}, f_{I} \in C^{1}(B_{I})\} \\ \tilde{\mathcal{H}}^{3/2+\varepsilon}(B_{I}) &= \{f_{I} \in L_{\text{per}}^{2}(B_{I}) : \hat{f}_{I}(j) = \mathcal{O}(\|j\|_{\infty}^{-2-\varepsilon})\}, \quad \varepsilon > 0. \end{aligned} $$

Of course, $\mathcal{H}^{3/2+\varepsilon}(B_I)$ and $\tilde{\mathcal{H}}^{3/2+\varepsilon}(B_I)$ are very similar. Now, we remind that for any $\varepsilon > 0$

$$ \tilde{\mathcal{H}}^{3/2+\varepsilon}(B_I) \subset C_{\text{per}}^{1}(B_I) \subset \mathcal{H}^{1}(B_I). \qquad (3.5) $$

Our main results expressed by Theorem 3.1 and Corollary 3.1, have the interest to focus on functions $f_I \in \mathcal{H}^s(B_I)$ with small regularity parameter $s > 0$, thus, according to (3.5), potentially to non-differentiable functions. In particular, we show that a CLT can be expected as soon as $s > 1/2$.

Finally, let us justify why we focus on periodic functions on $B_I$. For a non periodic function which is at least twice continuously differentiable it is a known fact that $|\hat{f}_I(j)| = \mathcal{O}(\|j\|_{\infty}^{-1})$. So, the summability condition in (3.4) for such a smooth function would be fulfilled only for $s < 1/2$, which constitutes a situation of no interest.

## 3.2 Case $l = d$

We first consider the $d$-dimensional case. Let $f \in L_{\text{per}}^2(B)$ with $B = [0,1]^d$. The first result shows how the variance given in (1.8) relates to the Fourier coefficients of $f$ in the case of an $(N, d)$-Dirichlet DPP.
**Proposition 3.1.** Let $\mathbf{X}$ be an $(N,d)$-Dirichlet DPP and $f \in L_{\text{per}}^2(B)$. For any $j \in \mathbb{Z}^d$, we let $c_j(f)$ denote the jth Fourier coefficient of $f$. Then, $\hat{\mu}_N(f)$ given by (3.1) is an unbiased estimator of $\mu(f) = \int_B f(u)du$ with variance given by

$$ \mathrm{Var}[\hat{\mu}_N(f)] = \frac{1}{N} \sum_{j \in \mathbb{Z}^d} |\hat{f}(j)|^2 - \frac{1}{N^2} \sum_{j,k \in E_N} |\hat{f}(j-k)|^2 \quad (3.6) $$

$$ = \frac{1}{N} \sum_{j \in \mathbb{Z}^d} |\hat{f}(j)|^2 - \frac{1}{N^2} \sum_{j \in F_N} \left[ \prod_{i=1}^{d} (n_i - |j_i|) \right] |\hat{f}(j)|^2 \quad (3.7) $$

where $F_N = \{j \in \mathbb{Z}^d : |j_i| \le n_i - 1, i = 1, \dots, d\}$.

This simple form of the variance of $\hat{\mu}_N(f)$ invites us to study its asymptotic behavior as $N \to \infty$. Given the fact that $N = \prod_i n_i$, we require a specific asymptotic. For $i = 1, \dots, d$, we assume that $(N_\nu)_{\nu \ge 1}$ and $(n_{i,\nu})_{\nu \ge 1}$ are integer sequences indexed by some $\nu \ge 1$. We assume that each sequence $n_{i,\nu}$ tends to $\infty$ as $\nu \to \infty$ and that there exist $\kappa_i > 0$ for $i = 1, \dots, d$ such that

$$ \lim_{\nu \to \infty} n_{i,\nu} N_{\nu}^{-1/d} = \kappa_i. \quad (3.8) $$

For the sake of conciseness, we skip the dependence in $\nu$. Similarly, when we write $N \to \infty$, we implicitly assume (3.8).

We can now obtain the following asymptotic behavior of the variance of $\hat{\mu}_N(f)$ and for some values of $s$ a central limit theorem. Regarding this asymptotic normality, when $d = 1$, as mentioned in Section 2, $\mathbf{X}$ corresponds, up to a normalization, to the joint distribution of the eigenvalues of a unitary matrix. Linear statistics for such a DPP have been deeply studied in the literature, see e.g. [28] and the references therein. When $d > 1$, Theorem 3.1 (iii) is therefore original, and relies upon [29, Theorem 1].

**Theorem 3.1.** Consider the asymptotic framework (3.8) and assume that $f \in \mathcal{H}^s(B)$. Then, we have the following statements.

(i) If $s \in (0,1/2)$, then as $N \to \infty$

$$ \mathrm{Var}[\hat{\mu}_N(f)] = \mathcal{O}\left(N^{-1-\frac{2s}{d}}\right). $$

(ii) If $s \ge 1/2$ for $d = 1$ or $s > 1/2$ for $d > 1$, then

$$ \lim_{N \to \infty} N^{1+1/d} \mathrm{Var}[\hat{\mu}_N(f)] = \sigma^2(f) = \sum_{j \in \mathbb{Z}^d} \left( \sum_{i=1}^{d} \frac{|j_i|}{\kappa_i} \right) |\hat{f}(j)|^2 \quad (3.9) $$

where $\hat{f}(j)$ is given by (3.3).
$$ (iii) \text{ If in addition, } \|f\|_{\infty} < \infty, \text{ then as } N \to \infty, \\ \sqrt{N^{1+1/d}} (\hat{\mu}_N(f) - \mu(f)) \to N(0, \sigma^2(f)) \quad (3.10) $$

in distribution.

Let us rephrase Theorem 3.1 (i): if $f \in \mathcal{H}^s(B)$ for some $s > 0$, then necessarily $\mathrm{Var}[\hat{\mu}_N(f)] = o(N^{-1})$. That is, the variance of the estimator proposed decreases to 0 faster than the standard Monte Carlo estimator, as soon as $|\hat{f}(j)| = \mathcal{O}(\|j\|_{\infty}^{-\epsilon})$ for some $\epsilon > 0$, which is a very weak assumption.

Theorem 3.1 (ii) shows that the variance $\sigma^2(f)$ can be consistently estimated by $\hat{\sigma}_N^2(f)$ given by

$$ \hat{\sigma}_N^2(f) = N^{1/d} \sum_{j \in F_N} \left( \sum_{i=1}^d \frac{|j_i|}{n_i} \right) |\hat{f}(j)|^2 \quad (3.11) $$

which has the interest to avoid the constants $\kappa_i$. Hence, using Slutsky's lemma, we also have the more practical central limit theorem

$$ \sqrt{N^{1+1/d}} \frac{\hat{\mu}_N(f) - \mu(f)}{\hat{\sigma}_N(f)} \to N(0, 1) \quad (3.12) $$

in distribution, as $N \to \infty$.

## 3.3 Case $I \subset \bar{d}$

We now consider the situation where we estimate $\mu(f_I)$ (on $B_I$) based on $\{u_1, \dots, u_N\}$ which is an $(N, d)$-Dirichlet DPP on $B$. In this section, we naturally assume that $d > 1$. The interest of the construction of our model is revealed by Corollary 3.1 which, briefly, states that Theorem 3.1 can be applied to functions of the form $f_I^\uparrow(x) = f_I(x_I)1[x_{I^c} \in B_{I^c}], x \in B$.

**Corollary 3.1.** Let $d > 1$ and $I \subset \{1, \dots, d\}$ with cardinality $\iota > 0$. Consider the asymptotic framework (3.8) assume that $f_I \in \mathcal{H}^s(B_I)$, then the following statements hold.

(i) If $s \in (0,1/2)$, then as $N \to \infty$

$$ \mathrm{Var}[\hat{\mu}_N(f_I)] = \mathcal{O}\left(N^{-1-\frac{2s}{d}}\right). $$

(ii) If $s > 1/2$, then

$$ \lim_{N \to \infty} N^{1+1/d} \mathrm{Var}[\hat{\mu}_N(f_I)] = \sigma^2(f_I) \quad (3.13) $$
where

$$ \sigma^2(f_I) = \sum_{j \in \mathbb{Z}^d} \left( \sum_{i \in I} \frac{|j_i|}{\kappa_i} \right) |\hat{f}_I(j)|^2 \quad (3.14) $$

and where $\hat{f}_I(j)$ is given by (3.3).

(iii) If, in addition $f_I$ is bounded, then as $N \to \infty$

$$ \sqrt{N^{1+1/d}} (\hat{\mu}_N(f_I) - \mu(f_I)) \to N(0, \sigma^2(f_I)) \quad (3.15) $$

in distribution.

The asymptotic constant $\sigma^2(f_I)$ can still be consistently estimated by (3.11). The proof of this result is a straightforward consequence of Theorem 3.1. Another approach using the fact that $X_I$ is distributed as an $\alpha$-DPP is proposed in Appendix B.

## 3.4 Multivariate central limit theorem

We can combine Theorem 3.1 and Corollary 3.1 to obtain a multivariate version of the central limit theorem.

**Corollary 3.2.** Let $p \ge 1$. For any $\ell = 1, \dots, p$, let $I_\ell \subseteq \{1, \dots, d\}$ with $\iota_\ell = |I_\ell|$ and assume that $f_{I_\ell} \in \mathcal{H}^s(B_{I_\ell})$ with $s \ge 1/2$ if $d=1$ or $s > 1/2$ if $d > 1$. Let $\hat{\mu}_{N,p} = (\hat{\mu}_N(f_{I_1}), \dots, \hat{\mu}_N(f_{I_p}))^\top$ and $\mu_p = (\mu(f_{I_1}, \dots, \mu(f_{I_p})))^\top$. Then, under the asymptotic framework (3.8), $\hat{\mu}_{N,p}$ is an unbiased estimator of $\mu_p$ and as $N \to \infty$

$$ \sqrt{N^{1+1/d}} (\hat{\mu}_{N,p} - \mu_p) \to N(0, \Sigma_p) $$

in distribution, where $\Sigma_p$ is the $(p, p)$ Hermitian matrix with entries

$$ (\Sigma_p)_{\ell\ell'} = \sum_{j \in \mathbb{Z}^d} \left( \sum_{i=1}^d \frac{|j_i|}{\kappa_i} \right) \hat{f}_{I_\ell}^\top(j) \overline{\hat{f}_{I_{\ell'}}^\top(j)} $$

where $\hat{f}_{I_\ell}^\top(x) = f_{I_\ell}(x_I)1[x_{I_\ell}^c \in B_{I_\ell}^c], x \in B$.

## 4 Simulation study

We propose now a simulation study to illustrate the results. We first consider the setting of Theorem 3.1 (ii)-(iii) and Corollary 3.1 (ii)-(iii), i.e. on the case $s > 1/2$ (when $d > 1$) and second the situation $s < 1/2$.
## 4.1 Case $s > 1/2$, illustration of Theorem 3.1 (ii)-(iii)

We consider three different functions with different regularity properties:

* Bump function

$$f_{\text{bump}}(x) = \prod_{i=1}^{d} \frac{\varphi_{\text{bump}}(x_i)}{\int_{0}^{1} \varphi_{\text{bump}}(t)dt}, \quad \varphi_{\text{bump}}(t) = \exp \left( -\frac{0.4}{1 - (t - 1/2)^2} \right). \quad (4.1)$$

* Mixture of cosines

$$f_{\text{mixcos}}(x) = \frac{1}{d} \sum_{i=1}^{d} \frac{\varphi_{\text{mc}}(x_i)}{\int_{0}^{1} \varphi_{\text{mc}}(t)dt}, \quad \varphi_{\text{mc}}(t) = 0.1 |\cos(5\pi(t - 1/2))| + (t - 1/2)^2. \quad (4.2)$$

* Normalized $L^\gamma$-norm: let $\gamma > 0$.

$$f_{\gamma}(x) = \frac{1}{d} \sum_{i=1}^{d} \frac{\varphi_{\gamma}(x_i)}{\int_{0}^{t} \varphi_{\gamma}(t)dt}, \quad \varphi_{\gamma}(t) = |t - 1/2|^{\gamma}. \quad (4.3)$$

It is worth mentioning that $f_{\text{bump}}$ is infinitely continuously differentiable, so $f_{\text{bump}} \in \mathcal{H}^s$, for any $s > 0$. The function $f_{\text{mixcos}}$ is a non-differentiable (with $4^d$ singularity points) which satisfies $\hat{f}_{\text{bump}}(j) = \mathcal{O}(\|j\|_{\infty}^{-2})$, so $f_{\text{mixcos}} \in \mathcal{H}^s(B)$ for any $s < 3/2$. Finally, $f_{\gamma}$ is also a non-differentiable squared integrable periodic function, which satisfies $\hat{f}_{\gamma}(j) = \mathcal{O}(\|j\|_{\infty}^{-1-\gamma})$. Hence, $f_{\gamma} \in \mathcal{H}^s(B)$ for any $s < 1/2 + \gamma$. In the following, we consider the cases $\gamma = 0.25$ and $\gamma = 0.75$. These four functions are depicted for $d = 2$ in Figure 1. We perform the following experiment. For $N = 50, 100, 150, \dots$,

500, 600, \dots, 1000 (i.e. 15 values for $N$) and for $d = 1, \dots, 6$, we generate 2500 realizations of the $(N, d)$-Dirichlet DPP. The factorization $N = \prod_{i=1}^{d} n_i$ is set such that the fluctuation of the $n_i$'s is minimized. For example, for $N = 100$ and $d = 2$, we set $n = (10, 10)$ while when $d = 6$ we choose $n = (5, 5, 2, 2, 1, 1)$.

We use the simulation algorithm provided in [11]. Basically, it relies upon a Gram-Schmidt orthogonalization of $N$ vectors with dimension $N$, with a cost of order $N^3$, and a rejection sampling step. Sampling DPPs with large $N$ or when $d > 4$ is very time consuming using the standard R package `spatstat` [2]. Therefore, we have reimplemented an R package based on C++ functions (available on https://github.com/AdriMaz/rcdpp/). We perform our experiments in very reasonable computing time, with a 2,3 GHz Intel Core i5 processor and only 8 Go (2133 MHz DDR4) of RAM. For example, sampling a DPP with $N = 1000$ in dimension 6 can be performed in a few seconds.
<figure><img src="image_2.png" /><figcaption>Figure 1: Test functions considered in the simulation study, given by (4.1)-(4.3) depicted in dimension $d = 2$.</figcaption></figure>
For each $d$, each function and each replication of the point pattern, we evaluate the estimator $\hat{\mu}_N(f_I)$ given by (3.1) with $I = \bar{d}$. To visualize the rate of convergence of the variance, we perform a linear regression of the logarithm of empirical variances in terms of $\log(N)$. According to Theorem 3.1 (ii), the expected slope is $-1 - 1/d$. For each test function, $d$ and $N$, we also test the normality of estimates using the Shapiro-Wilk test. After adjusting the p-values using Holm procedure for each function and each $d$, we find that none of the p-values are smaller than 5%. Therefore the normality assumption is never rejected.

Results are illustrated by Figure 2 and are in a clear agreement with our theoretical result. The asymptotic normality is not rejected even for small sample size $N$, and the slope of the logarithm of empirical variances is very close to $-1 - 1/d$. Note that the fourth columns of the summary tables expose the usual Student-t confidence interval for the slopes. For the sake of clarity of the plots, all curves have been translated in order to compare more easily the slopes. Thus for this figure, the y-axis has no meaning.

## 4.2 Case $s > 1/2$, illustration of Corollary 3.1 (ii)-(iii)

We perform similar experiments. We set the dimension to $d = 6$. For each configuration $\{x^{(1)}, \dots, x^{(N)}\}$ of $(N, 6)$-Dirichlet DPPs, we evaluate the estimator $\hat{\mu}_N(f_I)$ given by (3.1) with $I \subset \bar{d}$ and $|I| = \iota = 1, \dots, 6$. In other words, we use 6-dimensional configurations of points to estimate integrals of $\iota$-dimensional integrands. Let us precise that a single realization $(N, 6)$-Dirichlet is used for each value of $\iota$. However, the directions kept when projecting are chosen randomly. Results are illustrated in Figure 2, as for the previous case. This time, there are several values of $N$ for which the normality assumptions has been rejected, i.e. for which adjusted $p$-value is smaller than 0.05. These points are represented by crosses instead of regular dots. We still made the arbitrary choice to keep them when performing the regression lines. The conclusions are quite similar to previous case: points remain nicely aligned along regression lines and confidence intervals keep to be in agreement with theoretical slopes, which in this situation are all equal to $-1 - 1/6 \approx -1.17$. It seems that the normality is hardly hinted when a design is projected on low dimensional space: most of the values of $N$ are rejected by the Shapiro-Wilk tests for $\iota = 1, 2$.

## 4.3 Case $s < 1/2$

We now intend to illustrate Theorem 3.1 (i) for $d = 1$ and $s < 1/2$. We consider the following one-dimensional function defined on $[0, 1]$:

$$h_{\gamma}(x) = \sum_{j \ge 1} \frac{\cos(2\pi j(x - 1/2))}{2\pi j^{\gamma}}. \qquad (4.4)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Summary of experiments in which integrals of *d*-dimensional functions are estimated using an (*N*, *d*)-Dirichlet DPP.</figcaption></figure>

<table><thead><tr><th>d</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>1</td><td>-2.00</td><td>-2.00</td><td>[-2.01; -1.98]</td></tr><tr><td>2</td><td>-1.50</td><td>-1.50</td><td>[-1.53; -1.47]</td></tr><tr><td>3</td><td>-1.33</td><td>-1.35</td><td>[-1.37; -1.32]</td></tr><tr><td>4</td><td>-1.25</td><td>-1.26</td><td>[-1.28; -1.24]</td></tr><tr><td>5</td><td>-1.20</td><td>-1.18</td><td>[-1.21; -1.14]</td></tr><tr><td>6</td><td>-1.17</td><td>-1.12</td><td>[-1.15; -1.09]</td></tr></tbody></table>

<table><thead><tr><th>d</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>1</td><td>-2.00</td><td>-1.98</td><td>[-2; -1.96]</td></tr><tr><td>2</td><td>-1.50</td><td>-1.48</td><td>[-1.51; -1.46]</td></tr><tr><td>3</td><td>-1.33</td><td>-1.33</td><td>[-1.35; -1.31]</td></tr><tr><td>4</td><td>-1.25</td><td>-1.22</td><td>[-1.25; -1.2]</td></tr><tr><td>5</td><td>-1.20</td><td>-1.18</td><td>[-1.22; -1.15]</td></tr><tr><td>6</td><td>-1.17</td><td>-1.14</td><td>[-1.17; -1.1]</td></tr></tbody></table>

<table><thead><tr><th>d</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>1</td><td>-2.00</td><td>-2.01</td><td>[-2.03; -1.98]</td></tr><tr><td>2</td><td>-1.50</td><td>-1.51</td><td>[-1.53; -1.48]</td></tr><tr><td>3</td><td>-1.33</td><td>-1.34</td><td>[-1.37; -1.32]</td></tr><tr><td>4</td><td>-1.25</td><td>-1.31</td><td>[-1.33; -1.28]</td></tr><tr><td>5</td><td>-1.20</td><td>-1.25</td><td>[-1.29; -1.2]</td></tr><tr><td>6</td><td>-1.17</td><td>-1.19</td><td>[-1.24; -1.15]</td></tr></tbody></table>

<table><thead><tr><th>d</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>1</td><td>-2.00</td><td>-2.00</td><td>[-2.02; -1.97]</td></tr><tr><td>2</td><td>-1.50</td><td>-1.48</td><td>[-1.5; -1.45]</td></tr><tr><td>3</td><td>-1.33</td><td>-1.30</td><td>[-1.32; -1.28]</td></tr><tr><td>4</td><td>-1.25</td><td>-1.26</td><td>[-1.27; -1.24]</td></tr><tr><td>5</td><td>-1.20</td><td>-1.20</td><td>[-1.24; -1.16]</td></tr><tr><td>6</td><td>-1.17</td><td>-1.16</td><td>[-1.19; -1.12]</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: Summary of experiments in which integrals of $\iota$-dimensional functions are estimated by projecting a single $(N, 6)$-Dirichlet DPP ($\iota = 1...6$). A ● (resp. ×) indicates that the adjusted p-value of the Shapiro-Wilk test is not smaller (resp. smaller) than 5%.</figcaption></figure>

<table><thead><tr><th>l</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>6</td><td>-1.17</td><td>-1.12</td><td>[-1.15 ; -1.09]</td></tr><tr><td>5</td><td>-1.17</td><td>-1.13</td><td>[-1.15 ; -1.11]</td></tr><tr><td>4</td><td>-1.17</td><td>-1.13</td><td>[-1.16 ; -1.09]</td></tr><tr><td>3</td><td>-1.17</td><td>-1.14</td><td>[-1.16 ; -1.11]</td></tr><tr><td>2</td><td>-1.17</td><td>-1.13</td><td>[-1.16 ; -1.1]</td></tr><tr><td>1</td><td>-1.17</td><td>-1.14</td><td>[-1.17 ; -1.11]</td></tr></tbody></table>

<table><thead><tr><th>l</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>6</td><td>-1.17</td><td>-1.14</td><td>[-1.17 ; -1.1]</td></tr><tr><td>5</td><td>-1.17</td><td>-1.13</td><td>[-1.17 ; -1.1]</td></tr><tr><td>4</td><td>-1.17</td><td>-1.13</td><td>[-1.16 ; -1.09]</td></tr><tr><td>3</td><td>-1.17</td><td>-1.14</td><td>[-1.17 ; -1.11]</td></tr><tr><td>2</td><td>-1.17</td><td>-1.14</td><td>[-1.18 ; -1.1]</td></tr><tr><td>1</td><td>-1.17</td><td>-1.13</td><td>[-1.16 ; -1.1]</td></tr></tbody></table>

<table><thead><tr><th>l</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>6</td><td>-1.17</td><td>-1.19</td><td>[-1.24 ; -1.15]</td></tr><tr><td>5</td><td>-1.17</td><td>-1.20</td><td>[-1.25 ; -1.15]</td></tr><tr><td>4</td><td>-1.17</td><td>-1.20</td><td>[-1.25 ; -1.15]</td></tr><tr><td>3</td><td>-1.17</td><td>-1.20</td><td>[-1.24 ; -1.15]</td></tr><tr><td>2</td><td>-1.17</td><td>-1.19</td><td>[-1.23 ; -1.14]</td></tr><tr><td>1</td><td>-1.17</td><td>-1.17</td><td>[-1.21 ; -1.13]</td></tr></tbody></table>

<table><thead><tr><th>l</th><th>-1-1/d</th><th>Estimation</th><th>Conf. interval</th></tr></thead><tbody><tr><td>6</td><td>-1.17</td><td>-1.16</td><td>[-1.19 ; -1.12]</td></tr><tr><td>5</td><td>-1.17</td><td>-1.16</td><td>[-1.2 ; -1.13]</td></tr><tr><td>4</td><td>-1.17</td><td>-1.16</td><td>[-1.2 ; -1.12]</td></tr><tr><td>3</td><td>-1.17</td><td>-1.17</td><td>[-1.21 ; -1.13]</td></tr><tr><td>2</td><td>-1.17</td><td>-1.16</td><td>[-1.19 ; -1.12]</td></tr><tr><td>1</td><td>-1.17</td><td>-1.15</td><td>[-1.18 ; -1.11]</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 4: Summary of experiments in which integrals of function (4.4) are estimated using (N, 1)-Dirichlet DPP.</figcaption></figure>

It is clear that $h_{\gamma} \in L_{\text{per}}^{2}([0, 1])$ if $\gamma > 1/2$, and in such a case, $h_{\gamma} \in \mathcal{H}^{s}([0, 1])$ for $s < \gamma - 1/2$.

We repeat same experiments as in Section 4.1, but only for the case $d = 1$: getting an accurate evaluation of (4.4) becomes really expensive for higher dimension. We consider three values for $\gamma$: 0.625, 0.75 and 0.875. Empirical results are depicted in Figure 4. As expected, the normality is rejected for all values of $N$ (for more readability we keep the dot representation for the points). Surprisingly, the results remains very satisfactory. When $\gamma = 0.75, 0.875$, the points remain well-aligned and the estimated slopes and their confidence intervals are in good agreement with the theoretical $-1 - 2s/d$. The case $\gamma = 0.625$ is less convincing since the position of the points are more fluctuating around the regression lines and the estimated slope is further from the expected one. In its defence, $h_{0.625}$ is a very irregular function.

# 5 Conclusion

In this paper, we build a specific class of repulsive point process and use its realization as quadrature points to estimate integrals. The resulting Monte-Carlo estimator is unbiased with a variance scaling as $O(N^{-1-(2s^1)/d})$ when the integrand belongs to $\mathcal{H}^s([0, 1]^d)$. Our methodology and results have the interest to be applied to non differentiable functions, an assumption which is often considered for other methods such as grid-based methods, RQMC, scrambled nets [21]. We also show that the initial configuration of points can be used to estimate $\iota$-dimensional integrands ($\iota = 1, \dots, d$) with the same efficiency than in dimension $d$.

We consider in this paper integrals on $[0, 1]^d$ (or rectangles). Addressing the similar question for more complex domains or for improper integrals is of great interest. Similarly, considering integrands defined on manifolds such as torus or spheres is definitely an interesting perspective
In this paper we did not exploit the form of the integrand $f$. "Is it possible to improve the rate of convergence by designing an appropriate DPP which exploits the form of $f$?" seems to be a difficult but again an interesting question.

One of the limitations of the methodology lies in the simulation of a DPP. Current algorithms have a computational cost $O(N^3)$ and are not sequential. However this is a very active and productive research area. Should a faster algorithm come up, it will be directly plugged into this paper to generate $(N, d)$-DPPs. Finally, it is worth reminding the reader that we use homogeneous point patterns which can be used to estimate any integral. Therefore the 2500 replications from our model, for $N = 50, 100, ..., 1000$ and $d = 1, ..., 6$ can be used to evaluate any integral and are freely available upon request.

## References

[1] A. Baddeley. "Spatial Point Processes and their Applications". In: *Stochastic geometry: Lectures given at the C.I.M.E. Summer School held in Martina Franca, Italy, September 13–18, 2004*. Ed. by Wolfgang Weil. Berlin, Heidelberg: Springer, 2007, pp. 1–75.

[2] A. Baddeley, E. Rubak, and R. Turner. *Spatial point patterns: methodology and applications with R*. Chapman and Hall/CRC, 2015.

[3] R. Bardenet and A. Hardy. "Monte Carlo with determinantal point processes". In: *Ann. Appl. Probab.* (2020). (to appear).

[4] K. Basu and R. Mukherjee. "Asymptotic normality of scrambled geometric net quadrature". In: *The Annals of Statistics* 45.4 (2017), pp. 1759–1788.

[5] L. Debnath and P. Mikusiński. *Introduction to Hilbert spaces with applications*. 3rd. Elsevier Academic Press, 2005.

[6] L. Decreusefond et al. "Determinantal point processes". In: *Stochastic Analysis for Poisson Point Processes*. Cham: Springer International Publishing, 2016, pp. 311–342.

[7] J. Dick, F. Y. Kuo, and I. H. Sloan. "High-dimensional integration: the quasi-Monte Carlo way". In: *Acta Numerica* 22 (2013), pp. 133–288.

[8] R. Douc, E. Moulines, and D. Stoffer. *Nonlinear time series: Theory, methods and applications with R examples*. Chapman and Hall/CRC, 2014.

[9] Seymour Haber. "A modified Monte-Carlo quadrature". In: *Mathematics of Computation* 20.95 (1966), pp. 361–368.

[10] J. H. Halton. "Algorithm 247: Radical-inverse Quasi-random Point Sequence". In: *Commun. ACM* 7.12 (1964), pp. 701–702.
[11] J.B. Hough et al. "Determinantal Processes and Independence". In: *Probability Surveys* 3 (2006), pp. 206–229.

[12] J.B. Hough et al. *Zeros of Gaussian Analytic Functions and Determinantal Point Processes*. American Mathematical Soc., 2009.

[13] J. Illian et al. *Statistical Analysis and Modelling of Spatial Point Patterns*. Statistics in Practice. Chichester: Wiley, 2008.

[14] F. Lavancier, J. Møller, and E. Rubak. "Determinantal point process models and statistical inference: Extended version". In: *J. R. Stat. Soc. B* 77.4 (2015), pp. 853–877.

[15] A. Mazoyer, J.-F. Coeurjolly, and P.-O. Amblard. "Projections of determinantal point processes". (submitted). 2019.

[16] M.D. McKay, R. J. Beckman, and W.J. Conover. "A comparison of three methods for selecting values of input variables in the analysis of output from a computer code". In: *Technometrics* 21 (1979), pp. 239–245.

[17] J. Møller and R.P. Waagepetersen. *Statistical Inference and Simulation for Spatial Point Processes*. Boca Raton: Chapman and Hall/CRC, 2004.

[18] A. B. Owen. "Scrambled net variance for integrals of smooth functions". In: *The Annals of Statistics* 25.4 (1997), pp. 1541–1562.

[19] A. B. Owen. "Local antithetic sampling with scrambled nets". In: *The Annals of Statistics* 36.5 (2008), pp. 2319–2343.

[20] A.B. Owen. "Orthogonal arrays for computer experiments, integration and visualization". In: *Statisitica Sinica* 2 (1992), pp. 439–452.

[21] Art B. Owen. *Monte Carlo theory, methods and examples*. 2013.

[22] F. Riesz and B. Sz-Nagy. *Functional Analysis*. New York: Dover Publications, 1990.

[23] C.P. Robert and G. Casella. *Monte Carlo Statistical Methods*. New York: Springer-Verlag, 2004.

[24] T.J. Santner, B.J. Williams, and W.I. Notz. *The design and analysis of computer experiments*. Springer Science & Business Media, 2013.

[25] T. Shirai and Y. Takahashi. "Random point fields associated with certain Fredholm determinants I. Fermion, Poisson and boson point processes". In: *J. Funct. Anal.* 205 (2003), pp. 414–463.

[26] I.M Sobol. "On the distribution of points in a cube and the approximate evaluation of integrals". In: *USSR Computational Mathematics and Mathematical Physics* 7.4 (1967), pp. 86–112.
[27] A. Soshnikov. “Determinantal random point fields”. In: *Russ. Math. Surv+* 55.5 (2000), pp. 923–975.

[28] A. Soshnikov. “The central limit theorem for local linear statistics in classical compact groups and related combinatorial identities”. In: *The Annals of Probability* 28.3 (2000), pp. 1353–1370.

[29] A. Soshnikov. “Gaussian limit for determinantal random point fields”. In: *Annals of Probability* (2002), pp. 171–187.

[30] A. Zygmund. *Trigonometric Series (Volumes I & II combined)*. 3rd. Cambridge University Press, 2003.

# A Proofs

## A.1 Proof of Proposition 3.1

*Proof.* Since $X$ is a projection kernel, in (1.8), $\lambda_j = 0$ or 1, $\rho_Y = N$. Moreover, the trick is that the Fourier basis satisfies $\phi_j(u)\phi_k(u) = \phi_{j-k}(u)$ for any $j, k \in \mathbb{Z}^d$ and $u \in B$. These facts reduce (1.8) to

$$ \mathrm{Var}[\hat{\mu}_N(f)] = \frac{1}{N} \int_B f(u)^2 du - \sum_{j,k \in E_N} |\hat{f}(j-k)|^2, \quad (A.1) $$

whereby (3.6) is deduced using Parseval’s identity. Let us focus now on the second term of the right-hand side of (3.6). When $d = 1$, it is quite well-known that

$$ \sum_{j,k=0}^{N-1} |\hat{f}(j-k)|^2 = \sum_{|l| \le N-1} (N - |l|) |\hat{f}(l)|^2. \quad (A.2) $$

The result is easily deduced using the fact that $E_N$ is a rectangular subset of $\mathbb{Z}^d$. $\square$

## A.2 Proof of Theorem 3.1

*Proof.* Let $i = 1, \dots, d$ and $j_i \in \mathbb{Z}$ with $|j_i| \le n_i - 1$. Then,

$$ \begin{aligned} \prod_{i=1}^{d} (n_i - |j_i|) &= N \prod_{i=1}^{d} \left(1 - \frac{|j_i|}{n_i}\right) \\ &= N \sum_{k=0}^{d} (-1)^k \sigma_k \left( \frac{|j_1|}{n_1}, \dots, \frac{|j_d|}{n_d} \right) \\ &= N \left( 1 - \sum_{i=1}^{d} \frac{|j_i|}{n_i} + \sum_{1 \le i_1 < i_2 < d} \frac{|j_{i_1} j_{i_2}|}{n_{i_1} n_{i_2}} + \dots + (-1)^d \prod_{i=1}^{d} \frac{|j_i|}{n_i} \right) \end{aligned} $$
where, for $k = 1...d$, $\sigma_k(x_1, ..., x_d)$ denotes the k-th elementary symmetric polynomials in d variables:

$$ \sigma_k(x_1, ..., x_d) = \sum_{1 \le i_1 < \dots < i_k \le d} x_{i_1} \dots x_{i_k} $$

For $u, v \in \mathbb{R}^d$ (with $u_\ell \ne 0$, $\ell = 1, ..., d$), we define $1/u$ (resp. $v/u$) by $1/u = (1/u_1, ..., 1/u_d)^\top$ (resp. $v/u = (v_1/u_1, ..., v_d/u_d)^\top$). Similarly, $\|u\|_\infty$ stands for $\max_\ell |u_\ell|$. It is clear that for any $j \in F_N = \{j \in \mathbb{Z}^d : |j_i| \le n_i - 1, i = 1, ..., d\}$ and $k \ge 2$

$$ \sigma_k \left( \frac{|j_1|}{n_1}, \dots, \frac{|j_d|}{n_d} \right) = \mathcal{O} \left( \left\| \frac{j}{n} \right\|_{\infty}^2 \right). $$

Therefore, (3.7) can be rewritten as

$$ \begin{aligned} \mathrm{Var}[\hat{\mu}_N(f)] = & \frac{1}{N} \sum_{j \in F_N^c} |\hat{f}(j)|^2 + \frac{1}{N} \sum_{j \in F_N} \left\| \frac{j}{n} \right\|_1^2 |\hat{f}(j)|^2 \\ & + \frac{1}{N} \sum_{j \in F_N} R(j/n) |\hat{f}(j)|^2 \end{aligned} \quad (A.3) $$

where

$$ R(j/n) = \begin{cases} 0 & \text{if } d=1 \\ \mathcal{O}\left(\left\|\frac{j}{n}\right\|_{\infty}^{2}\right) & \text{if } d>1. \end{cases} $$

Under the asymptotic framework (3.8), we have that as $N \to \infty$

$$ \begin{aligned} \mathrm{Var}[\hat{\mu}_N(f)] \sim & \frac{1}{N} \sum_{j \in F_N^c} |\hat{f}(j)|^2 + \frac{1}{N^{1+1/d}} \sum_{j \in F_N} \left\| \frac{j}{\kappa} \right\|_1^2 |\hat{f}(j)|^2 \\ & + \frac{1}{N^{1+2/d}} \sum_{j \in F_N} R(j/\kappa) |\hat{f}(j)|^2. \end{aligned} \quad (A.4) $$

In the following, $\tau$ denotes a positive generic constant which may vary from line to line.

(i) Let $s \in (0, 1/2)$ and $j \in F_N$,

$$ \left\| \frac{j}{\kappa} \right\|_1 = \left\| \frac{j}{\kappa} \right\|_1^{2s} \left\| \frac{j}{\kappa} \right\|_1^{1-2s} \le \tau N^{\frac{1-2s}{d}} (1 + \|j\|_{\infty})^{2s}. \quad (A.5) $$

Hence,

$$ \begin{aligned} & \frac{1}{N^{1+1/d}} \sum_{j \in F_N} \left\| \frac{j}{\kappa} \right\|_1^2 |\hat{f}(j)|^2 \\ & \le \frac{\tau}{N^{1+2s/d}} \sum_{j \in F_N} (1 + \|j\|_{\infty})^{2s} |\hat{f}(j)|^2 \\ & = \mathcal{O} (N^{-1-2s/d}). \end{aligned} \quad (A.6) $$
Similarly,

$$ 
\begin{aligned}
\frac{1}{N^{1+2/d}} \sum_{j \in F_N} R(j/\kappa) |\hat{f}(j)|^2 &\le \frac{\tau N^{(2-2s)/d}}{N^{1+2/d}} \sum_{j \in F_N} (1 + \|j\|_{\infty})^{2s} |\hat{f}(j)|^2 \\
&= \mathcal{O}(N^{-1-2s/d}) \qquad (A.7)
\end{aligned}
$$

and

$$ 
\begin{aligned}
\frac{1}{N} \sum_{j \in F_N^c} |\hat{f}(j)|^2 &\le \frac{\tau}{N^{1+2s/d}} \sum_{j \in F_N^c} (1 + \|j\|_{\infty})^{2s} |\hat{f}(j)|^2 \\
&= o(N^{-1-2s/d}). \qquad (A.8)
\end{aligned}
$$

Combining (A.6)-(A.8) with (A.4) leads to the result.

(ii) We proceed similarly for the case $s > 1/2$. Let us focus on the case $d > 1$ as the other one is easily deduced. Since $f \in \mathcal{H}^s(B)$ for $s > 1/2$, we have that as $N \to \infty$

$$ \frac{1}{N^{1+1/d}} \sum_{j \in F_N} \left\| \frac{j}{\kappa} \right\|_1 |\hat{f}(j)|^2 \sim \frac{1}{N^{1+1/d}} \sum_{j \in \mathbb{Z}^d} \left\| \frac{j}{\kappa} \right\|_1 |\hat{f}(j)|^2. \qquad (A.9) $$

We also have that for any $j \in F_N$, $R(j/\kappa) \le \tau(1 + \|j\|_{\infty})^{\min(2s,2)}$ which leads to

$$ 
\begin{aligned}
\frac{1}{N^{1+2/d}} \sum_{j \in F_N} R(j/\kappa) &\le \frac{\tau}{N^{1+2/d}} \sum_{j \in F_N} (1 + \|j\|_{\infty})^{\min(2s,2)} |\hat{f}(j)|^2 \\
&= \mathcal{O}(N^{-1-2/d}). \qquad (A.10)
\end{aligned}
$$

Finally,

$$ 
\begin{aligned}
\frac{1}{N} \sum_{j \in F_N^c} |\hat{f}(j)|^2 &\le \frac{\tau}{N^{1+2s/d}} \sum_{j \in F_N^c} (1 + \|j\|_{\infty})^{2s} |\hat{f}(j)|^2 \\
&= \mathcal{O}(N^{-1-2s/d}) = o(N^{-1-1/d}). \qquad (A.11)
\end{aligned}
$$

Combining (A.9)-(A.11) with (A.4) leads again to the result.

(iii) We consider again the case $d > 1$. To achieve this step, we apply [29, Theorem 1] to the sequence of random variables $S_N(f) = \sum_{j=1}^N f(u_j)$. First, from (3.9)

$$ \mathrm{Var}(S_N(f)) \sim N^2 \mathrm{Var}(\hat{\mu}_N(f)) \sim N^{1-1/d} \sigma^2(f) \to \infty \qquad (A.12) $$

as $N \to \infty$ if $d > 1$. Second,

$$ \|f\|_{\infty} = \mathcal{O}(1) = o(N^{\tau(1-1/d)}) \qquad (A.13) $$
for any $\tau > 0$. Third, $\mathbb{E}(S_N(|f|)) = \mathcal{O}(N)$ and for $\delta = (1 + \delta')/(1 - 1/d)$ for some $\delta' > 0$, we have $\mathrm{Var}(S_N(f))^{\delta} = \mathcal{O}(N^{\delta(1-1/d)}) = \mathcal{O}(N^{1+\delta'})$, which implies that

$$ \mathbb{E}(S_N(|f|)) = \mathcal{O}(\mathrm{Var}(S_N(f))^{\delta}). \qquad (A.14) $$

Equations (A.12)-(A.14) are the key-ingredients of [29, Theorem 1], which proves that

$$ \frac{S_N(f) - N\mu(f)}{\sqrt{\mathrm{Var}(S_N(f))}} \to N(0, 1) $$

as $N \to \infty$ and yields the result. $\square$

### A.3 Proof of Corollary 3.1

*Proof.* Let $f_I^\uparrow : \mathbb{R}^d \to \mathbb{R}$ be the $d$-dimensional measurable function given by $f_I^\uparrow(x) = f_I(x_I)\mathbf{1}[x_{I^c} \in B_{I^c}]$. Then

$$ \hat{\mu}_N(f_I^\uparrow) = \frac{1}{N} \sum_{u \in \mathbf{X}} h(u) = \frac{1}{N} \sum_{u \in \mathbf{X}_I} f_I(u) = \hat{\mu}_N(f_I). $$

Since $f_I \in \mathcal{H}^s(B_I)$, it is straightforward to see that $f_I^\uparrow \in \mathcal{H}^s(B)$. In particular

$$ \begin{aligned} \sum_{j \in \mathbb{Z}^d} (1 + \|j\|_{\infty})^{2s} |\hat{f}_I^\uparrow(j)|^2 &= \sum_{j \in \mathbb{Z}^d} (1 + \|j\|_{\infty})^{2s} \mathbf{1}[j_{I^c} = \mathbf{0}_{d-\iota}] |\hat{f}_I(j_I)|^2 \\ &= \sum_{j \in \mathbb{Z}^d} (1 + \|j\|_{\infty})^{2s} |\hat{f}_I(j)|^2 \end{aligned} $$

where for some $p \ge 1$, $\mathbf{0}_p$ is the zero vector in $\mathbb{R}^p$. Corollary 3.1 is deduced by applying Theorem 3.1 to the function $f_I^\uparrow$. $\square$

### A.4 Proof of Corollary 3.2

*Proof.* Let us first note, as in the proof of Corollary 3.1 that $f_{I_\ell} \in \mathcal{H}^s(B_{I_\ell})$ is equivalent to say that $f_{I_\ell}^\uparrow \in \mathcal{H}^s(B)$. Moreover, following the proof of Theorem 3.1(i), it is shown that as $N \to \infty$

$$ N^{1+1/d} \mathrm{Cov} \left( \hat{\mu}_N(f_{I_\ell}^\uparrow), \hat{\mu}_N(f_{I_{\ell'}}^\uparrow) \right) \to (\Sigma_p)_{\ell\ell'}. $$

Now, we follow Cramèr-Wold device: let $a \in \mathbb{R}^p$, and let $Z_a = a^\top (\hat{\mu}_{N,p} - \mu_p)$. Then, $Z_a = \hat{\mu}_N(g_a) - \mu(g_a)$ with $g_a(u) = \sum_{\ell=1}^p a_\ell f_{I_\ell}^\uparrow(u)$. The result is therefore deduced since $g_a \in \mathcal{H}^s(B)$, $N^{1+1/d}\mathrm{Var}(\hat{\mu}_N(g_a)) \to a^\top \Sigma_p a$ and Theorem 3.1(ii) can be applied to $g_a$, that is as $N \to \infty$

$$ \sqrt{N^{1+1/d}} a^\top (\hat{\mu}_{N,p} - \mu_p) \to N(0, a^\top \Sigma_p a). $$

in distribution. $\square$
# B Alternative proof of Theorem 3.1

Here we propose an alternative proof of Corollary 3.1 (ii)-(iii) based on the characterization of the projected point process $X_I$.

*Proof.* (ii) From Proposition 2.1 and in particular the characterization of $(\alpha)$-DPPs as the union of independent particular DPPs (see [11]), we have that

$$ \hat{\mu}_N(f_I) = \frac{1}{N_{Ic}} \sum_{j=1}^{N_{Ic}} \hat{\mu}_{N_{I,j}}(f_I) \qquad (B.1) $$

where for $j = 1, \dots, N_{Ic}$

$$ \hat{\mu}_{N_{I,j}}(f_I) = \frac{1}{N_I} \sum_{u \in \mathbf{Y}_j} f_I(v) $$

and where $\mathbf{Y}_1, \dots, \mathbf{Y}_{N_{I,j}}$ are iid $(N_I, \iota)$-Dirichlet DPPs, that is $\hat{\mu}_{N_{I,j}}$ is nothing else than an average of unbiased estimators of $\mu(f_I)$ based on an $(N_I, \iota)$-Dirichlet DPP for which Theorem 3.1 can now be applied.

In particular, using Theorem 3.1 (ii), we have

$$ \mathrm{Var}(\hat{\mu}_N(f_I)) = \frac{1}{N_{Ic}} \mathrm{Var}(\hat{\mu}_{N_{I,1}}(f_I)) \sim \frac{1}{N_{Ic}} \frac{1}{(N_I)^{1+1/\iota}} \sigma^2(f_I) $$

as $N_I \to \infty$, with

$$ \sigma^2(f_I) = \sum_{j \in \mathbb{Z}^t} \left( \sum_{i \in I} \frac{|j_i|}{\gamma_i} \right) |\hat{f}_I(j)|^2 $$

where for any $i \in I$

$$ \gamma_i = \lim_{N_I \to \infty} n_i N_I^{-1/\iota}. $$

Now, since

$$ \gamma_i \underset{N_I \to \infty}{\sim} \kappa_i N^{1/d} N_I^{-1/\iota} \qquad (B.2) $$

where $\kappa_i$ are given by (3.8), we deduce that

$$ \mathrm{Var}(\hat{\mu}_N(f_I)) \sim \frac{1}{N^{1+1/d}} \sum_{j \in \mathbb{Z}^t} \sigma^2(f_I) $$

as $N_I \to \infty$ and where $\sigma^2(f_I)$ is given by (3.14), which yields the result.

(iii) We can observe from (B.1) that

$$ \hat{\mu}_N(f_I) - \mu(f_I) = \frac{1}{N_{Ic}} \sum_{j=1}^{N_{Ic}} (\hat{\mu}_{N_{I,j}}(f_I) - \mu(f_I)) $$
From Theorem 3.1 (ii) and for any $j \in I$

$$Z_{N_I, j} = \sqrt{N_I^{1+1/\iota} \frac{\hat{\mu}_{N_I, j}(f_I) - \mu(f_I)}{\zeta(f_I)}} \to N(0, 1)$$

in distribution. Then we apply Lindeberg-Feller theorem to establish that as $N_I^c \to \infty$

$$\frac{1}{\sqrt{N_I^c}} \sum_{j=1}^{N_I^c} Z_{N_I, j} \to N(0, 1)$$

in distribution. Therefore, in distribution as $N \to \infty$

$$\sqrt{\frac{N_I^{1+1/\iota}}{N_I^c}} \frac{\hat{\mu}_N(f_I) - \mu(f_I)}{\zeta(f_I)} \to N(0, 1)$$

whereby we deduce the result thanks to (B.2). □

<figure><img src="image_6.png" /></figure>