Source-linked AI summary
On two ways to use determinantal point processes for Monte Carlo integration
Guillaume Gautier, Rémi Bardenet, Michal Valko
TL;DR
Numerical integration is central to many machine-learning applications, and DPP-based estimators offer different ways to select quadrature nodes. The paper revisits two such estimators, supplies proofs and sampling tools, and reports that the Ermakov–Zolotukhin estimator converges surprisingly quickly in dimensions at most two but behaves more erratically in higher dimensions.
Problem
DPP-based integration estimators have differing properties, and the paper seeks a finer understanding of their behavior and design.
Source: §1 · pp. 1–2Method
The paper revisits the Bardenet–Hardy and Ermakov–Zolotukhin estimators, extends the latter’s analysis, and provides an exact multivariate Jacobi ensemble sampler for empirical investigation.
Source: §1 · p. 2, §3.2 · p. 5, §3.3 · pp. 6–7Results
The Ermakov–Zolotukhin estimator shows surprisingly fast empirical convergence for d ≤ 2, but more erratic behavior for d ≥ 3.
Source: §5 · p. 9Takeaways & Limitations
The estimator may merit further study for integration or interpolation when the function is decomposable in the kernel eigenbasis.
Source: §5 · p. 9A central limit theorem for the Ermakov–Zolotukhin estimator appears difficult to establish because its quadrature weights depend on all sampled nodes.
Source: §3.2 · p. 6
Abstract
from arXiv · showhide
The standard Monte Carlo estimator $\widehat{I}_N^{\mathrm{MC}}$ of $\int fd\omega$ relies on independent samples from $\omega$ and has variance of order $1/N$. Replacing the samples with a determinantal point process (DPP), a repulsive distribution, makes the estimator consistent, with variance rates that depend on how the DPP is adapted to $f$ and $\omega$. We examine two existing DPP-based estimators: one by Bardenet & Hardy (2020) with a rate of $\mathcal{O}(N^{-(1+1/d)})$ for smooth $f$, but relying on a fixed DPP. The other, by Ermakov & Zolotukhin (1960), is unbiased with rate of order $1/N$, like Monte Carlo, but its DPP is tailored to $f$. We revisit these estimators, generalize them to continuous settings, and provide sampling algorithms.
1 Introduction
Numerical integration underlies many machine-learning applications, and this paper studies how determinantal point processes can supply quadrature nodes. It revisits two DPP-based estimators and develops tools to investigate their properties.
Background
- DPPs are repulsive point distributions parameterized by a kernel, and this work belongs to integration methods that favor dissimilar nodes measured by kernels.Related approaches include kernel control variates, Bayesian quadrature, herding, and biased importance sampling.Source: §1 · p. 1
Prior DPP estimators
- Bardenet and Hardy [2019] developed a more natural DPP-based estimator and used random-matrix results to establish a fast central limit theorem; both approaches depend on efficient DPP sampling.Ermakov and Zolotukhin [1960] had earlier proposed determinantal quadrature nodes followed by a linear system for an unbiased estimator.Source: §1 · p. 2
Our contributions
- The paper connects Ermakov and Zolotukhin’s [1960] estimator to projection DPPs, gives a shorter proof, and shows that its Fourier-like coefficients can have zero-variance estimates when the integrand lies in the kernel eigenfunction span.In that case, one DPP sample permits perfect interpolation by solving a linear system.Source: §1 · p. 2
- The paper provides an exact sampler for the multivariate Jacobi ensemble and uses it to investigate both estimators in previously unexplored regimes.The implementation is available in the DPPy toolbox of Gautier et al. [2019].Source: §1 · p. 2
Scope
- The paper aims to characterize the two DPP estimators rather than compare them with the broader range of numerical integration algorithms, and presents empirical findings as motivation and guidance for further theory.The authors note that establishing a central limit theorem for the Ermakov–Zolotukhin estimator may be especially difficult.Source: §1 · p. 2
2 Quadrature, DPPs, and the multivariate Jacobi ensemble
This section introduces quadrature, projection DPPs, and the multivariate Jacobi ensemble used in the paper’s experimental work. It relates classical Gauss–Jacobi quadrature to a random, multivariate construction.
2.1 Standard quadrature
- A quadrature rule approximates an integral by a weighted sum of function evaluations at chosen nodes, and its weights need not be nonnegative or sum to one.For Gauss–Jacobi quadrature on [-1, 1], nodes are zeros of an orthonormal Jacobi polynomial; this rule exactly integrates polynomials through degree 2N − 1.Source: §2.1 · p. 2
2.2 Projection DPPs
- A projection DPP uses a projection kernel built from orthonormal functions and almost surely returns a fixed number of points equal to the kernel’s rank.Projection DPPs also serve as building blocks for general DPPs, which under stated assumptions are mixtures of them.Source: §2.2 · p. 3
- Projection DPP repulsion favors configurations whose feature vectors span a larger-volume parallelotope.The paper expresses the joint distribution through the determinant of the feature matrix.Source: §2.2 · p. 3
2.3 The multivariate Jacobi ensemble
- Multivariate orthonormal polynomials are formed by Gram–Schmidt, with multi-indices ordered first by maximum degree and then lexicographically.The ordering is needed because multivariate monomials have no natural order.Source: §2.3 · pp. 3–4
- For the multivariate Jacobi ensemble, the paper uses orthonormal polynomials on a separable base measure; on [-1, 1]^d with Jacobi weights, these define the ensemble.The univariate weights have the form (1 − z)^a_i(1 + z)^b_i, with a_i, b_i > −1.Source: §2.3 · p. 4
3 Monte Carlo integration with projection DPPs
The section presents two unbiased DPP-based Monte Carlo estimators and an exact sampling procedure for the multivariate Jacobi ensemble. The EZ estimator can perfectly integrate and interpolate functions that are linear combinations of the kernel eigenfunctions, while the BH estimator has faster-than-classical-Monte-Carlo decay under its stated smoothness condition.
3.1 A natural estimator
- For essentially C1 functions, the BH estimator has faster-than-classical-Monte-Carlo decay, with asymptotic variance governed by the decay of f’s Fourier coefficients.The multivariate Jacobi ensemble is described as a randomized multivariate counterpart of Gauss-Jacobi quadrature.Source: §3.1 · p. 4
3.2 The Ermakov-Zolotukhin estimator
- The EZ estimator links projection DPP sampling to an unbiased linear-system estimate of the integrand’s Fourier-like coefficients.The paper connects the estimator to projection DPPs and gives a simpler proof using a generalization of the Cauchy-Binet formula.Source: §3.2 · p. 5
- For f in H_N, the EZ coefficient estimates have zero variance, so one DPP sample perfectly reconstructs f.The estimates are uncorrelated and share a variance determined by the residual coefficients beyond N.Source: §3.2 · p. 5
3.3 How to sample from projection DPPs and the multivariate Jacobi ensemble
- The multivariate Jacobi sampler uses a two-layer rejection scheme with expected rejections of order 2^dN log N; sampling a 2D ensemble of 1,000 points takes under a minute rather than hours.For the univariate Jacobi ensemble, the paper also implements exact rejection-free sampling in O(N^2) using a random tridiagonal matrix.Source: §3.3 · p. 7
4 Empirical investigation
Across experiments with the multivariate Jacobi ensemble, EZ’s performance depends on dimension and how well the integrand aligns with the kernel eigenbasis. It can have very fast variance decay for eigenfunction sums, while BH is more regular in higher dimensions and Monte Carlo can compete when both DPP estimators are less effective.
4 Empirical investigation
- The experiments compare BH (7) and EZ (14) with vanilla Monte Carlo, build both estimators from the multivariate Jacobi ensemble, extend the smooth-bump CLT experiment to larger N, and examine eigenbasis-decomposable polynomials and nonsmooth functions.Source: §4 · p. 7, §4.1 · pp. 7–8
4.1 The bump experiment
- EZ dramatically outperforms BH for d ≤2 with surprisingly fast empirical convergence, but declines for d ≥3, where BH has faster and more regular variance decay.At N = 300, EZ Kolmogorov-Smirnov tests yielded small p-values across dimensions, from 0.03 to 0.24.Source: §4.1 · p. 8
- EZ estimates can have a few outliers, possibly from ill-conditioning, and are not guaranteed to preserve the sign of constant-sign integrands.The source of the outliers is left for future work.Source: §4.1 · p. 8
4.2 Integrating sums of eigenfunctions
- For each dimension, variance drops to zero when the DPP point count reaches N = M; when M = N + 1, it never reaches zero and decays at order 1/N^2.Source: §4.2 · p. 8
4.3 Further experiments
- For discontinuous or out-of-span functions, BH variance decay remains smooth but slower, EZ performance deteriorates with dimension, and vanilla Monte Carlo can compete at small N.For EZ, cross terms can introduce large-coefficient monomials outside H_N; more points reduce variance, while adapting the kernel to a sparse or fast-decaying eigenbasis is suggested.Source: §4.3 · p. 8
5 Conclusion
The paper characterizes EZ as an unbiased projection-DPP estimator whose linear system yields interpretable variance, and contrasts it with BH’s stronger theoretical properties for the multivariate Jacobi ensemble. Experiments suggest EZ is most promising when the integrand aligns with the kernel eigenbasis, while a generalized Cauchy-Binet formula supplies a modern proof argument.
5 Conclusion
- EZ’s linear system gives unbiased estimates of the integrand’s Fourier-like coefficients with identical variance, measuring approximation error from projection onto the kernel eigenfunctions.The paper also notes that BH has a less interpretable variance but stronger theoretical properties with the multivariate Jacobi ensemble.Source: §5 · pp. 8–9
- For the multivariate Jacobi ensemble, EZ converges surprisingly fast empirically for d ≤2 but behaves more erratically for d ≥3; ill-conditioning may cause outliers.Regularization may help but introduces a stability-bias trade-off; EZ is suggested for integration or interpolation when the function decomposes in the kernel eigenbasis.Source: §5 · p. 9
A.1 The generalized Cauchy-Binet formula: a modern argument
- The generalized Cauchy-Binet formula establishes that the joint distribution is normalized and supports a shorter proof of an improved Ermakov-Zolotukhin result.The paper presents it as a modern argument for Theorem 1.Source: §A.1 · p. 13
A.2 Proof of Theorem 1
The appendix proves properties of the Ermakov–Zolotukhin estimator and develops a sampling scheme for the multivariate Jacobi ensemble. It also clarifies the estimator’s quadrature interpretation and the sampling scheme’s computational trade-offs.
A.2 Proof of Theorem 1
- The EZ linear system has an almost-surely unique solution, whose coordinates estimate the corresponding Fourier coefficients with variance determined by the residual norm.The coordinate estimates are pairwise uncorrelated.Source: §A.2 · p. 13
- The coordinate estimates have zero covariance for every distinct pair.This improves the original result.Source: §A.2 · p. 13
A.3 The EZ estimator as a quadrature rule
- With a constant first eigenfunction, the EZ estimator is unbiased with known variance and forms a random quadrature rule whose weights sum to µ(X).The weights are not guaranteed to be nonnegative a priori.Source: §A.3 · pp. 14–15
A.4 Sampling from the multivariate Jacobi ensemble
- Sampling N points in dimension d requires an expected total of order 2dN log(N) rejections.The chain-rule sampler generates the points in sequence, after which their selection order can be forgotten.Source: §A.4 · p. 16
- The geometric sampling formulation avoids updating K^-1, is fully vectorized, and resembles sequential Gram–Schmidt orthogonalization.The authors present it as simpler and more interpretable than the alternative formulation.Source: §A.4 · p. 16
- Using the marginal distribution as the conditional proposal reduces costly acceptance-ratio evaluations.The marginal distribution is itself sampled using the arcsine proposal and rejection constant used by Bardenet and Hardy [2019].Source: §A.4 · pp. 16–17
A.4.1 Generate samples from the marginal distribution
The marginal sampler uses a mixture-based rejection procedure, while the timing discussion compares its computational costs across dimensions. Greater dimension shortens recurrence depth but increases rejection overhead.
A.4.1 Generate samples from the marginal distribution
- The marginal density is sampled by uniformly selecting a mixture component, then using rejection sampling with the arcsine proposal.The mixture has N components, each assigned weight 1/N.Source: §A.4.1 · p. 17
- The expected rejection count for the marginal-sampling step is of order 2^d, yielding total chain-rule rejections of order 2dN log(N).The rejection-bound argument draws on bounds for orthonormal polynomials under the stated Jacobi-parameter conditions.Source: §A.4.1 · p. 18
A.4.2 Empirical timing and number of rejections
- Greater dimension reduces recurrence depth but raises rejection rates, so rejection overhead dominates and sampling is not cheaper in larger dimensions.In d = 1, the rejection-free tridiagonal model provides substantial time savings.Source: §A.4.2 · p. 18, Figure A.2 · p. 18
- Using the marginal distribution as a proposal gives rejection constant N/(N−(n−1)) for the n-th conditional and applies in the general case.Sampling from this marginal still requires a separate procedure; when eigenfunctions are available, it can be treated as a mixture.Source: §A.4.2 · p. 18
B.1 Reproducing the bump example
The bump-function experiments contrast the estimators’ variance behavior and examine whether EZ estimates appear Gaussian. BH shows the expected decay and Gaussian fluctuations, while EZ can have near-zero variance when the bump is well approximated by a polynomial.
B.1 Reproducing the bump example
- BH variance decays at 1/N^(1+1/d), while Monte Carlo decays at 1/N and may compete at small N as d increases before BH becomes lower-variance.The figure describes a meeting point after which BH variance is smaller.Source: Figure B.1 · p. 19
- For d = 1, EZ has almost no variance for N ≥ 100 because the bump is extremely well approximated by a polynomial of degree N.The experiment uses ε = 0.05.Source: Figure B.1 · p. 19
- At N = 300, BH fluctuations are consistent with the Gaussian CLT, while outliers in EZ estimates appear to challenge Gaussianity.The bump is extremely well approximated by a degree-300 polynomial, leaving EZ with almost no variance.Source: §B.1 · p. 20
- The EZ estimator does not preserve the sign of the integrand.Source: §B.1 · p. 20
B.2 Integrating sums of eigenfunctions
The eigenfunction experiments test two regimes: a finite sum whose size is within the DPP sample count, and a sum with one additional eigenfunction. EZ’s variance behavior differs across these settings, as do the comparisons with BH and Monte Carlo.
B.2 Integrating sums of eigenfunctions
- EZ has zero variance once N ≥ 70 for a sum of 70 kernel eigenfunctions, a performance neither BH nor vanilla Monte Carlo reaches.This is the setting where the number of eigenfunctions does not exceed the number of DPP sample points.Source: Figure B.3 · p. 21, §B.2 · p. 21
- When the sum has N + 1 eigenfunctions, BH and vanilla Monte Carlo are comparable, while EZ variance decays as 1/N^2 and eventually falls below Monte Carlo.The BH CLT conditions are not met, so its 1/N^(1+1/d) variance-decay guarantee does not apply in this setting.Source: §B.2 · pp. 21–22, Figure B.4 · p. 22
B.3 Integrating absolute value
The absolute-value experiment compares the BH and EZ estimators with vanilla Monte Carlo. EZ performs well in one dimension but becomes more erratic as dimension increases, while BH performs smoothly for d ≤ 2.
B.3 Integrating absolute value
- The comparison uses BH (7) and EZ (14) estimators with a_i and b_i drawn i.i.d. uniformly from [−1/2, 1/2].Source: §B.3 · p. 23
- For d ≤ 2, BH is smoother than vanilla Monte Carlo, and its 1/N^(1+1/d) rate appears to hold despite unmet CLT conditions.For larger dimensions, that rate no longer appears to hold.Source: §B.3 · p. 23, Figure B.5 · p. 23
- In d = 1, the absolute value is well approximated by a low-order truncated Taylor series, and EZ performs very well.EZ becomes more erratic as dimension increases.Source: §B.3 · p. 23
B.4 Integrating Heaviside
The Heaviside experiment compares BH and EZ on a discontinuous function. EZ is erratic, while BH retains a smooth variance decay and the expected rate in the tested dimensions despite unmet CLT conditions.
B.4 Integrating Heaviside
- The experiment applies BH and EZ using independently uniform coefficients a_i and b_i from [−1/2, 1/2].Figure B.6 is identified as a comparison of the estimators.Source: §B.4 · p. 24, Figure B.6 · p. 24
- EZ behaves erratically on the introduced discontinuity, consistent with large squared Taylor-series coefficients as the smoothing parameter tends to zero.The text suggests performance might improve for very large N.Source: §B.4 · p. 24
- Despite unmet CLT conditions, BH has smooth performance and its 1/N^(1+1/d) rate appears to hold.The comparison is reported in Figure B.6.Source: §B.4 · p. 24, Figure B.6 · p. 24
B.5 Integrating cosine
For cosine, EZ performs well through dimension two but deteriorates from dimension three, while BH and vanilla Monte Carlo behave similarly there. The next subsection tests a mixture of smooth and nonsmooth functions.
B.5 Integrating cosine
- For d ≥ 3, BH and vanilla Monte Carlo behave similarly; for d ≤ 2, BH's 1/N^(1+1/d) rate appears to hold despite unmet CLT conditions.Figure B.7 compares the estimators for cosine.Source: §B.5 · p. 25, Figure B.7 · p. 25
- For cosine, EZ performs well for d ≤ 2 but deteriorates for d ≥ 3 as Taylor cross terms introduce large-coefficient monomials outside H_N.The text suggests performance might improve for very large N.Source: §B.5 · p. 25
B.6 Integrating a mixture of smooth and non smooth functions
- The next test uses H(x)(cos(πx) + cos(2πx) + sin(5πx)) and compares BH and EZ.Figure B.8 is identified as a comparison for this subsection.Source: §B.6 · p. 26, Figure B.8 · p. 26