Determinantal point processes (DPP) have been introduced in the seventies [24] to model fermionic particles with repulsion like electrons. They recently regained interest since they represent locations of eigenvalues of some random matrices. A determinantal point process is characterized by an integral operator with kernel K and a reference measure m. The integral operator is compact and symmetric and thus characterized by its eigenfunctions and its eigenvalues. Following [18], the eigenvalues are not measurable functions of the realizations of the point process so it is difficult to devise how a modification of the eigenfunctions, respectively of the eigenvalues or of the reference measure, modifies the random configurations of a DPP.

A careful analysis of the simulation algorithm given in [18] yields several answers to these questions. For instance, it is clear that the eigenvalues control the distribution of the number of points and the eigenfunctions determine the positions of the atoms once their number is known. The above mentioned algorithm is a beautiful piece of work but requires to draw points according to distributions whose densities are not expressed as combinations of classical functions, hence the necessity to use rejection sampling method. Unfortunately, as the number of drawn points increases, the densities present high peaks and deep valleys which is the most adverse situation for rejection sampling, see Figure 1.

As a consequence, it is hardly feasible to simulate a DPP with more than one thousand points in less than a few hours. As some DPP arise as locations of eigenvalues of some matrix ensembles, it may seem faster and simpler to draw random matrices and compute their eigenvalues with the optimized libraries to do so. However, there are several drawbacks to this approach: a) we cannot control the domain into which the points fall, and for some applications it may be important to simulate DPP restricted to some compact sets, b) as eigenvalues belong to R or C, we cannot imagine DPP in higher dimensions with this approach, c) for a stationary DPP, it is often useful to simulate under the Palm measure (see below) which no longer corresponds to a known random matrix ensemble.

Several refinements of Algorithm 1 have been proposed along the years but the most advanced contributions have been made for DPP on lattices, which are of a totally different nature than continuous DPP. We here propose to speed up the simulation of a DPP by reducing the number of eigenvalues considered and approximating the eigenfunctions by functions whose quadrature can be easily inversed to get rid of the rejection part.

We evaluate the impact of these approximations by bounding the distances between the original distribution of the DPP to be simulated and the distribution according to which the points are drawn. The computations of the error margin are specific to the Ginibre DPP but could be done for other processes with radial symmetry like polyanalytic Ginibre ensembles [17] or the Bergman DPP [19].

Actually, there are several notions of distances between the distributions of point processes (see [11] and references therein). We focus here on the total variation distance and on the quadratic Wasserstein distance. The former counts the difference of the number of points in an optimal coupling between two distributions. The latter evaluates the matching distance between two realizations of an optimal coupling provided that it exists.

The paper is organized as follows. We first recall the definition and salient properties of DPP. In Section 3, we briefly introduce the optimal transportation problem in its full generality and give some elements dedicated to point processes. In Section 4, we show how the eigenvalues and eigenfunctions do appear in the evaluation of the distances under scrutiny. In Section 5, we apply these results to the simulation of the Ginibre process.

Let E be a Polish space, O(E) the family of all nonempty open subsets of E and B denotes the corresponding Borel σ-algebra. In the sequel, m is a Radon measure on (E,B) . Let NE be the space of locally finite subsets of E, also called the configuration space: NE={ξ⊂E:|Λ∩ξ|<∞for any compact setΛ⊂E}, equipped with the topology of the vague convergence – the topology generated by the seminorms pf(ξ)=∑x∈ξf(x) for f:E→R continuous with compact support. This makes NE a Polish space and we denote by FNE its Borelean σ-field (see [20]). We call elements of NE configurations and identify a locally finite configuration ξ with the atomic Radon measure ∑y∈ξδy , where we have written δy for the Dirac measure at y∈E .

Next, let NEf={ξ∈NE:ξ(E)<∞} , the space of all finite configurations on E. It is naturally equipped with the trace σ-algebra FNEf=F|NEf . Note that NEf= ⋃n=0∞NE(n) where NE(n)={ξ∈NEf,ξ(E)=n}. Since NE(n) can be identified with En/Sn where Sn is the group of permutations over n elements, every function f:NEf→R is in fact equivalent to a family of functions (fn,n≥1) where fn is a symmetric function from En to R.

A random point process μ is characterized by its Laplace transform, which is defined for any measurable nonnegative function f on E as Lμ(f)=∫NEe−∑x∈ξf(x)dμ(ξ). A point process is simple if with probability 1, for any x∈E , Ξμ({x})≤1 . In the sequel, we consider only simple point processes so that we remove the word simple. Point processes are also often characterized via their correlation functions defined as follows.

We see that ρ1 is the mean density of particles with respect to m, and ρn(x1,…,xn)dm(x1)…dm(xn) is the probability of finding at least one particle in the vicinity of each xi , i=1,…,n .

For details about the relationships between correlation functions and Janossy densities, see [7].

For details on optimal transport in Rd and in general Polish spaces, we refer to [33, 32]. For two Polish spaces X and Y, for μ (respectively ν) a probability measure on X (respectively Y), Σ(μ,ν) is the set of probability measures on X×Y whose first marginal is μ and second marginal is ν. We also need to consider a lower semicontinuous function c from X×Y to R+ . The Monge–Kantorovitch problem associated to μ, ν and c, denoted by MKP(μ , ν, c) for short, consists in finding (6) infγ∈Σ(μ,ν)∫X×Yc(x,y)dγ(x,y). More precisely, since X and Y are Polish and c is l.s.c., it is known from the general theory of optimal transportation, that there exists an optimal measure γ∈Σ(μ,ν) and that the minimum coincides with sup(F,G)∈Φc∫XFdμ+∫YGdν, where (F,G)∈Φc if and only if F∈L1(dμ) , G∈L1(dν) and F(x)+G(y)≤c(x,y) . We will denote by Wc(μ,ν) the value of the infimum in (6). In the sequel, we need the following theorem of Brenier [32, Chapter 2] or [25].

When X=Y and c is a distance on X, Wc also defines a distance on M1(Rk) , often called Kantorovitch–Rubinstein or Wasserstein-1 distance. It admits the alternative characterization. Theorem 9

Let (X,c) be a Polish space. For μ and ν, two probability measures on X, distKR(μ,ν):=Wc(μ,ν)=supf∈Lip1(X,c)∫Xfdμ−∫Xfdν, where Lip1(X,c)=f:X→R,∀x,y∈X,|f(x)−f(y)|≤c(x,y).

There are several ways to define a distance between point processes on the same underlying space E. We here focus on two of them. They are constructed similarly: Choose a cost function c on NE and then consider Wc defined by the solution of MKP(μ,ν,c) for μ and ν, two elements of M1(NE) .

In the sequel, we assume that E=Rk and the reference measure m is the Lebesgue measure on Rk . For the quadratic distance, the cost function is the squared Euclidean distance (the 1/2 factor is purely cosmetic but traditional) ce=12‖x−y‖2 and we define a cost between configurations (see also [3, 2]) as the lifting of ce on NE : For two configurations ξ1 , ξ2 , let c(ξ1,ξ2)=inf∫ce(x,y)dβ(x,y),β∈Σ(ξ1,ξ2), where Σ(ξ1,ξ2) denotes the set of β∈NE×E , i.e. configurations on E×E , having marginals ξ1 and ξ2 . First remark that when ξ1(E) is finite, the cost is finite only if ξ1(E)=ξ2(E) , otherwise Σ(ξ1,ξ2) is empty and then, by convention, the cost is infinite.

Moreover, the optimum cost is attained at the permutation of {1,…,ξ1(E)} which minimizes the sum of the squared distances c(ξ1,ξ2)=12minσ∈Sξ1(E) ∑j=1ξ1(E)‖xi−yσ(i)‖2 where ξ1={xj,1≤j≤ξ1(E)} and ξ2={yj,1≤j≤ξ1(E)} . For infinite configurations, it is not immediate that the cost function so defined has the minimum regularity required to consider an optimal transport problem. According to [29], this is indeed true as c is lower semicontinuous on NE×NE . We can then consider the Monge–Kantorovitch problem MKP(μ,ν,c) on M1(NE) . The main theorem of [8] is the following (see Definition 1 for the notations).

This means that whenever the distance between μ and ν is finite, there exists a strong coupling which works as follows: 1) draw a discrete random variable with the distribution of ξ(Λ) , and let ι be the obtained value 2) draw the points of ξ according to μι and then 3) apply the map φι to ξ. The obtained configuration is distributed according to νι .

It is shown in [8] that for two Poisson point processes of respective intensities σ1 and σ2 , the distance defined above is finite if and only if σ1(E)=σ2(E) and ζ=∑x∈ξδt(x),γopt-a.s. where t is the optimal transport map between σ1/σ1(E) and σ2/σ2(E) for the Euclidean cost as defined in Theorem 8. Note that the optimal map is a transformation which is applied to each atom independently of the others. This amounts to say that φn(x1,…,xn)=(t(x1),…,t(xn)), instead of φn(x1,…,xn)=(t1(x1,…,xn),…,tn(x1,…,xn)), in non-Poissonian cases.