Rémi Bardenet
email: remi.bardenet@gmail.com

Michalis K Titsias
email: mtitsias@aueb.gr

Inference for determinantal point processes without spectral knowledge

Determinantal point processes (DPPs) are point process models that naturally encode diversity between the points of a given realization, through a positive definite kernel K. DPPs possess desirable properties, such as exact sampling or analyticity of the moments, but learning the parameters of kernel K through likelihood-based inference is not straightforward. First, the kernel that appears in the likelihood is not K, but another kernel L related to K through an often intractable spectral decomposition. This issue is typically bypassed in machine learning by directly parametrizing the kernel L, at the price of some interpretability of the model parameters. We follow this approach here. Second, the likelihood has an intractable normalizing constant, which takes the form of a large determinant in the case of a DPP over a finite set of objects, and the form of a Fredholm determinant in the case of a DPP over a continuous domain. Our main contribution is to derive bounds on the likelihood of a DPP, both for finite and continuous domains. Unlike previous work, our bounds are cheap to evaluate since they do not rely on approximating the spectrum of a large matrix or an operator. Through usual arguments, these bounds thus yield cheap variational inference and moderately expensive exact Markov chain Monte Carlo inference methods for DPPs.

Introduction

Determinantal point processes (DPPs) are point processes [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF] that encode repulsiveness using algebraic arguments. They first appeared in [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF], and have since then received much attention, as they arise in many fields, e.g. random matrix theory, combinatorics, quantum physics. We refer the reader to [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF][START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF][START_REF] Hough | Determinantal processes and independence[END_REF] for detailed tutorial reviews, respectively aimed at audiences of machine learners, statisticians, and probabilists. More recently, DPPs have been considered as a modelling tool, see e.g. [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF][START_REF] Kulesza | Determinantal point processes for machine learning[END_REF][START_REF] Zou | Priors for diversity in generative latent variable models[END_REF]: DPPs appear to be a natural alternative to Poisson processes when realizations should exhibit repulsiveness.

In [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF], for example, DPPs are used to model diversity among summary timelines in a large news corpus. In [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], DPPs model diversity among the results of a search engine for a given query. In [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], DPPs model the spatial repartition of trees in a forest, as similar trees compete for nutrients in the ground, and thus tend to grow away from each other. With these modelling applications comes the question of learning a DPP from data, either through a parametrized form [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF][START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], or non-parametrically [START_REF] Gillenwater | Expectation-maximization for learning determinantal point processes[END_REF][START_REF] Mariet | Fixed-point algorithms for learning determinantal point processes[END_REF]. We focus in this paper on parametric inference.

Similarly to the correlation between the function values in a Gaussian process [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], the repulsiveness in a DPP is defined through a kernel K, which measures how much two points in a realization repel each other. The likelihood of a DPP involves the evaluation and the spectral decomposition of an operator L defined through a kernel L that is related to K. There are two main issues that arise when performing likelihood-based inference for a DPP. First, the likelihood involves evaluating the kernel L, while it is more natural to parametrize K instead, and there is no easy link between the parameters of these two kernels. The second issue is that the spectral decomposition of the operator L required in the likelihood evaluation is rarely available in practice, for computational or analytical reasons. For example, in the case of a large finite set of objects, as in the news corpus application [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF], evaluating the likelihood once requires the eigendecomposition of a large matrix. Similarly, in the case of a continuous domain, as for the forest application [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], the spectral decomposition of the operator L may not be analytically tractable for nontrivial choices of kernel L. In this paper, we focus on the second issue, i.e., we provide likelihood-based inference methods that assume the kernel L is parametrized, but that do not require any eigendecomposition, unlike [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF]. More specifically, our main contribution is to provide bounds on the likelihood of a DPP that do not depend on the spectral decomposition of the operator L. For the finite case, we draw inspiration from bounds used for variational inference of Gaussian processes [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], and we extend these bounds to DPPs over continuous domains.

For ease of presentation, we first consider DPPs over finite sets of objects in Section 2, and we derive bounds on the likelihood. In Section 3, we plug these bounds into known inference paradigms: variational inference and Markov chain Monte Carlo inference. In Section 4, we extend our results to the case of a DPP over a continuous domain. Readers who are only interested in the finite case, or who are unfamiliar with operator theory, can safely skip Section 4 without missing our main points. In Section 5, we experimentally validate our results, before discussing their breadth in Section 6.

DPPs over finite sets

Definition and likelihood

Consider a discrete set of items Y = {x 1 , . . . , x n }, where x i ⊂ R d is a vector of attributes that describes item i. Let K be a symmetric positive definite kernel [START_REF] Cristianini | Kernel methods for pattern recognition[END_REF] on R d , and let K = ((K(x i , x j ))) be the Gram matrix of K. The DPP of kernel K is defined as the probability distribution over all possible 2 n subsets Y ⊆ Y such that

P(A ⊂ Y ) = det(K A ), (1)
where K A denotes the sub-matrix of K indexed by the elements of A. This distribution exists and is unique if and only if the eigenvalues of K are in [0, 1] [START_REF] Hough | Determinantal processes and independence[END_REF]. Intuitively, we can think of K(x, y) as encoding the amount of negative correlation, or "repulsiveness" between x and y. Indeed, as remarked in [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF], (1) first yields that diagonal elements of K are marginal probabilities:

P(x i ∈ Y ) = K ii .
Equation (1) then entails that x i and x j are likely to co-occur in a realization of Y if and only if

det K {xi,xj } = K(x i , x i )K(y i , y i ) -K(x i , x j ) 2 = P(x i ∈ Y )P(x j ∈ Y ) -K 2 ij
is large: off-diagonal terms in K indicate whether points tend to co-occur. Providing the eigenvalues of K are further restricted to be in [0, 1), the DPP of kernel K has a likelihood [START_REF] Daley | An Introduction to the Theory of Point Processes[END_REF]. More specifically, writing Y 1 for a realization of Y ,

P(Y = Y 1 ) = det L Y1 det(L + I) , ( 2
)
where L = (I -K) -1 K, I is the n × n identity matrix, and L Y1 denotes the sub-matrix of L indexed by the elements of Y 1 . Now, given a realization Y 1 , we would like to infer the parameters of kernel K, say the parameters

θ K = (a K , σ K ) ∈ (0, ∞) 2 of a squared exponential kernel [10] K(x, y) = a K exp - x -y 2 2σ 2 K .
Since the trace of K is the expected number of points in Y [START_REF] Hough | Determinantal processes and independence[END_REF], one can estimate a K by the number of points in the data divided by n [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF]. But σ K , the parameter governing the repulsiveness, has to be fitted. If the number of items n is large, likelihood-based methods such as maximum likelihood are too costly: each evaluation of (2) requires O(n 2 ) storage and O(n 3 ) time. Furthermore, valid choices of θ K are constrained, since one needs to make sure the eigenvalues of K remain in [0, 1). A partial work-around is to note that given any symmetric positive definite kernel L, the likelihood (2) with matrix L = ((L(x i , x j ))) corresponds to a valid choice of K, since the corresponding matrix K = L(I + L) -1 necessarily has eigenvalues in [0, 1], which makes sure the DPP exists [START_REF] Hough | Determinantal processes and independence[END_REF]. The work-around consists in directly parametrizing and inferring the kernel L instead of K, so that the numerator of ( 2) is cheap to evaluate, and parameters are less constrained. Note that this step favours tractability over interpretability of the inferred parameters: if we assume

L(x, y) = a L exp - x -y 2 2σ 2 L ,
the number of points and the repulsiveness of the points in Y do not decouple as nicely as when K is squared exponential. For example, the expected number of items in Y depends on a L and σ L now, and both parameters also affect repulsiveness. There is some work investigating approximations to K to retain the more interpretable parametrization [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF], but the machine learning literature [START_REF] Kulesza | Determinantal point processes for machine learning[END_REF][START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF] almost exclusively adopts the more tractable parametrization of L. In this paper, we also make this choice of parametrizing L directly. Now, the only expensive step in the evaluation of (2) is the computation of det(L + I). While this still prevents the application of maximum likelihood, bounds on this determinant can be used in a variational approach or an MCMC algorithm, for example. In [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], bounds on det(L + I) are proposed, requiring only the first m eigenvalues of L, where m is chosen adaptively at each MCMC iteration to make the acceptance decision possible. This still requires the application of power iteration methods, which are limited to the finite domain case, require storing the whole n × n matrix L, and are prohibitively slow when the number of required eigenvalues m is large.

Nonspectral bounds on the likelihood

Let us denote by L AB the submatrix of L where row indices correspond to the elements of A, and column indices to those of B. When A = B, we simply write L A for L AA , and we drop the subscript when A = Y. Drawing inspiration from sparse approximations to Gaussian processes using inducing variables [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], we let Z = {z 1 , . . . , z m } be an arbitrary set of points in R d , and we approximate L by Q = L YZ [L Z ] -1 L ZY . Note that we do not constrain Z to belong to Y, so that our bounds do not rely on a Nyström-type approximation [START_REF] Affandi | Nyström approximation for large-scale determinantal processes[END_REF]. We term Z "pseudo-inputs", or "inducing inputs".

Proposition 1. 1 det(Q + I) e -tr(L-Q) ≤ 1 det(L + I) ≤ 1 det(Q + I) . ( 3
)
Proof. The right inequality is a straightforward consequence of the Schur complement L -Q being positive semidefinite. For instance, once could remark that for v ∈ R n ,

v T Lv ≤ 1 ⇒ v T Qv ≤ 1, so that R n 1 {v T (L+I)v≤1} dv ≤ R n 1 {v T (Q+I)v≤1} dv.
A change of variables using the Cholesky decompositions of L+I and Q+I yields the desired inequality.

The left inequality in (3) can be proved along the lines of [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], using variational arguments (as will be discussed in detail in Section 3.1.1). Next we give an alternative, more direct proof based on an inequality on determinants [14, Theorem 1]. For any real symmetric matrix A = P diag(λ i )P T , define its absolute value as |A| = P diag(|λ i |)P T . In particular, for a positive semidefinite A, |A| = A. Applying [14, Theorem 1] and noting that L, Q and L -Q are positive semidefinite, it comes

det(L + I) = det(L -Q + Q + I) ≤ det(|L -Q| + I)det(|Q| + I) = det(L -Q + I)det(Q + I) (4)
Now, denote by λi the eigenvalues of L -Q, which are all nonnegative. It comes

det(L -Q + I) = n i=1 (1 + λi ) ≤ n i=1 e λi = e tr(L-Q) , ( 5
)
where we used the inequality 1 + x ≤ e x . Plugging ( 5) into (4) yields the left part of (3).

Learning a DPP using bounds

In this section, we explain how to run variational inference and Markov chain Monte Carlo methods using the bounds in Proposition 1. In this section, we also make connections with variational sparse Gaussian processes more explicit.

Variational inference

The lower bound in Proposition 1 can be used for variational inference. Assume we have T point process realizations Y 1 , . . . , Y T , and we fit a DPP with kernel L = L θ . The log likelihood can be expressed using (2)

(θ) = T i=1 log det(L Yt ) -T log det(L + I). (6)
Let Z be an arbitrary set of m points in R d . Proposition 1 then yields a lower bound

F(θ, Z) T t=1 log det(L Yt ) -T log det(Q + I) + T tr (L -Q) ≤ (θ). (7)
The lower bound F(θ, Z) can be computed efficiently in O(nm2 ) time. Instead of maximizing (θ), we can maximize F(θ, Z) jointly with respect to the kernel parameters θ and the variational parameters Z.

To maximize [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], one can e.g. implement an EM-like scheme, alternately optimizing in Z and θ. Kernels are often differentiable with respect to θ, and sometimes F will also be differentiable with respect to the pseudo-inputs Z, so that gradient-based methods can help. In the general case, black-box optimizers such as CMA-ES [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], can also be employed.

Connections with variational inference in sparse GPs

We now discuss the connection between the variational lower bound in Proposition 1 and the variational lower bound used in sparse Gaussian process (GP; [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) models [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]. Although we do not explore this in the current paper, this connection could extend the repertoire of variational inference algorithms for DPPs by including, for instance, stochastic optimization variants.

Assume function f follows a GP distribution with zero mean function and kernel function L, so that the vector f of function values evaluated at Y follows the Gaussian distribution N (f |0, L). Then, the standard Gaussian integral yields1

1 det(L + I) = N (f |0, L)e -1 2 f T f df 2 . ( 8
)
Following this interpretation, we augment the vector f with a vector of m extra auxiliary function values u, referred to as inducing variables, evaluated at the inducing points Z so that jointly (f , u) follows

p(f , u) = N 0, L L YZ L ZY L Z , = N (f |L YZ [L Z ] -1 u, L -Q)N (u|0, L Z ). (9)
Now by using the fact that N (f |0, L) = p(f , u)du, the integral in (8) can be expanded so that

1 det(L + I) = N (f |L YZ [L Z ] -1 u, L -Q)N (u|0, L Z )e -1 2 f T f df du 2 . ( 10
)
We can bound the above integral using Jensen's inequality and the variational distribution

q(f , u) = N (f |L YZ [L Z ] -1 u, L -Q)q(u)
, where q(u) is a marginal variational distribution over the inducing variables u. This form of variational distribution is exactly the one used for sparse GPs [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], and by treating the factor q(u) optimally we can recover the left lower bound in Proposition 1, following the lines of [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]. We provide details in Appendix A.

The above connection suggests that much of the technology developed for speeding up GPs can be transferred to DPPs. For instance, if we explicitly represent the q(u) variational distribution in the above formulation, then we can develop stochastic variational inference variants for learning DPPs based on data subsampling [START_REF] Hoffman | Stochastic variational inference[END_REF]. In other words, we can apply to DPPs stochastic variational inference algorithms for sparse GPs such as [START_REF] Hensman | Gaussian processes for big data[END_REF].

Markov chain Monte Carlo inference

If approximate inference is not suitable, we can use the bounds in Proposition 1 to build a more expensive Markov chain Monte Carlo [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] sampler. Given a prior distribution p(θ) on the parameters θ of L, Bayesian inference relies on the posterior distribution π(θ) ∝ exp( (θ))p(θ), where the log likelihood (θ) is defined in [START_REF] Zou | Priors for diversity in generative latent variable models[END_REF]. A standard approach to sample approximately from π(θ) is the Metropolis-Hastings algorithm (MH; [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]Chapter 7.3]). MH consists in building an ergodic Markov chain of invariant distribution π(θ). Given a proposal q(θ |θ), the MH algorithm starts its chain at a user-defined θ 0 , then at iteration k + 1 it proposes a candidate state θ ∼ q(•|θ k ) and sets θ k+1 to θ with probability

α(θ k , θ ) = min 1, e (θ ) p(θ ) e (θ k ) p(θ k ) q(θ k |θ ) q(θ |θ k ) (11)
while θ k+1 is otherwise set to θ k . The core of the algorithm is thus to draw a Bernoulli variable with parameter α = α(θ, θ ) for θ, θ ∈ R d . This is typically implemented by drawing a uniform u ∼ U [0,1] and checking whether u < α. In our DPP application, we cannot evaluate α. But we can use Proposition 1 to build a lower and an upper bound (θ)

∈ [b -(θ, Z), b + (θ, Z)],
which can be arbitrarily refined by increasing the cardinality of Z and optimizing over Z. We can thus build a lower and upper bound for

α b -(θ , Z ) -b + (θ, Z) + log p(θ ) p(θ) ≤ log α ≤ b + (θ , Z ) -b -(θ, Z) + log p(θ ) p(θ) . ( 12
)
Now, another way to draw a Bernoulli variable with parameter α is to first draw u ∼ U [0,1] , and then refine the bounds in [START_REF] Cristianini | Kernel methods for pattern recognition[END_REF], by augmenting the numbers |Z|, |Z | of inducing variables and optimizing over Z, Z , until log u is out of the interval formed by the bounds in [START_REF] Cristianini | Kernel methods for pattern recognition[END_REF]. Then one can decide whether u < α. This Bernoulli trick is sometimes named retrospective sampling and has been suggested as early as [START_REF] Devroye | Non-uniform random variate generation[END_REF]. It has been used within MH for inference on DPPs with spectral bounds in [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], we simply adapt it to our non-spectral bounds.

The case of continuous DPPs

DPPs can be defined over very general spaces [START_REF] Hough | Determinantal processes and independence[END_REF]. We limit ourselves here to point processes on X ⊂ R d such that one can extend the notion of likelihood. In particular, we define here a DPP on X as in [1, Example 5.4(c)], by defining its Janossy density. For definitions of traces and determinants of operators, we follow [20, Section VII].

Definition

Let µ be a measure on (R d , B(R d )) that is continuous with respect to the Lebesgue measure, with density µ . Let L be a symmetric positive definite kernel. L defines a self-adjoint operator on L2 (µ) through L(f ) L(x, y)f (y)dµ(y). Assume L is trace-class, and

tr(L) = X L(x, x)dµ(x). (13)
We assume [START_REF] Affandi | Nyström approximation for large-scale determinantal processes[END_REF]

) = det((L(x i , x j )) det(I + L) µ (x 1 ) . . . µ (x n ), (14)
where B(x, r) is the open ball of center x and radius r, and where det(I + L) ∞ i=1 (λ i + 1) is the Fredholm determinant of operator L [20, Section VII]. Such a process is called the determinantal point process associated to kernel L and base measure µ. 3 Equation ( 14) is the continuous equivalent of (2). Our bounds require Ψ to be computable. This is the case for the popular Gaussian kernel with Gaussian base measure.

Nonspectral bounds on the likelihood

In this section, we derive bounds on the likelihood ( 14) that do not require to compute the Fredholm determinant det(I + L).

Proposition 2. Let

Z = {z 1 , . . . , z m } ⊂ R d , then det L Z det(L Z + Ψ) e -L(x,x)dµ(x)+tr(L -1 Z Ψ) ≤ 1 det(I + L) ≤ det L Z det(L Z + Ψ) , ( 15
)
where L Z = ((L(z i , z j )) and Ψ ij = L(z i , x)L(x, z j )dµ(x).

To see how similar [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF] is to (3), we define Q Z to be the operator on L 2 (µ) associated to the kernel

Q Z (x, x ) = L(x, Z)L -1 Z L(Z, x ), (16)
where

L(x, Z) = L(Z, x) T L(x, z 1 ) . . . L(x, z m ) .
Then the following extension of the matrix-determinant lemma shows that the common factor in the left and right hand side of ( 15) is the inverse of det(I + Q Z ), as in (3).

Lemma 4.1. With the notation of Section 4, it holds

det(I + Q Z ) = det(L Z + Ψ) det L Z .
Proof. First note that Q Z has finite rank since for f ∈ L 2 (µ),

Q Z f = M u,v=1 [L -1 Z ] uv L(z u , •) L(z v , y)f (y) dµ(y) ∈ S with S = Span (L(z i , •); 1 ≤ i ≤ M ) .
Note also that the L(z i , •)'s are linearly independent since L is a positive definite kernel. Now let (φ i ) 1≤i≤M be an orthonormal basis of S, i.e. Span(φ i ; 1 ≤ i ≤ M ) = S and

φ i φ j dµ = δ ij ,
and define the matrix W by

W ij = L(•, z i ),
(I + Q Z ) = det δ jk + Q Z φ j , φ i 1≤i,j≤n . Since Q Z φ j , φ i = M m,n=1 W nj W mi [L -1 Z ] mn it comes det(I + Q Z ) = det(I + W T L -1 Z W). Applying the classical matrix determinant lemma, it comes det(I + Q Z ) = det(L Z + WW T ) det L Z .
We finally remark that

[WW T ] ij = M k=1 L(z i , •), φ k L(z j , •), φ k = M k=1 L(z i , •), φ k φ k , L(z j , •) = L(z i , •), L(z j , •) .
Proof. (of Proposition 2) We first prove the right inequality in [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF]. From [20, Section VII.7], using [START_REF] Affandi | Nyström approximation for large-scale determinantal processes[END_REF], it holds

det(I + L) = 1 + ∞ k=1 1 k! det L(x i , x j ) dµ(x 1 ) . . . dµ(x k ). ( 17
)
We now apply the same argument as in the proof of the finite case (proof of Proposition 1). Denoting L X Y = det((L(x i , y i ))) and L X = ((det L(x i , x j )), we know from the positive definiteness of the kernel

L that L X -L X Z L -1 Z L ZX is positive semidefinite, which yields det L X ≥ det L X Z L -1 Z L ZX .
Plugging this into (17) yields the right inequality in [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF]. Upon noting that

tr(Q Z ) = [L -1 Z ] uv L(x, z u )L(z v , x)dµ(x) = tr(L -1 Z Ψ),
the proof of the left inequality in [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF] follows the lines of the proof of Proposition 1, since the main tool [14, Theorem 1] is valid for any trace-class operators.

Experiments

A toy Gaussian continuous experiment

In this section, we consider a DPP on R, so that the bounds derived in Section 4 apply. As in [7, Section 5.1], we take the base measure to be proportional to a Gaussian, i.e. its density is µ (x) = κN (x|0, (2α) -2 ). We consider a squared exponential kernel L(x, y) = exp -2 x -y 2 . In this particular case, the spectral decomposition of operator L is known [START_REF] Fasshauer | Stable evaluation of gaussian radial basis function interpolants[END_REF] 4 : the eigenfunctions of L are scaled Hermite polynomials, while the eigenvalues are a geometrically decreasing sequence. This 1D Gaussian-Gaussian example is interesting for two reasons: first, the spectral decomposition of L is known, so that we can sample exactly from the corresponding DPP [START_REF] Hough | Determinantal processes and independence[END_REF] and thus generate synthetic datasets. Second, the Fredholm determinant det(I +L) in this special case is a q-Pochhammer symbol, and can thus be efficiently computed 5 , which allows for comparison with "ideal" likelihood-based methods, to check the validity of our MCMC sampler, for instance. We emphasize that these special properties are not needed for the inference methods in Section 3, they are simply useful to demonstrate their correctness. We sample a synthetic dataset using (κ, α, ) = (1000, 0.5, 1), resulting in 13 points shown in red in Figure 1(a). Applying the variational inference method of Section 3.1, jointly optimizing in Z and θ = (κ, α, ) using the CMA-ES optimizer [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], yields poorly consistent results: κ varies over several orders of magnitude from one run to the other, and relative errors for α and go up to 100% (not shown). We thus investigate the identifiability of the parameters with the retrospective MH of Section 3.2. To limit the range of κ, we choose for (log κ, log α, log ) a wide uniform prior over We use a Gaussian proposal, the covariance matrix of which is adapted on-the-fly [START_REF] Haario | An adaptive Metropolis algorithm[END_REF] so as to reach 25% of acceptance. We start each iteration with m = 20 pseudo-inputs, and increase it by 10 and reoptimize when the acceptance decision cannot be made. Most iterations could be made with m = 20, and the maximum number of inducing inputs required in our run was 80. We show the results of a run of length 10 000 in histograms in Figures 1(b), 1(c), and 1(d). Retrospective MH and the ideal MH agree. The prior pdf is in green. The posterior marginals of α and are centered around the values used for simulation, and are very different from the prior, showing that the likelihood contains information about α and . However, as expected, almost nothing is learnt about κ, as posterior and prior roughly coincide. This is an example of the issues that come with parametrizing L directly, as mentioned in Section 1.

To conclude, we show a set of optimized pseudo-inputs Z in black in Figure 1(a). We also superimpose the marginal of any single point in the realization, which is available through the spectral decomposition of L here [START_REF] Hough | Determinantal processes and independence[END_REF]. In this particular case, this marginal is a Gaussian. Interestingly, the pseudo-inputs look like evenly spread samples from this marginal. Intuitively, they are likely to make the denominator in the likelihood ( 14) small, as they represent an ideal sample of the Gaussian-Gaussian DPP.

Diabetic neuropathy dataset

Here, we consider a real dataset of spatial patterns of nerve fibers in diabetic patients. These nerve fibers become more clustered as diabetes progresses [START_REF] Waller | Second-order spatial analysis of epidermal nerve fibers[END_REF]. The dataset consists of seven samples collected from diabetic patients at different stages of diabetic neuropathy and one healthy subject. We follow the experimental setup used in [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF] and we split the total samples into two classes: Normal/Mildly Diabetic and Moderately/Severely Diabetic. The first class contains three samples and the second one the remaining four. Figure 2 displays the point process data, which contain on average 90 points per sample in the Normal/Mildly class and 67 for the Moderately/Severely class. We wish to investigate the differences between these classes by fitting a separate DPP to each class and then quantify the differences of the repulsion or overdispersion of the point process data through the inferred kernel parameters. Paraphrasing [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], we consider a continuous DPP on R 2 , with kernel function

L(x i , x j ) = exp - 2 d=1 (x i,d -x j,d ) 2 2σ 2 d , ( 18
)
and base measure proportional to a Gaussian µ (x) = κ

2 d=1 N (x d |µ d , ρ 2 d ).
As in [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF], we quantify the overdispersion of realizations of such a Gaussian-Gaussian DPP through the quantities γ d = σ d /ρ d , which are invariant to the scaling of x. Note however that, strictly speaking, κ also mildly influences repulsion.

We investigate the ability of the variational method in Section 3.1 to perform approximate maximum likelihood training over the kernel parameters θ = (κ, σ 1 , σ 2 , ρ 1 , ρ 2 ). Specifically, we wish to fit a separate continuous DPP to each class by jointly maximizing the variational lower bound over θ and the inducing inputs Z using gradient-based optimization. Given that the number of inducing variables determines the amount of the approximation, or compression of the DPP model, we examine different settings for this number and see whether the corresponding trained models provide similar estimates for the overdispersion measures. Thus, we train the DPPs under different approximations having m ∈ {50, 100, 200, 400, 800, 1200} inducing variables and display the estimated overdispersion measures in Figures 3(a) and 3(b). These estimated measures converge to coherent values as m increases. They show a clear separation between the two classes, as also found in [START_REF] Affandi | Learning the parameters of determinantal point processes[END_REF][START_REF] Waller | Second-order spatial analysis of epidermal nerve fibers[END_REF]. Furthermore, Figures 3(c) and 3(d) show the values of the upper and lower bounds on the log likelihood, which as expected, converge to the same limit as m increases. We point out that the overall optimization of the variational lower is relatively fast in our MATLAB implementation. For instance, it takes 24 minutes for the most expensive run where m = 1200 to perform 1 000 iterations until convergence. Smaller values of m yield significantly smaller times.

Finally, as in Section 5.1, we comment on the optimized pseudo-inputs. Figure 4 displays the inducing points at the end of a converged run of variational inference for various values of m. Similarly to Figure 1(a), these pseudo-inputs are placed in remarkably neat grids and depart significantly from their initial locations.

Discussion

We have proposed novel, cheap-to-evaluate, nonspectral bounds on the determinants arising in the likelihoods of DPPs, both finite and continuous. We have shown how to use these bounds to infer the parameters of a DPP, and demonstrated their use for expensive-but-exact MCMC and cheap-butapproximate variational inference. In particular, these bounds have some degree of freedom -the pseudo-inputs -, which we optimize so as to tighten the bounds. This optimization step is crucial for likelihood-based inference of parametric DPP models, where bounds have to adapt to the point where the likelihood is evaluated to yield decisions which are consistent with the ideal underlying algorithms.

In future work, we plan to investigate connections of our bounds with the quadrature-based bounds for Fredholm determinants of [START_REF] Bornemann | On the numerical evaluation of Fredholm determinants[END_REF]. We also plan to consider variants of DPPs that condition on the number of points in the realization, to put joint priors over the within-class distributions of the features in classification problems, in a manner related to [START_REF] Zou | Priors for diversity in generative latent variable models[END_REF]. In the long term, we will investigate connections between kernels L and K that could be made without spectral knowledge, to address the issue of replacing L by K.

A On the connection to variational sparse GPs

Here, we provide further details about how the variational lower bound for DPPs over finite sets in Proposition 1 can be obtained by the variational approach to sparse GPs [START_REF] Michalis | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]. As mentioned in Section 3. A lower bound to the likelihood (2) can thus be obtained if we bound

F = log N (f |L YZ [L Z ] -1 u, L -Q)N (u|0, L Z )e -1 2 f T f df du.
This has a similar functional form with the marginal likelihood in a standard GP regression model: e -1 2 f T f plays the role of an unnormalized Gaussian likelihood where the observation vector is equal to zero and the noise variance is equal to one. To lower bound the above we can consider the variational distribution q(f , u) = N (f |L YZ [L Z ] -1 u, L -Q)q(u) and apply Jensen's inequality so that

F ≥ N (f |L YZ [L Z ] -1 u, L -Q)q(u) log N (u|0, L Z )e -1 2 f T f q(u) df du,
where the term N (f |L YZ [L Z ] -1 u, L -Q) cancels out inside the logarithm. This can be written as

F ≥ q(u) - 1 2 N (f |L YZ [L Z ] -1 u, L -Q)f T f df + log N (u|0, L Z ) q(u) du.
Further, given that

N (f |L YZ [L Z ] -1 u, L -Q)f T f df = α T α + tr(L -Q),
where α = L YZ [L Z ] -1 u, the bound can be written as

F ≥ q(u) log N (u|0, L Z )e -1 2 α T α q(u) du - 1 2 tr(L -Q).
Now if we analytically maximize w.r.t. q(u), under the constraint that q(u) is a distribution, we obtain q(u) = N (u|0, L Z )e -1 2 α T α N (u|0, L Z )e -1 2 α T α du .

Plugging this optimal q back into the bound, we obtain F ≥ log N (u|0, L Z )e -1 2 α T α du -1 2 tr(L -Q).

After computing the Gaussian integral w.r.t. u, the r.h.s. reduces to the logarithm of the DPP bound for the finite case, see Proposition 1.

B Ψ matrix for Gaussian kernels

We give here more details on the Gaussian kernel with Gaussian base measure used in the experimental Section 5. We use the notation of Section 5. .

In this Gaussian-Gaussian case, the Ψ matrix defined in Proposition 2 can be analytically computed: the ij-th element is given by

, ( 19
)
where zd = z i,d +z j,d 2

.

[

  200, 2000] × [-10, 10] × [-10, 10].

Figure 5 . 1 .

 51 Removing a burn-in sample of size 1000, we show the resulting marginal

Figure 1 :

 1 Figure 1: Results of running adaptive Metropolis-Hastings in the 1D Gaussian continuous experiment of Section 5.1. Figure 1(a) shows data in red, a set of optimized pseudo-inputs in black for θ set to the value used in the generation of the synthetic dataset, and the marginal of one point in the realization in blue. Figures 1(b), 1(c), and 1(d) show marginal histograms of κ, α, .

Figure 2 :Figure 3 :

 23 Figure 2: Six out of the seven nerve fiber samples. The first three samples (from left to right) correspond to a Normal Subject and two Mildly Diabetic Subjects, respectively. The remaining three samples correspond to a Moderately Diabetic Subject and two Severely Diabetic Subjects.

Figure 4 :

 4 Figure 4: We illustrate the optimization over the inducing inputs Z for several values of m ∈ {50, 100, 200, 400, 800, 1200} in the DPP of Section 5.2. We consider the Normal/Mildly diabetic class. The panels in the top row show the initial inducing input locations for various values of m, while the corresponding panels in the bottom row show the optimized locations.

1 . 1 , 2 .

 112 it holds that1 det(L + I) = N (f |L YZ [L Z ] -1 u, L -Q)N (u|0, L Z )e -1 2 f T f df du N (f |L YZ [L Z ] -1 u, L -Q)N (u|0, L Z )e -1 2 f T f df du.

2 .

 2 The kernel isL(x i , x j ) = e -D d=1 ( x i,d -x j,d) -µ d ) 2

[

  Ψ] ij = R D L(z i , x)L(x, z j )dµ(x) = κ i,d -z j,d ) 2 -

  to avoid technicalities. Proving[START_REF] Affandi | Nyström approximation for large-scale determinantal processes[END_REF] can be done by requiring various assumptions on L and µ. Under the assumptions of Mercer's theorem, for instance, (13) will be satisfied[20, Section VII, Theorem 2.3]. More generally, the assumptions of [21, Theorem 2.12] apply to kernels over noncompact domains, in particular the Gaussian kernel with Gaussian base measure that is often used in practice. We denote by λ i the eigenvalues of the compact operator L. There exists[1, Example 5.4(c)] a simple 2 point process on R d such that P There are n particles, one in each of the infinitesimal balls B(x i , dx i

Notice that N (f |0, L)e -1

f T f df = (2π) n/2 N (f |0, L)N (0|f , I)df = (2π) n/2 N (0|0, L + I),and (8) follows.

i.e., for which all points in a realization are distinct.

There is a notion of kernel K for general DPPs[START_REF] Hough | Determinantal processes and independence[END_REF], but we define L directly here, for the sake of simplicity. The interpretability issues of using L instead of K are the same as for the finite case, see Sections 2 and 5.

We follow the parametrization of[START_REF] Fasshauer | Stable evaluation of gaussian radial basis function interpolants[END_REF] for ease of reference.

http://docs.sympy.org/latest/modules/mpmath/functions/qfunctions.html#q-pochhammer-symbol

Acknowledgments

The authors would like to thank Adrien Hardy for useful discussions and Emily Fox for kindly providing access to the diabetic neuropathy dataset.