Smooth Monte-Carlo Estimators

• Smooth Monte-Carlo estimators are simulation methods that embed smoothing techniques directly into the estimator design to reduce variance and mitigate outlier effects.
• They employ strategies such as smoothed median-of-means, control functionals, and structured point sets (e.g., DPPs and Frolov lattices) to enhance accuracy and convergence.
• These estimators achieve faster convergence rates and more stable error assessments across applications like randomized approximation schemes, MCMC diagnostics, and density estimation.

Smooth Monte-Carlo estimators are Monte Carlo procedures in which smoothness is not incidental but part of the estimator design. In the literature represented here, smoothness enters in several distinct ways: by smoothing noisy batch means before robust aggregation, by fitting smooth surrogate functions whose integrals are analytically tractable, by replacing i.i.d. sampling with repulsive or lattice-based point sets adapted to smooth integrands, by smoothing discontinuous simulation outputs through conditioning or kernels, by stabilizing noisy autocovariance sequences in MCMC, and by using randomized Taylor truncations to estimate smooth functions of expectations without bias (Huber, 2014, Oates et al., 2016, Gautier et al., 21 Apr 2026, Ullrich, 2016, Abdellah et al., 2018, L'Ecuyer et al., 2019, Dai et al., 2017, Chopin et al., 2024, Salaün et al., 2022).

1. Conceptual scope

A useful way to organize this area is by asking what is being smoothed. In some constructions, the estimator itself is smoothed: the “smooth” median-of-means of randomized approximation schemes multiplies each batch mean by an independent $R \sim \mathrm{Unif}([1-\epsilon,1+\epsilon])$, so that outliers are softened before median aggregation (Huber, 2014). In others, the target integrand is approximated by a smooth control functional or regression surrogate, and Monte Carlo is applied only to the residual (Oates et al., 2016, Salaün et al., 2022). A third class smooths the sampling design rather than the function, using structured point sets such as projection DPPs or randomly shifted and dilated lattices to exploit spectral or Sobolev regularity (Gautier et al., 21 Apr 2026, Ullrich, 2016). A fourth class smooths discontinuous outputs or diagnostics, for example by conditioning to turn indicator-based density estimators into continuous conditional densities, or by imposing PSD and monotonicity constraints on empirical autocovariance sequences in MCMC (L'Ecuyer et al., 2019, L'Ecuyer et al., 2021, Dai et al., 2017).

This literature also suggests that “smoothness” has two mathematically different roles. In one role, smoothness is an assumption on $f$, $\pi$, or $\widehat f$, and sharper convergence follows because Fourier tails, RKHS approximation error, or polynomial coefficients decay rapidly. In the other role, smoothness is a device for variance control even when only weak moments are available, as in relative-error estimation under $\mathrm{SD}(X)\le c\mu$, or in conditional Monte Carlo constructions that replace discontinuous indicators by continuous conditional expectations (Huber, 2014, Oates et al., 2016, Ullrich, 2016, L'Ecuyer et al., 2019).

2. Smoothed relative-error estimators in randomized approximation schemes

A central problem in randomized approximation schemes is to estimate a strictly positive mean $\mu$ from i.i.d. samples $X_1,X_2,\ldots$ under the coefficient-of-variation bound $\sigma \le c\mu$, with guarantee

$\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$

for $\epsilon\in(0,1/3)$ and $f$0. The construction in “Improving Monte Carlo randomized approximation schemes” replaces the classic block-median argument by a smoothed batch statistic. For batch size

$f$1

one forms batch means $f$2, draws $f$3 independently, sets $f$4, and returns the median of $f$5 over

$f$6

batches (Huber, 2014).

The key smoothing lemma states that if $f$7 has mean $f$8 and standard deviation at most $f$9, then with independent $\pi$0,

$\pi$1

This per-batch tail probability $\pi$2 is then propagated through a median concentration bound for $\pi$3 independent batches, giving

$\pi$4

with total sample size

$\pi$5

In leading order,

$\pi$6

improving on the classic median-of-means constant $\pi$7 (Huber, 2014).

The significance of this result is not only the reduced constant. The smoothing step yields the two-sided $\pi$8 batch failure probability without the larger block sizes that Chebyshev-based arguments would require. The paper explicitly situates the method in randomized approximation schemes for the volume of a convex body, the permanent of a nonnegative matrix, the number of linear extensions of a poset, and partition functions of the Ising model. The guarantee depends critically on $\pi$9, independence of samples, and correct specification of the bound $\widehat f$0; if $\widehat f$1, if $\widehat f$2 is extremely small, or if the samples are dependent, the relative-error formulation becomes ill-posed or the concentration argument can fail (Huber, 2014).

3. Control functionals and regression surrogates

A major line of work smooths Monte Carlo by replacing a rough integrand with a learned surrogate whose integral is known. In “Convergence Rates for a Class of Estimators Based on Stein’s Method”, the surrogate has the form

$\widehat f$3

with $\widehat f$4 in an RKHS and $\widehat f$5. Under sample splitting, $\widehat f$6 is used to fit $\widehat f$7 by minimum-norm interpolation in the augmented RKHS $\widehat f$8, and $\widehat f$9 is used to compute

$\mathrm{SD}(X)\le c\mu$0

The estimator is unbiased, and its conditional variance is $\mathrm{SD}(X)\le c\mu$1 (Oates et al., 2016).

The principal rate theorem is a scattered-data approximation statement. If $\mathrm{SD}(X)\le c\mu$2 has smoothness $\mathrm{SD}(X)\le c\mu$3, the kernel has smoothness $\mathrm{SD}(X)\le c\mu$4, the domain satisfies an interior cone condition, and $\mathrm{SD}(X)\le c\mu$5, then in the independent case

$\mathrm{SD}(X)\le c\mu$6

When $\mathrm{SD}(X)\le c\mu$7, this yields

$\mathrm{SD}(X)\le c\mu$8

strictly faster than standard $\mathrm{SD}(X)\le c\mu$9 Monte Carlo whenever $\mu$0. The same MSE rate extends to uniformly ergodic reversible Markov chains via variance-bounding arguments. At the same time, the gain deteriorates with dimension through the factor $\mu$1, which the paper identifies as an inherent curse of dimensionality (Oates et al., 2016).

A related but computationally simpler approach appears in “Regression-based Monte Carlo Integration”. There, Monte Carlo integration is reinterpreted as estimation of a constant function, and then generalized to an analytically integrable surrogate $\mu$2. Writing $\mu$3 and

$\mu$4

one obtains an unbiased control-variate estimator whose variance is $\mu$5. When $\mu$6 is fitted by least squares over a basis that includes the constant function, the resulting estimator is provably better than or equal to conventional Monte Carlo, because the constant model is always available as a fallback. In the rendering experiments reported in the paper, low-order polynomial surrogates acted as a practical drop-in replacement for standard Monte Carlo and combined effectively with MIS (Salaün et al., 2022).

Taken together, these papers suggest a common structural principle: learn a smooth or low-complexity approximation whose integral is available in closed form, then use stochastic sampling only for the residual. The distinction is that Stein control functionals obtain the zero-mean correction from $\mu$7 and RKHS structure, whereas regression-based Monte Carlo relies on basis functions with analytically known expectations (Oates et al., 2016, Salaün et al., 2022).

4. Structured point sets: determinantal processes and randomized lattices

Another family of smooth Monte-Carlo estimators leaves the integrand unchanged and instead regularizes the point set. In “On two ways to use determinantal point processes for Monte Carlo integration”, projection DPPs are used to generate exactly $\mu$8 repulsive nodes from a kernel

$\mu$9

The Bardenet–Hardy estimator is

$X_1,X_2,\ldots$0

which is unbiased. For the multivariate Jacobi ensemble on $X_1,X_2,\ldots$1 and $X_1,X_2,\ldots$2 essentially $X_1,X_2,\ldots$3, the cited CLT implies variance of order $X_1,X_2,\ldots$4, hence RMSE $X_1,X_2,\ldots$5. The Ermakov–Zolotukhin estimator uses the same DPP nodes but solves a linear system in the basis $X_1,X_2,\ldots$6; when $X_1,X_2,\ldots$7 is constant, it yields an unbiased integral estimator whose variance is exactly the residual $X_1,X_2,\ldots$8 energy of $X_1,X_2,\ldots$9 outside the span $\sigma \le c\mu$0. If $\sigma \le c\mu$1, the variance is zero (Gautier et al., 21 Apr 2026).

The same paper emphasizes that the two estimators exploit smoothness differently. BH uses a fixed DPP and benefits when the smoothness of $\sigma \le c\mu$2 yields small Fourier coefficients relative to the arcsine equilibrium measure. EZ is tailored to $\sigma \le c\mu$3 through the spectral residual and can integrate finite expansions exactly, but it may be numerically erratic because the interpolation matrix can be ill-conditioned and the resulting weights need not be nonnegative. Exact chain-rule sampling algorithms are provided for multivariate Jacobi ensembles, with expected rejection count approximately $\sigma \le c\mu$4, and a rejection-free $\sigma \le c\mu$5 tridiagonal sampler is available in one dimension (Gautier et al., 21 Apr 2026).

A different structured-sampling mechanism appears in “A Monte Carlo method for integration of multivariate smooth functions”. There the estimator is a randomly shifted and dilated Frolov lattice rule

$\sigma \le c\mu$6

with $\sigma \le c\mu$7, $\sigma \le c\mu$8, and $\sigma \le c\mu$9. For compactly supported $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$0,

$\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$1

so the RMS error is the high-frequency $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$2-mass of $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$3 outside a growing neighborhood of the origin. This immediately converts Fourier decay into Sobolev rates. For mixed smoothness $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$4 with $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$5, the optimal randomized order on the cube is $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$6, independent of dimension. For isotropic smoothness, the rate becomes $\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$7 (Ullrich, 2016).

The conceptual contrast is instructive. DPP estimators obtain variance reduction through repulsion and spectral projection, while randomized Frolov rules obtain it through controlled aliasing in Fourier space. Both methods, however, convert smoothness into decay of a spectral tail, and both depart from classical Monte Carlo by imposing strong structure on the point set rather than on the estimator algebra (Gautier et al., 21 Apr 2026, Ullrich, 2016).

5. Density estimation and conditioning-based smoothers

Simulation-based density estimation makes the role of smoothing especially explicit. In the scalar KDE setting

$\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$8

standard Monte Carlo gives

$\mathbb P(|\hat\mu-\mu|>\epsilon\mu)\le \delta,$9

so $\epsilon\in(0,1/3)$0 at bandwidth $\epsilon\in(0,1/3)$1. Replacing MC by RQMC or stratification changes only the variance, not the bias. “Density estimation by Randomized Quasi-Monte Carlo” shows that Koksma–Hlawka-type bounds lead to

$\epsilon\in(0,1/3)$2

which is too loose in moderate or high dimension because the Hardy–Krause variation grows like $\epsilon\in(0,1/3)$3. For nested uniform scrambling, however,

$\epsilon\in(0,1/3)$4

so the optimal bandwidth scaling remains $\epsilon\in(0,1/3)$5 while constants can improve substantially. Under monotone $\epsilon\in(0,1/3)$6, stratification yields

$\epsilon\in(0,1/3)$7

with $\epsilon\in(0,1/3)$8 (Abdellah et al., 2018).

The review “Density Estimation by Monte Carlo and Quasi-Monte Carlo” broadens this picture by comparing KDE with conditional density estimators, smoothed perturbation analysis, LR, and GLR-U constructions. The central identity is $\epsilon\in(0,1/3)$9, $f$00. The problem with the empirical cdf is that it is discontinuous in the underlying uniforms and its derivative is a sum of deltas. Conditioning replaces the indicator by a smooth conditional cdf or conditional density. In “Monte Carlo and Quasi-Monte Carlo Density Estimation via Conditioning”, if $f$01 is continuous in $f$02, differentiable except on a countable set, and $f$03 is uniformly bounded by an $f$04 random variable $f$05, then

$f$06

is unbiased, with $f$07. Under bounded Hardy–Krause variation of the derivative integrand, RQMC improves the pointwise MSE to $f$08. The paper reports especially large gains for the cantilever beam, stochastic activity network, and Asian option examples, with empirical log–log slopes near $f$09 for several CDE+RQMC configurations (L'Ecuyer et al., 2021, L'Ecuyer et al., 2019).

At a more general level, “A Monte Carlo Method to Approximate Conditional Expectations based on a Theorem of Besicovitch” justifies shrinking-neighborhood Monte Carlo approximations of conditional expectations by a differentiation-of-measures theorem. For $f$10 and integrable $f$11, the paper shows that for $f$12-almost every $f$13,

$f$14

providing a nonparametric conditional-expectation mechanism even when densities are unavailable (Nogales et al., 2013).

Across these density-estimation papers, smoothing serves a precise purpose: it transforms discontinuous, high-variation simulation outputs into continuous functions of the driving uniforms, thereby making RQMC, conditioning, and derivative-based arguments effective (Abdellah et al., 2018, L'Ecuyer et al., 2021, L'Ecuyer et al., 2019, Nogales et al., 2013).

6. MCMC covariance smoothing and stable Monte Carlo error assessment

Not all smooth Monte-Carlo estimators target means or integrals directly. In “Multivariate initial sequence estimators in Markov chain Monte Carlo”, smoothness refers to stabilization of the estimated long-run covariance matrix in the multivariate MCMC CLT

$f$15

Writing $f$16 for adjacent-sum autocovariances, Proposition 1 shows that in the population sequence each $f$17 is positive definite, $f$18 is positive definite, and $f$19. The empirical counterparts $f$20 are noisy, and naive truncation can lead to unstable or non-PSD covariance estimates (Dai et al., 2017).

The multivariate initial-sequence estimator mIS addresses this by enforcing shape constraints. One forms

$f$21

lets $f$22 be the first index for which $f$23 is positive definite, and then takes $f$24 as the last index before the determinant sequence $f$25 ceases to increase. The estimator is $f$26. The adjusted version mISadj replaces each $f$27 beyond $f$28 by its PSD projection $f$29, ensuring that increments cannot shrink the confidence ellipsoid in negative-eigenvalue directions. Theorems 2 and 3 show that both mIS and mISadj conservatively estimate the generalized variance, in the sense that

$f$30

almost surely (Dai et al., 2017).

The practical role of smoothing here is shape-constrained regularization of a noisy time-domain spectral estimate. The paper contrasts mIS and mISadj with Kosorok’s multivariate estimator mK, which truncates at the first PSD violation and can underestimate generalized variance, and with componentwise univariate initial-sequence estimators, which ignore cross-correlation and yield cube-shaped rather than ellipsoidal uncertainty regions. The result is not differentiability but stability: positive-definite covariance estimates, more reliable Monte Carlo standard errors, and confidence ellipsoids with empirically improved coverage (Dai et al., 2017).

7. Unbiased nonlinear transforms of expectations and recurrent limitations

A distinct use of smoothness appears in “Towards a turnkey approach to unbiased Monte Carlo estimation of smooth functions of expectations”. Here the objective is not $f$31 itself but $f$32, where $f$33 and $f$34 is analytic, with special emphasis on $f$35 and $f$36. Expanding around $f$37,

$f$38

one can unbiasedly estimate $f$39 by products of independent terms $f$40. The paper studies both the simple estimator $f$41 and a variance-reduced cycling estimator $f$42 that averages cyclic products across all $f$43 samples. With an independent truncation $f$44 and weights $f$45, the sum estimator

$f$46

is unbiased for $f$47 (Chopin et al., 2024).

The finite-variance theory is driven by

$f$48

For geometric truncation $f$49, the random-truncation variance behaves like $f$50, but the sampling-noise term behaves very differently for the two estimators: the simple estimator does not have vanishing noise variance as $f$51, while the cycling estimator satisfies an $f$52 bound, i.e. near Monte Carlo rate up to a logarithmic factor. The paper proposes automatic tuning via a pilot sample, with

$f$53

and a one-sided bootstrap safeguard for the minimum admissible $f$54 (Chopin et al., 2024).

This line of work is especially relevant to latent-variable likelihood estimation and unnormalized models, because $f$55 and $f$56 arise naturally in those settings. Its limitations are equally explicit: for $f$57 and $f$58, one needs $f$59, $f$60, $f$61, and geometric rather than sub-geometric truncation tails; heavy-tailed or high-variance $f$62 can make $f$63 close to $f$64, forcing very small $f$65 and therefore large expected truncation levels (Chopin et al., 2024).

Viewed together, these papers reveal recurrent trade-offs. Gains from smoothness often deteriorate with dimension through scattered-data fill distance, polynomial or Fourier residuals, or DPP sampling cost (Oates et al., 2016, Gautier et al., 21 Apr 2026, Ullrich, 2016). Some methods require access to $f$66, compact domains, or boundary mollifiers (Oates et al., 2016). Others depend on positivity and correctly specified variance surrogates, as in relative-error randomized approximation schemes and unbiased reciprocal or logarithmic transforms (Huber, 2014, Chopin et al., 2024). In RQMC density estimation, classical Koksma–Hlawka bounds become unusable because variation blows up as $f$67, even though empirical performance can still be strong (Abdellah et al., 2018). In MCMC covariance smoothing, reversibility, a multivariate CLT, and reliable autocovariance estimation remain essential (Dai et al., 2017).

The unifying implication is not that there is a single canonical “smooth Monte-Carlo estimator”, but that multiple Monte Carlo subfields repeatedly exploit the same principle: replace a discontinuous, noisy, or spectrally diffuse object by one with more regular structure, then use stochastic sampling on what remains. The specific regularity may be uniform scaling of batch means, RKHS smoothness, polynomial or Fourier sparsity, repulsive node geometry, conditional differentiability, PSD shape constraints, or analyticity of $f$68. The estimator class changes, but the operational objective is stable: lower variance, stronger nonasymptotic concentration, or more reliable Monte Carlo error quantification from the same underlying simulation budget (Huber, 2014, Oates et al., 2016, Gautier et al., 21 Apr 2026, L'Ecuyer et al., 2019, Dai et al., 2017, Chopin et al., 2024).