On two ways to use determinantal point processes for Monte Carlo integration

arXiv:2604.19698 · cs.LG, math.ST, stat.ML, stat.TH · Submitted 2026-04-21 · Read on arXiv

Listen

Radio episode about this paper

Transcript

Introduction to the show: ident: AI Radio. Generated commentary on the latest Artificial Intelligence papers.

Tom: Today's paper: "On two ways to use determinantal point processes for Monte Carlo integration".

Jane: When approximating integrals by a weighted sum of function evaluations, determinantal point processes (DPPs) provide a method to enforce repulsion between evaluation points,

Tom: First, who's behind it and why it matters.

Title and authors: Tom: So, let's talk about who wrote this paper and what exactly they are proposing in "On two ways to use determinantal point processes for Monte Carlo integration." The authors include Guillaume Gautier, Rémi Bardenet, Michal Valko, and others from places like Inria Lille-Nord Europe and DeepMind Paris.

Jane: It’s interesting that the work is coming from such established research groups; it suggests this isn't just a casual idea but something built on a strong mathematical foundation in random matrix theory. The core concept they introduce is using determinantal point processes to control the sampling of function evaluation points for integral approximation.

Lu: The authors are revealing a close link between DPPs and the approach of Ermakov and Zolotukhin from one thousand nine hundred sixty which is a significant theoretical connection because it validates their historical intuition with modern machinery <ref:2604.19698#pg0>. They use this to show that these estimators yield unbiased estimates of Fourier-like coefficients with identical variance when projected onto kernel eigenfunctions.

Meng: So, they're linking a deterministic sampling method—the DPP—to an existing integration technique, which is helpful for understanding the underlying structure of the error we might be introducing in our computations. I’m wondering how this linkage translates into something tangible for high-dimensional data we deal with daily.

Lalam: This connection means we can potentially design sampling strategies that are inherently more robust against poor point selection, which could lead to much more stable training procedures across various complex models.

The paper's summary: Tom: Moving on to the actual summary of the paper, what they’re saying is that they focus on two specific Monte Carlo estimators derived from these DPPs: the Bardenet and Hardy estimator and the Ermakov-Zolotukhin estimator.

Jane: These are two distinct approaches, and both are shown to have strong theoretical properties. The Bardenet and Hardy estimator is associated with the multivariate Jacobi ensemble, which they show has a Central Limit Theorem with a convergence rate that is faster than classical Monte Carlo for smooth functions.

Lu: The Ermakov-Zolotukhin estimator is described as a "perfect integrator" and an "perfect interpolator of functions that are linear combinations of eigenfunctions of the associated kernel." Crucially, they find that it yields zero variance when the integrand fits within the span of the first M eigenfunctions, which is if f belongs to H N.

Meng: Zero variance sounds incredible on paper. If we can hit that condition, it means our approximation is exact for that subspace, which simplifies things immensely when we are trying to estimate expected values in high-dimensional models. What about the other estimator?

Lalam: The paper explains that solving a linear system derived from Equation (ten) with points drawn from a projection DPP yields unbiased estimates of the N Fourier-like coefficients for k=zero to N-one and these estimates are uncorrelated and share the same variance.

The paper's improvements: Tom: Now, let’s look at what improvements or specific properties they highlight. They focus on how the estimators behave in terms of their statistical properties and how they relate to standard methods like importance sampling or Quasi-Monte Carlo methods.

Jane: One key improvement is that the variance of the Ermakov-Zolotukhin estimator clearly reflects exactly how accurate the approximation of f is by its projection onto H N. The covariance between any two distinct coordinates in that system is zero, which simplifies things greatly.

Lu: Furthermore, when one eigenfunction, like phi zero is constant, the Ermakov-Zolotukhin estimator can be viewed as a quadrature rule with weights summing to mu(X), and its variance is exactly mu(X) times A <ref:2604.19698#pg2>.six. This gives us a concrete way to interpret the estimation through an existing numerical method structure.

Meng: That connection to quadrature rules suggests we can use this framework not just for Monte Carlo, but also for more structured numerical integration where we already have a known basis. It’s practical because it grounds the abstract DPP theory in something computable.

Lalam: This structure is powerful because it moves us beyond just sampling; it allows us to understand the variance as a direct measure of how well our function fits the model we're using, which is incredibly useful for debugging our AI approximations.

Conclusion: Tom: So, to wrap up this discussion on "On two ways to use determinantal point processes for Monte Carlo integration," the main points are that they established a theoretical link between DPPs and classic integrators like Ermakov and Zolotukhin, introducing two powerful estimators with distinct advantages.

Jane: They show that we have a Bardenet and Hardy estimator with good convergence properties for smooth functions, alongside the Ermakov-Zolotukhin estimator, which offers exact results under specific conditions related to function approximation in the eigenbasis.

Lu: The implications are deep; they give us tools to precisely measure the accuracy of our function approximations within certain subspaces using variance analysis, and they provide a way to sample in a controlled manner that enforces repulsion among evaluation points.

Meng: For practical application, this means we can select the right tool for the job; if our function is well-behaved in a specific subspace, we use EZ for zero variance, and if it's more general, we might lean towards BH for better convergence rates.

Lalam: I think this whole paper really points toward a future where our AI models can dynamically choose the sampling method based on the mathematical structure of the function they are trying to integrate, making our entire workflow much more efficient and reliable.

Guillaume Gautier, Rémi Bardenet, Michal Valko

Univ. Lille, CNRS, Centrale Lille

cs.LG, math.ST, stat.ML, stat.TH

Submitted: 2026-04-21

Updated: 2026-04-21

Code: https://github.com/guilgautier/DPPy

Importance score: 77/100

The gist: When approximating integrals by a weighted sum of function evaluations, determinantal point processes (DPPs) provide a method to enforce repulsion between evaluation points, which is crucial for

Key concepts

Determinantal Point Processes (DPPs)
DPPs are a mathematical framework used in Monte Carlo integration where the probability of selecting a set of points depends on the determinant of a kernel matrix. This structure naturally enforces repulsion between the chosen sampling points, which is vital for improving sampling efficiency in integration tasks.
Bardenet and Hardy (BH) Estimator
This estimator is derived from a multivariate Jacobi ensemble DPP. It is shown to have a Central Limit Theorem convergence rate faster than classical Monte Carlo for smooth functions, specifically achieving a rate proportional to p^(1+1/d) when dealing with essential C¹ functions.
Ermakov-Zolotukhin (EZ) Estimator
This estimator solves a linear system derived from Equation (10). It acts as a 'perfect integrator' for functions that are linear combinations of the kernel's eigenfunctions, meaning its variance is zero if the function belongs to a specific subspace defined by those eigenfunctions.

Terminology

Summary

When approximating integrals by a weighted sum of function evaluations, determinantal point processes (DPPs) provide a method to enforce repulsion between evaluation points, which is crucial for sampling in Monte Carlo integration. This work investigates two approaches—the Ermakov-Zolotukhin estimator and the Bardenet-Hardy estimator—using projection DPPs, focusing on their theoretical properties and empirical performance across various function classes.

The gist

The authors reveal the close link between DPPs and the approach of Ermakov and Zolotukhin [1960], providing a modern proof of their result using a generalization of the Cauchy-Binet formula that shows these estimators provide unbiased estimates of Fourier-like coefficients with identical variance, which measures the accuracy of approximation by projection onto kernel eigenfunctions.

How it works

The paper focuses on two primary Monte Carlo estimators derived from DPPs:

  1. The Bardenet and Hardy (BH) estimator, defined as:

IbBH N (f) ≜ X N n=1 f(xn) KN (xn, xn), which is associated with the multivariate Jacobi ensemble. This estimator is shown to have a Central Limit Theorem (CLT) with a faster rate than classical Monte Carlo for smooth functions, specifically proving that for essential C1 functions, the convergence rate is proportional to p(1+1/d) where p is related to N.

  1. The Ermakov-Zolotukhin (EZ) estimator, which utilizes the linear system derived from solving Equation (10):

ϕ0(x1)... ϕN-1(x1)

·

·

φ0(xN)... ϕN-1(xN)



y1... yN

)= (f(x1) · … · f(xN))

The solution vector coordinates yk are given by Equation (A.5): yk = det Φϕk-1,f (x1:N) / det Φ(x1:N). This estimator is shown to be a perfect integrator and perfect interpolator of functions that are linear combinations of eigenfunctions of the associated kernel, as it yields zero variance when the integrand is in the span of the first M eigenfunctions, i.e., if f ∈ HN.

Theoretical Connections and Properties

The authors establish several key theoretical connections between these estimators and DPP theory:

(1) Link to Fourier Coefficients:

Solving Equation (10) with points drawn from a projection DPP yields unbiased estimates of the N Fourier-like coefficients ⟨f, ϕk⟩ for k=0 to N-1. These estimates are uncorrelated and have the same variance, Var[yk] = f squared - (N X−1 l=0 ⟨f, ϕl⟩ 2).

(2) Variance Interpretation:

The variance of the EZ estimator clearly reflects the accuracy of the approximation of f by its projection onto HN. The covariance between any two distinct coordinates is zero, Cov[yj, yk] = 0, for all 1 ≤ j ≠ k ≤ N.

(3) Quadrature Rule Interpretation:

When one eigenfunction (say ϕ0) is constant, the EZ estimator IbEZ N (f) can be viewed as a quadrature rule with weights summing to µ(X), and its variance is exactly µ(X) × (A.6).

Sampling from the Multivariate Jacobi Ensemble

The paper details an efficient implementation for sampling from a specific DPP called the multivariate Jacobi ensemble, which arises when the base measure is separable, i.e., ω(x) = ω 1(x 1) × · · · × ω d(x d).

(1) Sampling Procedure:

The authors propose a two-layer rejection sampling scheme for exact sampling from the multivariate Jacobi ensemble. This involves:

a. Selecting a multi-index k = b−1(n) with n drawn uniformly at random in the set of indices.

b. Sampling from ϕk(x) 2ω(x) dx using rejection sampling with proposal distribution ωeq(x) dx, which is the limiting marginal distribution as N goes to infinity.

(2) Efficiency:

The expected total number of rejections for this procedure is of order 2 d N log N. For the one-dimensional case (d=1), a rejection-free method based on the random tridiagonal matrix model of Killip and Nenciu [2004, Theorem 2] is implemented, offering tremendous time savings.

Empirical Investigation

The paper performs extensive experiments comparing the two estimators against baseline vanilla Monte Carlo methods across various dimensions (d=1 to d=4)

Improvements for AI systems

Based on the provided scientific paper, here are specific improvements for AI systems that could be derived from its findings:

Improvement 1: Integration of Determinantal Point Process (DPP) Sampling into Bayesian Optimization and Model Selection.

The paper provides an efficient implementation for sampling from a specific DPP, the multivariate Jacobi ensemble, which is related to projection DPPs. This can be leveraged in AI systems that require sampling from complex, high-dimensional probability distributions or models where repulsion between samples is a key constraint (e.g., Gaussian Processes or Markov Chain Monte Carlo).

The improved AI system could:

  1. Implement a specialized sampler for Bayesian inference by using the exact sampling procedure detailed in Section 3.3 for the multivariate Jacobi ensemble. This offers an expected total number of rejections of order 2dN log(N), which is significantly faster than traditional rejection sampling methods for complex, high-dimensional kernels.

  2. Apply this DPP-based sampling to design more efficient acquisition functions in Bayesian Optimization by constraining the set of candidate points (nodes) based on a kernel that enforces repulsive constraints, leading to better exploration of the search space compared to i.i.d. node sampling (Importance Sampling).

Improvement 2: Development of a New Class of Unbiased Monte Carlo Estimators for High-Dimensional Integration.

The paper introduces two unbiased estimators derived from DPPs—the Bardenet and Hardy (BH) estimator and the Ermakov and Zolotukhin (EZ) estimator. The EZ estimator is particularly powerful when the integrand is a linear combination of kernel eigenfunctions, offering zero variance in that specific case.

The improved AI system could:

  1. Develop a library of Kernel-Adapted Integrators. When an AI task involves integrating functions defined in a basis related to a specific kernel (e.g., Fourier or wavelet bases), the system should automatically select and apply the EZ estimator for maximal accuracy (zero variance) instead of relying on more general methods like vanilla Monte Carlo or the BH estimator.

  2. For tasks involving functional approximation (e.g., estimating expected values of complex neural network outputs), use the EZ estimator to achieve perfect interpolation of functions within a chosen feature subspace, providing a highly accurate and interpretable estimate whose variance directly reflects the function's residual error in that subspace.

Improvement 3: Adaptive Sampling Strategy Based on Integrand Sparsity and Smoothness.

The paper explicitly states that for certain basis functions (where the integrand is sparse or has fast-decaying coefficients), the EZ estimator performs better, while for others, it can become erratic. This suggests an adaptive strategy based on prior knowledge of the function structure.

The improved AI system could:

  1. In a machine learning pipeline, analyze the structure of a loss function or target distribution using kernel methods to identify regions where the integrand (or its approximation) is sparse or has fast-decaying coefficients.

  2. Dynamically switch between estimators: use EZ when the structure suggests good behavior (e.g., when integrating linear combinations of kernel eigenfunctions) and fall back to BH or vanilla Monte Carlo otherwise, optimizing computational cost versus accuracy in real-time for integration tasks embedded within the AI workflow.

Improvement 4: Robustness Testing and Outlier Detection for Complex Estimators.

The paper provides extensive empirical results comparing BH and EZ estimators across various function classes (smooth functions, polynomials, absolute values, Heaviside functions). It notes that the ill-conditioning of the linear system in the EZ estimator can lead to outliers.

The improved AI system could:

  1. Implement a pre-estimation diagnostic routine that analyzes the condition number of the linear system required by Theorem 1 (Equation 10) before running an integration task using EZ.

  2. If ill-conditioning is detected, the system should automatically switch to a regularized version of the EZ estimator or revert to BH, mitigating the risk of erratic behavior and ensuring stability in high-dimensional or complex non-smooth function scenarios.

Sources

Related papers