Nonparametric estimation of continuous DPPs with kernel methods
Year: 2021

Abstract:
Determinantal Point Process (DPPs) are statistical models for repulsive point patterns. Both sampling and inference are tractable for DPPs, a rare feature among models with negative dependence that explains their popularity in machine learning and spatial statistics. Parametric and nonparametric inference methods have been proposed in the finite case, i.e. when the point patterns live in a finite ground set. In the continuous case, only parametric methods have been investigated, while nonparametric maximum likelihood for DPPs – an optimization problem over trace-class operators – has remained an open question. In this paper, we show that a restricted version of this maximum likelihood (MLE) problem falls within the scope of a recent representer theorem for nonnegative functions in an RKHS. This leads to a finite-dimensional problem, with strong statistical ties to the original MLE. Moreover, we propose, analyze, and demonstrate a fixed point algorithm to solve this finite-dimensional problem. Finally, we also provide a controlled estimate of the correlation kernel of the DPP, thus providing more interpretability.

heading: 1 Introduction
text: Determinantal point processes (DPPs) are a tractable family of models for repulsive point patterns, where interaction between points is parametrized by a positive semi-definite kernel. They were introduced by Macchi [1975] in the context of fermionic optics, and have gained a lot of interest since the 2000s, in particular in probability [Hough et al., 2006], spatial statistics [Lavancier et al., 2014], and machine learning [Kulesza and Taskar, 2012]. In machine learning at large, DPPs have been used essentially for two purposes: as statistical models for diverse subsets of items, like in recommendation systems [Gartrell et al., 2019], and as subsampling tools, like in experimental design [Derezinski et al., 2020], column subset selection [Belhadji et al., 2020b], or Monte Carlo integration [Gautier et al., 2019; Belhadji et al., 2019]. In this paper, we are concerned with DPPs used as statistical models for repulsion, and more specifically with inference for continuous DPPs.
DPP models in Machine Learning (ML) have so far mostly been finite DPPs: they are distributions over subsets of a (large) finite ground set, like subsets of sentences from a large corpus of documents [Kulesza and Taskar, 2012]. Since Affandi et al. [2014], a lot of effort has been put into designing efficient inference procedures for finite DPPs. In particular, the fixed point algorithm of Mariet and Sra [2015] allows for nonparametric inference of a finite DPP kernel, thus learning the features used for modelling diversity from the data. DPP models on infinite ground sets, say Rd, while mathematically and algorithmically very similar to finite DPPs, have been less popular in ML than in spatial statistics. It is thus natural that work on inference for continuous DPPs has happened mostly in the latter community; see e.g. the seminal paper [Lavancier et al., 2014]. Inference for continuous DPPs has however focused on the parametric setting, where a handful of interpretable parameters are learned. Relatedly, spatial statisticians typically learn the correlation kernel of a DPP, which is more
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
interpretable, while machine learners focus on the likelihood kernel, with structural assumptions to make learning scale to large ground sets.
In this paper, we tackle nonparametric inference for continuous DPPs using recent results on kernel methods. More precisely, maximum likelihood estimation (MLE) for continuous DPPs is an optimization problem over trace-class operators. Our first contribution is to show that a suitable modification of this problem is amenable to the representer theorem of Marteau-Ferey, Bach, and Rudi [2020]. Further drawing inspiration from the follow-up work [Rudi, Marteau-Ferey, and Bach, 2020], we derive an optimization problem over matrices, and we prove that its solution has a near optimal objective in the original MLE problem. We then propose, analyze, and demonstrate a fixed point algorithm for the resulting finite problem, in the spirit [Mariet and Sra, 2015] of nonparametric inference for finite DPPs. While our optimization pipeline focuses on the so-called likelihood kernel of a DPP, we also provide a controlled sampling approximation to its correlation kernel, thus providing more interpretability of our estimated kernel operator. A by-product contribution of independent interest is an analysis of a sampling approximation for Fredholm determinants.
The rest of the paper is organized as follows. Since the paper is notation-heavy, we first summarize our notation and give standard definitions in Section 1.1. In Section 2, we introduce DPPs and prior work on inference. In Section 3, we introduce our constrained maximum likelihood problem, and study its empirical counterpart. We analyze an algorithm to solve the latter in Section 4. Statistical guarantees are stated in Section 5, while Section 6 is devoted to numerically validating the whole pipeline. Our code is freely available1.

heading: 1.1 Notation and background
text: Sets. It is customary to define DPPs on a compact Polish space X endowed with a Radon measure µ, so that we can define the space of square integrable functions L2(X ) for this measure [Hough et al., 2009]. Outside of generalities in Section 2, we consider a compact X ⊂ Rd and µ the uniform probability measure on X . Let (H, 〈·, ·〉) be an RKHS of functions on X with a bounded continuous kernel kH(x, y), and let κ2 = supx∈X kH(x, x). Denote by φ(x) = kH(x, ·) ∈ H the canonical feature map. For a Hilbert space F , denote by S+(F) the space of symmetric and positive semi-definite trace-class operators on F . By a slight abuse of notation, we denote by S+(Rn) the space of n× n real positive semi-definite matrices. Finally, all sets are denoted by calligraphic letters (e.g. C, I).
Operators and matrices. Trace-class endomorphisms of L2(X ), seen as integral operators, are typeset as uppercase sans-serif (e.g. A,K), and the corresponding integral kernels as lowercase sans-serif (e.g. a, k). Notice that k(x, y) and kH(x, y) are distinct functions. Other operators are written in standard fonts (e.g. A,S), while we write matrices and finite-dimensional vectors in bold (e.g. K,C,v). The identity operator is written commonly as I, whereas the n× n identity matrix is denoted by In. When C is a subset of {1, . . . , n} and K is an n× n matrix, the matrix KCC is the square submatrix obtained by selecting the rows and columns of K indexed by C.
Restriction and reconstruction operators. Following Rosasco et al. [2010, Section 3], we define the restriction operator S : H → L2(X ) as (Sg)(x) = g(x). Its adjoint S∗ : L2(X ) → H is the reconstruction operator S∗h = ∫ X h(x)φ(x)dµ(x). The classical integral operator given
by TkHh = ∫ X kH(·, x)h(x)dµ(x), seen as an endomorphism of L
2(X ), thus takes the simple expression TkH = SS
∗. Similarly, the so-called covariance operator C : H → H, defined by C = ∫ X φ(x)⊗ φ(x)dµ(x), writes C = S
∗S. In the tensor product notation defining C, φ(x) is an element of the dual ofH and φ(x)⊗φ(x) is the endomorphism ofH defined by ((φ(x)⊗φ(x))(g) = g(x)φ(x); see e.g. Sterge et al. [2020]. Finally, for convenience, given a finite set {x1, . . . , xn} ⊂ X , we also define discrete restriction and reconstruction operators, respectively, as Sn : H → Rn such that Sng = (1/ √ n)[g(x1), . . . , g(xn)] >, and S∗nv = (1/ √ n) ∑n i=1 viφ(xi) for any v ∈ Rn. In particular, we have SnS∗n = (1/n)K where K = [kH(xi, xj)]1≤i,j≤n is a kernel matrix, which is defined for a given ordering of the set {x1, . . . , xn}. To avoid cumbersome expressions, when several discrete sets of different cardinalities, say n and p, are used, we simply write the respective sampling operators as Sn and Sp.
1https://github.com/mrfanuel/LearningContinuousDPPs.jl

heading: 2 Determinantal point processes and inference
text: Determinantal point processes and L-ensembles. Consider a simple point process Y on X , that is, a random discrete subset of X . For D ⊂ X , we denote by Y(D) the number of points of this process that fall within D. Letting m be a positive integer, we say that X has m-point correlation function %m w.r.t. to the reference measure µ if, for any mutually disjoint subsets D1, . . . ,Dm ⊂ X ,
E [ m∏ i=1 Y(Di) ] = ∫ ∏m i=1Di %m(x1, . . . , xm)dµ(x1) . . . dµ(xm).
In most cases, a point process is characterized by its correlation functions (ρm)m≥1. In particular, a determinantal point process (DPP) is defined as having correlation functions in the form of a determinant of a Gram matrix, i.e. %m(x1, . . . , xm) = det[k(xi, xj)] for all m ≥ 1. We then say that k is the correlation kernel of the DPP. Not all kernels yield a DPP: if k(x, y) is the integral kernel of an operator K ∈ S+(L2(X )), the Macchi-Soshnikov theorem [Macchi, 1975; Soshnikov, 2000] states that the corresponding DPP exists if and only if the eigenvalues of K are within [0, 1]. In particular, for a finite ground set X , taking the reference measure to be the counting measure leads to conditions on the kernel matrix; see Kulesza and Taskar [2012].
A particular class of DPPs is formed by the so-called L-ensembles, for which the correlation kernel writes K = A(A + I)−1, (1) with the likelihood operator A ∈ S+(L2(X )) taken to be of the form Af(x) = ∫ X a(x, y)f(y)dµ(y). The kernel a of A is sometimes called the likelihood kernel of the L-ensemble, to distinguish it from its correlation kernel k. The interest of L-ensembles is that their Janossy densities can be computed in closed form. Informally, the m-Janossy density describes the probability that the point process has cardinality m, and that the points are located around a given set of distinct points x1, . . . , xm ∈ X . For the rest of the paper, we assume that X ⊂ Rd is compact, and that µ is the uniform probability measure on X ; our results straightforwardly extend to other densities w.r.t. Lebesgue. With these assumptions, the m-Janossy density is proportional to
det(I + A)−1 · det[a(xi, xj)]1≤i,j≤m, (2)
where the normalization constant is a Fredholm deteminant that implicitly depends on X . Assume now that we are given s i.i.d. samples of a DPP, denoted by the sets C1, . . . , Cs in the window X . The associated maximum log-likelihood estimation problem reads
max A∈S+(L2(X ))
1
s s∑ `=1 log det [a(xi, xj)]i,j∈C` − log det(I + A). (3)
Solving (3) is a nontrivial problem. First, it is difficult to calculate the Fredholm determinant. Second, it is not straightforward to optimize over the space of operators S+(L2(X)) in a nonparametric setting. However, we shall see that the problem becomes tractable if we restrict the domain of (3) and impose regularity assumptions on the integral kernel a(x, y) of the operator A. For more details on DPPs, we refer the interested reader to Hough et al. [2006]; Kulesza and Taskar [2012]; Lavancier et al. [2014].
Previous work on learning DPPs. While continuous DPPs have been used in ML as sampling tools [Belhadji et al., 2019, 2020a] or models [Bardenet and Titsias, 2015; Ghosh and Rigollet, 2020], their systematic parametric estimation has been the work of spatial statisticians; see Lavancier et al. [2015]; Biscio and Lavancier [2017]; Poinas and Lavancier [2021] for general parametric estimation through (3) or so-called minimum-contrast inference. Still for the parametric case, a two-step estimation was recently proposed for non-stationary processes by Lavancier et al. [2021]. In a more general context, non-asymptotic risk bounds for estimating a DPP density are given in Baraud [2013].
Discrete DPPs have been more common in ML, and the study of their estimation has started some years ago [Affandi et al., 2014]. Unlike continuous DPPs, nonparametric estimation procedures have been investigated for finite DPPs by Mariet and Sra [2015], who proposed a fixed point algorithm. Moreover, the statistical properties of maximum likelihood estimation of discrete L-ensembles were studied by Brunel et al. [2017b]. We can also cite low-rank approaches [Dupuy and Bach, 2018;
Gartrell et al., 2017], learning with negative sampling [Mariet et al., 2019], learning with moments and cycles [Urschel et al., 2017], or learning with Wasserstein training [Anquetil et al., 2020]. Learning non-symmetric finite DPPs [Gartrell et al., 2019, 2021] has also been proposed, motivated by recommender systems.
At a high level, our paper is a continuous counterpart to the nonparametric learning of finite DPPs with symmetric kernels in Mariet and Sra [2015]. Our treatment of the continuous case is made possible by recent advances in kernel methods .

heading: 3 A sampling approximation to a constrained MLE
text: Using the machinery of kernel methods, we develop a controlled approximation of the MLE problem (3). Let us outline the main landmarks of our approach. First, we restrict the domain of the MLE problem (3) to smooth operators. On the one hand, this restriction allows us to develop a sampling approximation of the Fredholm determinant. On the other hand, the new optimization problem now admits a finite rank solution that can be obtained by solving a finite-dimensional problem. This procedure is described in Algorithm 1 and yields an estimator for the likelihood kernel. Finally, we use another sampling approximation and solve a linear system to estimate the correlation kernel of the fitted DPP; see Algorithm 2.
Restricting to smooth operators. In order to later apply the representer theorem of MarteauFerey et al. [2020], we restrict the original maximum likelihood problem (3) to “smooth” operators A = SAS∗, with A ∈ S+(H) and S the restriction operator introduced in Section 1.1. Note that the kernel of A now writes a(x, y) = 〈φ(x), Aφ(y)〉. (4) With this restriction on its domain, the optimization problem (3) now reads
min A∈S+(H) f(A) = −1 s s∑ `=1 log det [〈φ(xi), Aφ(xj)〉]i,j∈C` + log det(I + SAS ∗). (5)
Approximating the Fredholm determinant. We use a sampling approach to approximate the normalization constant in (5). We sample a set of points I = {x′i : 1 ≤ i ≤ n} i.i.d. from the ambient probability measure µ. For definiteness, we place ourselves on an event happening with probability one where all the points in I and C` for 1 ≤ ` ≤ s are distinct. We define the sample version of f(A) as
fn(A) = − 1
s s∑ `=1 log det [〈φ(xi), Aφ(xj)〉]i,j∈C` + log det(In + SnAS ∗ n),
where the Fredholm determinant of A = SAS∗ has been replaced by the determinant of a finite matrix involving SnAS∗n = [〈φ(x′i), Aφ(x′j)〉]1≤i,j≤n. Theorem 1 (Approximation of the Fredholm determinant). Let δ ∈ (0, 1/2). With probability at least 1− 2δ,
|log det(In + SnAS∗n)− log det(I + SAS∗)| ≤ log det(I + cnA),
with
cn = 4κ2 log
( 2κ2
`δ ) 3n + √ 2κ2` log ( 2κ2 `δ ) n ,
where ` = λmax(TkH) and κ 2 = supx∈X kH(x, x) <∞.
The proof of Theorem 1 is given in Supplementary Material in Section S3.2. Several remarks are in order. First, the high probability2 in the statement of Theorem 1 is that of the event {‖S∗S − S∗nSn‖op . cn}. Importantly, all the results given in what follows for the approximation of the solution of (5) only depend on this event, so that we do not need any union bound. Second, we emphasize that TkH , defined in Section 1.1, should not be confused with the correlation kernel (1). Third, to interpret the bound in Theorem 1, it is useful to recall that log det(I + cnA) ≤ cn Tr(A),
2We write a . b if there exists a constant c > 0 such that a ≤ cb.
since A ∈ S+(H). Thus, by penalizing Tr(A), one also improves the upper bound on the Fredholm determinant approximation error. This remark motivates the following infinite dimensional problem
min A∈S+(H) fn(A) + λTr(A), (6)
for some λ > 0. The penalty on Tr(A) is also intuitively needed so that the optimization problem selects a smooth solution, i.e., such a trace regularizer promotes a fast decay of eigenvalues ofA. Note that this problem depends both on the data C1, . . . , Cn ⊂ X and the subset I used for approximating the Fredholm determinant.
Finite-dimensional representatives. In an RHKS, there is a natural mapping between finite rank operators and matrices. For the sake of completeness, let K = [kH(zi, zj)]1≤i,j≤m be a kernel matrix and let K = R>R be a Cholesky factorization. Throughout the paper, kernel matrices are always assumed to be invertible. This is not a strong assumption: if kH is the Laplace, Gaussian or Sobolev kernel, this is true almost surely if zi for 1 ≤ i ≤ m are sampled e.g. w.r.t. the Lebesgue measure; see Bochner’s classical theorem [Wendland, 2004, Theorem 6.6 and Corollary 6.9]. In this case, we can define a partial isometry V : H → Rm as V = √ m(R−1)>Sm. It satisfies V V ∗ = I, and V ∗V is the orthogonal projector onto the span of φ(zi) for all 1 ≤ i ≤ m. This is helpful to define
Φi = V φ(zi) ∈ Rn, (7)
the finite-dimensional representative of φ(zi) ∈ H for all 1 ≤ i ≤ m. This construction yields a useful mapping between an operator in S+(H) and a finite matrix, which is instrumental for obtaining our results.
Lemma 2 (Finite dimensional representatives, extension of Lemma 3 in Rudi et al. [2020]). Let A ∈ S+(H). Then, the matrix B̄ = V AV ∗ is such that Φ>i B̄Φj = 〈φ(zi), Aφ(zj)〉 for all 1 ≤ i, j ≤ m, and log det(I + B̄) ≤ log det(I +A), as well as Tr(B̄) ≤ Tr(A).
The proof of Lemma 2 is given in Section S3.1. Notice that the partial isometry V also helps to map a matrix in S+(Rm) to an operator in S+(H), as B 7→ V ∗BV, in such a way that we have the matrix element matching 〈φ(zi), V ∗BV φ(zj)〉 = Φ>i BΦj for all 1 ≤ i, j ≤ m.
Finite rank solution thanks to a representer theorem. The sampling approximation of the Fredholm determinant also yields a finite rank solution for (6). For simplicity, we define C , ∪s`=1C` and recall I = {x′1, . . . , x′n}. Then, write the set of points Z , C∪I as {z1, . . . , zm}, withm = |C|+n, and denote the corresponding restriction operator Sm : H → Rm. Consider the kernel matrix K = [kH(zi, zj)]1≤i,j≤m. In particular, since we used a trace regularizer, the representer theorem of Marteau-Ferey et al. [2020, Theorem 1] holds: the optimal operator is of the form
A = m∑ i,j=1 Cijφ(zi)⊗ φ(zj) with C ∈ S+(Rm). (8)
In this paper, we call C the representer matrix of the operator A. If we do the change of variables B = RCR>, we have the following identities: A = mSmCS∗m = V
∗BV and Tr(A) = Tr(KC) = Tr(B), thanks to Lemma S1 in Supplementary Material. Therefore, the problem (6) boils down to the finite non-convex problem:
min B 0
fn(V ∗BV ) + λTr(B), (9)
where fn(V ∗BV ) = − 1s ∑s `=1 log det [ Φ>i BΦj ] i,j∈C` + log det [ δij + Φ > i BΦj/|I| ] i,j∈I . We assume that there is a global minimizer of (9) that we denote by B?. The final estimator of the integral kernel of the likelihood A depends on C? = R−1B?R−1> and reads â(x, y) =∑m i,j=1 C?ijkH(zi, x)kH(zj , y). The numerical strategy is summarized in Algorithm 1.
Algorithm 1 Estimation of the integral kernel a(x, y) of the DPP likelihood kernel A. procedure ESTIMATEA(λ, C1, . . . , Cs)
Sample I = {x′1, . . . , x′n} i.i.d. from µ . Sample n points for Fredhom det. approx. Define Z , ∪s`=1C` ∪ I . Collect all samples Compute K = R>R with K = [kH(zi, zj)]i,j∈Z . Cholesky of kernel matrix Solve (9) with iteration (14) to obtain B? . Regularized Picard iteration Compute C? = R−1B?R−1> . Representer matrix of â(x, y) return â(x, y) = ∑m i,j=1 C?ijkH(zi, x)kH(zj , y) . Likelihood kernel
end procedure
Estimation of the correlation kernel. The exact computation of the correlation kernel of the L-ensemble DPP K(γ) = A(A + γI)−1, (10) requires the exact diagonalization of A = SAS∗. For more flexibility, we introduced a scale parameter γ > 0 which often takes the value γ = 1. It is instructive to approximate K in order to easily express the correlation functions of the estimated point process. We propose here an approximation scheme based once again on sampling. Recall the form of the solution A = mS∗mCSm of (6), and consider the factorization C = Λ>Λ with Λ = FR−1> where F>F = B? is the Cholesky factorization of B?. Let {x′′1 , . . . , x′′p} ⊆ X be sampled i.i.d. from the probability measure µ and denote by Sp : H → Rp the corresponding restriction operator. The following integral operator
K̂ = mSS∗mΛ >(mΛSmS ∗ pSpS ∗ mΛ > + γIm) −1ΛSmS ∗, (11)
gives an approximation of K. The numerical approach for solving (11) relies on the computation of Kmp = √ mpSmS ∗ p = [kH(zi, x ′′ j )] with 1 ≤ i ≤ m and 1 ≤ j ≤ p is a rectangular kernel matrix, associated to a fixed ordering of Z = {z1, . . . , zm} and {x′′1 , . . . , x′′p}. Our strategy is described in Algorithm 2.
Algorithm 2 Estimation of the integral kernel k(x, y) of the DPP correlation kernel K = A(A+ I)−1.
procedure ESTIMATEK(Z,C?) Compute C? = Λ>Λ . Factorization of representer matrix Sample {x′′1 , . . . , x′′p} ⊆ X i.i.d. from µ . Sample p points Compute Kmp = [kH(zi, x′′j )] ∈ Rm×p . Cross kernel matrix Compute Ω = Λ>(ΛKmp 1pK > mpΛ > + Im) −1Λ . Representer matrix of k̂(x, y)
return k̂(x, y) = ∑m i,j=1 ΩijkH(zi, x)kH(zj , y) . Correlation kernel
end procedure

heading: 4 Implementation
text: We propose an algorithm for solving the discrete problem (9) associated to (6). To simplify the discussion and relate it to Mariet and Sra [2015], we define the objective g(X) = fn(V ∗B(X)V ) + λTr(B(X)) with the change of variables B(X) = R−1>XR−1. Then we can rephrase (9) as
min X 0 g(X) = −1 s s∑ `=1 log det(XC`C`) + log det ( I|I| + 1 n XII ) + λTr(XK−1), (12)
where we recall that n = |I|. Define for convenience U` as the matrix obtained by selecting the columns of the identity matrix which are indexed by C`, so that, we have in particular XC`C` = U>` XU`. Similarly, define a sampling matrix UI associated to the subset I. Recall the Cholesky decomposition K = R>R. To minimize (12), we start at some X0 0 and use the following iteration
Xk+1 = 1
2λ R>
(( Im + 4λR −1>p(Xk)R −1)1/2 R− Im)R, (13)
where p(X) = X + X∆X and ∆(X) = 1s ∑s `=1 U`X −1 C`C`U > ` − UI(XII + nI|I|)−1U>I . We dub this sequence a regularized Picard iteration, as it is a generalization of the Picard iteration which was introduced by Mariet and Sra [2015] in the context of learning discrete L-ensemble DPPs. In Mariet and Sra [2015], the Picard iteration, defined as Xk+1 = p(Xk), is shown to be appropriate for minimizing a different objective given by: − 1s ∑s `=1 log det(XC`C`) + log det(I + X). The following theorem indicates that the iteration (13) is a good candidate for minimizing g(X). Theorem 3. Let Xk for integer k be the sequence generated by (13) and initialized with X0 0. Then, the sequence g(Xk) is monotonically decreasing.
For a proof, we refer to Section S3.4. In practice, we use the iteration (13) with the inverse change of variables X(B) = R>BR and solve
Bk+1 = 1
2λ
( (Im + 4λq(Bk)) 1/2 − Im ) , with q(B) = B + BR∆ ( X(B) ) R>B, (14)
where ∆(X) is given hereabove. For the stopping criterion, we monitor the objective values of (9) and stop if the relative variation of two consecutive objectives is less than a predefined precision threshold tol. Contrary to (12), the objective (9) does not include the inverse of K, which might be ill-conditioned. The interplay between λ and n is best understood by considering (9) with the change of variables B′ = B/n, yielding the equivalent problem
min B′ 0 −1 s s∑ `=1 log det ( Φ>B′Φ ) C`C` + log det ( I + Φ>B′Φ ) II + λnTr(B ′),
where Φ = R is a matrix whose i-th column is Φi for 1 ≤ i ≤ m as defined in (7). Notice that, up to a 1/n factor, Φ>B′Φ is the in-sample Gram matrix of â(x, y) evaluated on the data set Z; see Algorithm 1. Thus, in the limit λ→ 0, the above expression corresponds to the MLE estimation of a finite DPP if ∪`C` ⊆ I. This is the intuitive connection with finite DPP: the continuous DPP is well approximated by a finite DPP if the ground set I is a dense enough sampling within X .

heading: 5 Theoretical guarantees
text: We now describe the guarantees coming with the approximations presented in the previous section.
Statistical guarantees for approximating the maximum likelihood problem Next, we give a statistical guarantee for the approximation of the log-likelihood by its sample version. Theorem 4 (Discrete optimal objective approximates full MLE objective). Let B? be the solution of (9). Let A? be the solution of (5). Let δ ∈ (0, 1/2). If λ ≥ 2cn(δ), then with probability at least 1− 2δ, it holds that
|f(A?)− fn(V ∗B?V )| ≤ 3
2 λTr(A?),
with 0 < cn . 1/ √ n given in Theorem 1.
The above result, proved in Section S3.3, shows that, with high probability, the optimal objective value of the discrete problem is not far from the optimal log-likelihood provided n is large enough. As a simple consequence, the discrete solution also yields a finite rank operator VB?V ∗ whose likelihood is not far from the optimal likelihood f(A?), as it can be shown by using a triangle inequality. Corollary 5 (Approximation of the full MLE optimizer by a finite rank operator). Under the assumptions of Theorem 4, if λ ≥ 2cn(δ), with probability at least 1− 2δ, it holds
|f(A?)− f(V ∗B?V )| ≤ 3λTr(A?) with cn . 1/ √ n given in Theorem 1.
The proof of Corollary 5 is also provided in Section S3.3.
Approximation of the correlation kernel An important quantity for the control of the amount of points necessary to approximate well the correlation kernel is the so-called effective dimension deff(γ) = Tr ( A(A + γI)−1 ) , which is the expected sample size under the DPP with correlation kernel K = A(A + γI)−1.
Theorem 6 (Correlation kernel approximation). Let δ ∈ (0, 1) be a failure probability, let ∈ (0, 1) be an acurracy parameter and let γ > 0 be a scale factor. Let K(γ) be the correlation kernel (10) defined with A = SAS∗. Consider K̂(γ) defined in (11) with i.i.d. sampling of p points in X wrt µ. If we take p ≥ 8κ
2‖A‖op γ 2 log ( 4deff (γ) δ‖K‖op ) , then, with probability 1− δ, the following multiplicative error
bound holds 11+ K(γ) K̂(γ) 1 1− K(γ).
The proof of Theorem 6, given in Section S3.5, mainly relies on a matrix Bernstein inequality. Let us make a few comments. First, we can simply take γ = 1 in Theorem 6 to recover the common definition of the correlation kernel (1). Second, the presence of deff(γ) in the logarithm is welcome since it is the expected subset size of the L-ensemble. Third, the quantity ‖A‖op directly influences the typical sample size to get an accurate approximation. A worst case bound is ‖A‖op ≤ λmax(C)λmax(K) with K = [kH(zi, zj)]1≤i,j≤m and where we used that A = mSmCS∗m in the light of (8). Thus, the lower bound on p may be large in practice. Although probably more costly, an approach inspired from approximate ridge leverage score sampling [Rudi et al., 2018] is likely to allow lower p’s. We leave this to future work.

heading: 6 Empirical evaluation
text: We consider an L-ensemble with correlation kernel k(x, y) = ρ exp(−‖x− y‖22/α2) defined on Rd with α = 0.05. Following Lavancier et al. [2014], this is a valid kernel if ρ < ( √ πα)−d. Note that the intensity, defined as x 7→ k(x, x), is constant equal to ρ; we shall check that the fitted kernel recovers that property. We draw samples3 from this continuous DPP in the window X = [0, 1]2. Two such samples are shown as hexagons in the first column of Figure 1, with respective intensity ρ = 50 and ρ = 100. For the estimation, we use a Gaussian kernel kH(x, y) = exp ( −‖x− y‖22/(2σ2)
) with σ > 0. The computation of the correlation kernel always uses p = 1000 uniform samples. Iteration (14) is run until the precision threshold tol = 10−5 is achieved. For stability, we add 10−10
3We used the code of Poinas and Lavancier [2021], available at https://github.com/APoinas/MLEDPP. It relies on the R package spatstat [Baddeley et al., 2015], available under the GPL-2 / GPL-3 licence.
to the diagonal of the Gram matrix K. The remaining parameter values are given in captions. We empirically observe that the regularized Picard iteration returns a matrix B? such that Φ>B?Φ is low rank; see Figure 2 (left). A lesson from Figure 1 is that the sample size of the DPP has to be large enough to retrieve a constant intensity k̂(x, x). In particular, the top row of this figure illustrates a case where σ is too small. Also, due to the large regularization λ = 0.1 and the use of only one DPP sample, the scale of ρ is clearly underestimated in this example. On the contrary, in Figure 3, for a smaller regularization parameter λ = 0.01, the intensity scale estimate is larger. We also observe that a large regularization parameter tends to smooth out the local variations of the intensity, which is not surprising. A comparison between a Gram matrix of k̂(x, y) and k(x, y) is given in Figure 2 corresponding to the example of Figure 1. The short-range diagonal structure is recovered, while some long-range structures are smoothed out. More illustrative simulations are given in Section S4, with a study of the influence of the hyperparameters, including the use of s > 1 DPP samples. In particular, the estimation of the intensity is improved if several DPP samples are used with a smaller value of λ.

heading: 7 Discussion
text: We leveraged recent progress on kernel methods to propose a nonparametric approach to learning continuous DPPs. We see three major limitations of our procedure. First, our final objective function is nonconvex, and our algorithm is only guaranteed to increase its objective function. Experimental evidence suggests that our approach recovers the synthetic kernel, but more work is needed to study the maximizers of the likelihood, in the spirit of Brunel et al. [2017a] for finite DPPs, and the properties of our fixed point algorithm. Second, the estimated integral kernel does not have any explicit structure, other than being implicitly forced to be low-rank because of the trace penalty. Adding structural assumptions might be desirable, either for modelling or learning purposes. For modelling, it is not uncommon to assume that the underlying continuous DPP is stationary, for example, which implies that the correlation kernel k(x, y) depends only on x − y. For learning, structural assumptions on the kernel may reduce the computational cost, or reduce the number of
maximizers of the likelihood. The third limitation of our pipeline is that, like most nonparametric methods, it still requires to tune a handful of hyperparameters, and, in our experience, the final performance varies significantly with the lengthscale of the RKHS kernel or the coefficient of the trace penalty. An automatic tuning procedure with guarantees would make the pipeline turn-key.
Maybe unexpectedly, future work could also deal with transferring our pipeline to the finite DPP setting. Indeed, we saw in Section 4 that in some asymptotic regime, our regularized MLE objective is close to a regularized version of the MLE objective for a finite DPP. Investigating this maximum a posteriori inference problem may shed some new light on nonparametric inference for finite DPPs. Intuitively, regularization should improve learning and prediction when data is scarce.

heading: Acknowledgements
text: We thank Arnaud Poinas and Guillaume Gautier for useful discussions on the manuscript. We acknowledge support from ERC grant BLACKJACK (ERC-2019-STG-851866) and ANR AI chair BACCARAT (ANR-20-CHIA-0002).

Review summaries:
This paper studies nonparametric MLE inference of continuous determinantal point processes (DPPs). The authors analyze various representations between discrete DPPs and continuous DPPs. In particular, they observe that the difference of normalization of DPP is bounded by the trace of a matrix. Then, they consider a new optimization problem penalizing the trace and propose a method of estimating likelihood/correlation DPP kernels. Moreover, they study an iterative method (a.k.a Picard iteration) for solving the MLE problem with the regularization term. Statistical guarantees of the MLE objective and correlation kernel are also studied. 

This paper introduces an algorithm to learn non-parametric, continuous determinantal point processes (DPPs). This algorithm relies on recent representer results for nonnegative functions in a RKHS and on a previous, related algorithm to learn discrete DPP parameterizations.   The authors prove that the proposed algorithm is guaranteed to increase the MLE objective function at each iteration, as well as prove statistical guarantees for the resulting solution: distance to an optimal parameterization, sample complexity.

This paper presents a new approach for maximum likelihood estimation (MLE) for continuous determinantal point processes (DPPs), which leverages recent work on kernel methods.  This approach involves nonparametric learning for a restricted class of the MLE problem, which is shown to fall within the scope of recent work on a representer theorem for nonnegative functions in an RKHS.  An optimization problem over matrices is characterized, and then a fixed point algorithm for solving the finite-dimensional optimization problem is proposed.  Statistical guarantees for the MLE approximation are presented.  Finally, the authors provide an empirical validation of their proposed method using experiments run on synthetic data sampled from a DPP.