Nonparametric estimation of continuous DPPs
with kernel methods

Consider a simple point process $\mathcal{Y}$ on $\mathcal{X}$, that is, a random discrete subset of $\mathcal{X}$. For $\mathcal{D}\subset\mathcal{X}$, we denote by $\mathcal{Y}(\mathcal{D})$ the number of points of this process that fall within $\mathcal{D}$. Letting $m$ be a positive integer, we say that $\mathcal{X}$ has $m$-point correlation function $\varrho_{m}$ w.r.t. to the reference measure $\mu$ if, for any mutually disjoint subsets $\mathcal{D}_{1},\dots,\mathcal{D}_{m}\subset\mathcal{X}$,

In most cases, a point process is characterized by its correlation functions $(\rho_{m})_{m\geq 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. $\varrho_{m}(x_{1},\dots,x_{m})=\det[\mathsf{k}(x_{i},x_{j})]$ for all $m\geq 1$. We then say that $\mathsf{k}$ is the correlation kernel of the DPP. Not all kernels yield a DPP: if $\mathsf{k}(x,y)$ is the integral kernel of an operator $\mathsf{K}\in\mathcal{S}_{+}(L^{2}(\mathcal{X}))$, the Macchi-Soshnikov theorem [Macchi, 1975, Soshnikov, 2000] states that the corresponding DPP exists if and only if the eigenvalues of $\mathsf{K}$ are within $[0,1]$. In particular, for a finite ground set $\mathcal{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

with the likelihood operator $\mathsf{A}\in\mathcal{S}_{+}(L^{2}(\mathcal{X}))$ taken to be of the form $\mathsf{A}f(x)=\int_{\mathcal{X}}\mathsf{a}(x,y)f(y)\mathrm{d}\mu(y)$. The kernel $\mathsf{a}$ of $\mathsf{A}$ is sometimes called the likelihood kernel of the L-ensemble, to distinguish it from its correlation kernel $\mathsf{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 $x_{1},\dots,x_{m}\in\mathcal{X}$. For the rest of the paper, we assume that $\mathcal{X}\subset\mathbb{R}^{d}$ is compact, and that $\mu$ is the uniform probability measure on $\mathcal{X}$; our results straightforwardly extend to other densities w.r.t. Lebesgue. With these assumptions, the $m$-Janossy density is proportional to

where the normalization constant is a Fredholm deteminant that implicitly depends on $\mathcal{X}$. Assume now that we are given $s$ i.i.d. samples of a DPP, denoted by the sets $\mathcal{C}_{1},\dots,\mathcal{C}_{s}$ in the window $\mathcal{X}$. The associated maximum log-likelihood estimation problem reads

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 $\mathcal{S}_{+}(L^{2}(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 $\mathsf{a}(x,y)$ of the operator $\mathsf{A}$. For more details on DPPs, we refer the interested reader to Hough et al. [2006], Kulesza and Taskar [2012], Lavancier et al. [2014].

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.

In order to later apply the representer theorem of Marteau-Ferey et al. [2020], we restrict the original maximum likelihood problem (3) to “smooth” operators $\mathsf{A}=SAS^{*}$, with $A\in\mathcal{S}_{+}(\mathcal{H})$ and $S$ the restriction operator introduced in Section 1.1. Note that the kernel of $\mathsf{A}$ now writes

With this restriction on its domain, the optimization problem (3) now reads

We use a sampling approach to approximate the normalization constant in (5). We sample a set of points $\mathcal{I}=\{x^{\prime}_{i}:1\leq i\leq n\}$ i.i.d. from the ambient probability measure $\mu$. For definiteness, we place ourselves on an event happening with probability one where all the points in $\mathcal{I}$ and $\mathcal{C}_{\ell}$ for $1\leq\ell\leq s$ are distinct. We define the sample version of $f(A)$ as

where the Fredholm determinant of $\mathsf{A}=SAS^{*}$ has been replaced by the determinant of a finite matrix involving $S_{n}AS_{n}^{*}=[\langle\phi(x^{\prime}_{i}),A\phi(x^{\prime}_{j})\rangle]_{1\leq i,j\leq n}$.

The proof of Theorem 1 is given in Supplementary Material in Section C.2. Several remarks are in order. First, the high probability²²2We write $a\lesssim b$ if there exists a constant $c>0$ such that $a\leq cb$. in the statement of Theorem 1 is that of the event $\{\|S^{*}S-S^{*}_{n}S_{n}\|_{op}\lesssim c_{n}\}.$ 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 $\mathsf{T}_{k_{\mathcal{H}}}$, 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(\mathbb{I}+c_{n}A)\leq c_{n}\operatorname{Tr}(A)$, since $A\in\mathcal{S}_{+}(\mathcal{H})$. Thus, by penalizing $\operatorname{Tr}(A)$, one also improves the upper bound on the Fredholm determinant approximation error. This remark motivates the following infinite dimensional problem

for some $\lambda>0$. The penalty on $\operatorname{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 of $A$. Note that this problem depends both on the data $\mathcal{C}_{1},\dots,\mathcal{C}_{n}\subset\mathcal{X}$ and the subset $\mathcal{I}$ used for approximating the Fredholm determinant.

In an RHKS, there is a natural mapping between finite rank operators and matrices. For the sake of completeness, let $\mathbf{K}=[k_{\mathcal{H}}(z_{i},z_{j})]_{1\leq i,j\leq m}$ be a kernel matrix and let $\mathbf{K}=\mathbf{R}^{\top}\mathbf{R}$ be a Cholesky factorization. Throughout the paper, kernel matrices are always assumed to be invertible. This is not a strong assumption: if $k_{\mathcal{H}}$ is the Laplace, Gaussian or Sobolev kernel, this is true almost surely if $z_{i}$ for $1\leq i\leq 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:\mathcal{H}\to\mathbb{R}^{m}$ as $V=\sqrt{m}(\mathbf{R}^{-1})^{\top}S_{m}.$ It satisfies $VV^{*}=\mathbf{I}$, and $V^{*}V$ is the orthogonal projector onto the span of $\phi(z_{i})$ for all $1\leq i\leq m$. This is helpful to define

the finite-dimensional representative of $\phi(z_{i})\in\mathcal{H}$ for all $1\leq i\leq m$. This construction yields a useful mapping between an operator in $\mathcal{S}_{+}(\mathcal{H})$ and a finite matrix, which is instrumental for obtaining our results.

The proof of Lemma 2 is given in Section C.1. Notice that the partial isometry $V$ also helps to map a matrix in $\mathcal{S}_{+}(\mathbb{R}^{m})$ to an operator in $\mathcal{S}_{+}(\mathcal{H})$, as $\mathbf{B}\mapsto V^{*}\mathbf{B}V,$ in such a way that we have the matrix element matching $\langle\phi(z_{i}),V^{*}\mathbf{B}V\phi(z_{j})\rangle=\bm{\Phi}_{i}^{\top}\mathbf{B}\bm{\Phi}_{j}$ for all $1\leq i,j\leq m$.

The sampling approximation of the Fredholm determinant also yields a finite rank solution for (6). For simplicity, we define $\mathcal{C}\triangleq\cup_{\ell=1}^{s}\mathcal{C}_{\ell}$ and recall $\mathcal{I}=\{x^{\prime}_{1},\dots,x^{\prime}_{n}\}$. Then, write the set of points $\mathcal{Z}\triangleq\mathcal{C}\cup\mathcal{I}$ as $\{z_{1},\dots,z_{m}\}$, with $m=|\mathcal{C}|+n$, and denote the corresponding restriction operator $S_{m}:\mathcal{H}\to\mathbb{R}^{m}$. Consider the kernel matrix $\mathbf{K}=[k_{\mathcal{H}}(z_{i},z_{j})]_{1\leq i,j\leq 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

In this paper, we call $\mathbf{C}$ the representer matrix of the operator $A$. If we do the change of variables $\mathbf{B}=\mathbf{R}\mathbf{C}\mathbf{R}^{\top}$, we have the following identities: $A=mS_{m}\mathbf{C}S_{m}^{*}=V^{*}\mathbf{B}V$ and $\operatorname{Tr}(A)=\operatorname{Tr}(\mathbf{K}\mathbf{C})=\operatorname{Tr}(\mathbf{B})$, thanks to Lemma 7 in Supplementary Material. Therefore, the problem (6) boils down to the finite non-convex problem:

where $f_{n}(V^{*}\mathbf{B}V)=-\frac{1}{s}\sum_{\ell=1}^{s}\log\det\left[\bm{\Phi}_{i}^{\top}\mathbf{B}\bm{\Phi}_{j}\right]_{i,j\in\mathcal{C}_{\ell}}+\log\det\left[\delta_{ij}+\bm{\Phi}_{i}^{\top}\mathbf{B}\bm{\Phi}_{j}/|\mathcal{I}|\right]_{i,j\in\mathcal{I}}$. We assume that there is a global minimizer of (9) that we denote by $\mathbf{B}_{\star}$. The final estimator of the integral kernel of the likelihood $\mathsf{A}$ depends on $\mathbf{C}_{\star}=\mathbf{R}^{-1}\mathbf{B}_{\star}\mathbf{R}^{-1\top}$ and reads $\hat{\mathsf{a}}(x,y)=\sum_{i,j=1}^{m}\mathbf{C}_{\star ij}k_{\mathcal{H}}(z_{i},x)k_{\mathcal{H}}(z_{j},y).$ The numerical strategy is summarized in Algorithm 1.

The exact computation of the correlation kernel of the L-ensemble DPP

requires the exact diagonalization of $\mathsf{A}=SAS^{*}$. For more flexibility, we introduced a scale parameter $\gamma>0$ which often takes the value $\gamma=1$. It is instructive to approximate $\mathsf{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_{m}^{*}\mathbf{C}S_{m}$ of (6), and consider the factorization $\mathbf{C}=\mathbf{\Lambda}^{\top}\mathbf{\Lambda}$ with $\mathbf{\Lambda}=\mathbf{F}\mathbf{R}^{-1\top}$ where $\mathbf{F}^{\top}\mathbf{F}=\mathbf{B}_{\star}$ is the Cholesky factorization of $\mathbf{B}_{\star}$. Let $\{x^{\prime\prime}_{1},\dots,x^{\prime\prime}_{p}\}\subseteq\mathcal{X}$ be sampled i.i.d. from the probability measure $\mu$ and denote by $S_{p}:\mathcal{H}\to\mathbb{R}^{p}$ the corresponding restriction operator. The following integral operator

gives an approximation of $\mathsf{K}$. The numerical approach for solving (11) relies on the computation of $\mathbf{K}_{mp}=\sqrt{mp}S_{m}S^{*}_{p}=[k_{\mathcal{H}}(z_{i},x^{\prime\prime}_{j})]$ with $1\leq i\leq m$ and $1\leq j\leq p$ is a rectangular kernel matrix, associated to a fixed ordering of $\mathcal{Z}=\{z_{1},\dots,z_{m}\}$ and $\{x^{\prime\prime}_{1},\dots,x^{\prime\prime}_{p}\}$. Our strategy is described in Algorithm 2.

Next, we give a statistical guarantee for the approximation of the log-likelihood by its sample version.

The above result, proved in Section C.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 $V\mathbf{B}_{\star}V^{*}$ whose likelihood is not far from the optimal likelihood $f(A_{\star})$, as it can be shown by using a triangle inequality.

The proof of Corollary 5 is also provided in Section C.3.

An important quantity for the control of the amount of points necessary to approximate well the correlation kernel is the so-called effective dimension $d_{\rm eff}(\gamma)=\operatorname{Tr}\left(\mathsf{A}(\mathsf{A}+\gamma\mathbb{I})^{-1}\right),$ which is the expected sample size under the DPP with correlation kernel $\mathsf{K}=\mathsf{A}(\mathsf{A}+\gamma\mathbb{I})^{-1}$.

The proof of Theorem 6, given in Section C.5, mainly relies on a matrix Bernstein inequality. Let us make a few comments. First, we can simply take $\gamma=1$ in Theorem 6 to recover the common definition of the correlation kernel (1). Second, the presence of $d_{\rm eff}(\gamma)$ 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}\leq\lambda_{\max}(\mathbf{C})\lambda_{\max}(\mathbf{K})$ with $\mathbf{K}=[k_{\mathcal{H}}(z_{i},z_{j})]_{1\leq i,j\leq m}$ and where we used that $A=mS_{m}\mathbf{C}S_{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.

In Section A, we present useful technical results. Next, in a specific case, we show that the discrete problem (6) admits a closed-form solution that we discuss in Section B. Noticeably, this special case allows the understanding of the behaviour of the estimated DPP kernel in both the small and large regularization ($\lambda$) limits. In Section C, all the deferred proofs are given. Finally, Section D provides a finer analysis of the empirical results of Section 6 as well as a description of the convergence of the regularized Picard algorithm to the closed-form solution described in Section B.