Adaptive stratified Monte Carlo using decision trees

Hejin Wang; Mathieu Gerber; Nicolas Chopin

REVIEW 2 major objections 4 minor 17 references

Adaptive stratified Monte Carlo using decision trees

T0 review · 2 major / 4 minor · reviewed 2026-08-10 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Adaptive stratified Monte Carlo via decision trees provably converges faster than plain Monte Carlo — at rate $O(N^{-1/2-r'})$ for monotone integrands — while keeping $O(N \log N)$ cost.

desk verdict Good adaptive stratification idea, but the proof of the main rate claim contains a reversed inequality and the headline result is unproven as written. read the letter →

arxiv 2501.04842 v1 pith:725BRMVD submitted 2025-01-08 stat.CO

Nicolas Chopin , Hejin Wang , Mathieu Gerber This is my paper · ORCID

classification stat.CO MSC 65C0562G0865D30

keywords adaptivestratifiedMonteCarlodecisiontreesintegrationvariancereductionconvergencerateunbiasedestimationsparsity-adaptivehigh-dimensional

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper proposes AdaStrat, an unbiased Monte Carlo estimator for integrals over $[0,1]^s$ whose strata are learned from the integrand. The idea is to take a first $N$-point sample, build a decision tree that splits rectangles in half along the direction that most reduces empirical variance, and then use the tree's leaves as strata for a second independent sample. The paper claims that for integrands that are strictly increasing in every coordinate and obey a two-sided Lipschitz condition, this adaptive estimator converges at rate $O(N^{-1/2-r'})$ for some $r' > 0$, improving on the $O(N^{-1/2})$ rate of plain Monte Carlo. Because the tree ignores coordinates that do not affect the integrand, the method is sparsity-adaptive, and it stays practical in dimensions where classical stratified estimators such as Haber's require $k^s$ cells. A reader should care because it offers a cheap, unbiased way to beat plain Monte Carlo when the integrand is modestly regular and the dimension is moderate to high.

What carries the argument

The load-bearing mechanism is the variance-reduction factor of a single binary split. For a rectangle $R$ split in half along axis $j$, the stratified variance is multiplied by $(\Delta(R_+^{(j)})+\Delta(R_-^{(j)}))/(2\Delta(R))$, where $\Delta$ is the conditional variance of $f$. Under monotonicity and the two-sided Lipschitz bounds, the best coordinate $j^\star$ makes this factor at most $\frac{\frac43\sum_i \beta_i^2/\alpha_i^2}{1+\frac43\sum_i \beta_i^2/\alpha_i^2}$, uniformly over rectangles, so the $k$ splitting levels compound into a geometric variance decay. To pass from the oracle tree to the feasible one, the paper introduces an $\varepsilon$-oracle tree, which may choose any near-optimal direction when two directions are nearly tied, and uses a concentration inequality for empirical variances to show the estimated tree matches the $\varepsilon$-oracle up to depth $L_k$ with probability $1-O(N^{-2})$.

What would settle it

Take $f(x)=\exp(\sum_{i=1}^s \lambda_i x_i)$ with $\lambda_i=i^{-2}$, satisfying Theorem 2, compute $r(s)$ from the implied $\alpha_i,\beta_i$, and estimate the RMSE exponent from $\log\mathrm{RMSE}$ versus $\log N$ over $N=2^{10},2^{12},\dots,2^{20}$ with many independent runs. If the exponent is not strictly better than $-1/2$ (or fails to reach $-(1/2+r')$ for an $r'<r(s)$), the paper's central rate claim is false; alternatively, for $f=\lambda x_1+\sin(2\pi x_2)$, the paper's counterexample, the predicted exponent should be no better than $-1/2$.

Watch

Extended reading notes

Core claim

The paper's central claim is Theorem 3: under the assumptions of Theorem 2 — $f$ strictly increasing in each coordinate and satisfying coordinate-wise lower and upper Lipschitz bounds — the AdaStrat estimator has RMSE $O(N^{-1/2-r'})$ for any $0<r'<r(s)$. The oracle rate is $r(s)=\frac{1}{2\log 2}\log\left(\frac{\frac34+\sum_i \beta_i^2/\alpha_i^2}{\sum_i \beta_i^2/\alpha_i^2}\right)$. The proof first analyses an oracle tree that splits using true conditional variances and shows that each split reduces variance by a factor bounded away from one; it then shows that the estimated tree agrees with an $\varepsilon$-oracle tree up to depth $L_k$ with probability $1-O(N^{-2})$, so the same rate survives estimation of the tree. The paper also proves a sparsity-adaptive rate $O(N^{-1/2-r(s_0)})$ for linear functions supported on $s_0$ coordinates, and gives Example 1 showing that without monotonicity no such rate can hold in general.

Load-bearing premise

The rate guarantee rests on the integrand being strictly increasing in every coordinate and satisfying a lower Lipschitz bound in every coordinate; if any coordinate is non-monotone, the oracle tree may never split it and the proof of a faster-than-$N^{-1/2}$ rate collapses, as the paper's sine example shows.

Editorial extensions

If this is right

For any integrand in the paper's monotone Lipschitz class, AdaStrat's RMSE beats the $N^{-1/2}$ Monte Carlo benchmark by a positive power of $N$.
Because only influential coordinates are split, the rate depends on the effective dimension $s_0$ for linear functions; the method automatically ignores irrelevant dimensions where fixed-grid stratification wastes evaluations.
The estimator is unbiased and costs $O(N \log N)$ to build, so it can be plugged into pseudo-marginal samplers or stochastic approximation routines that require unbiased estimates of an integral.
The extension to arbitrary $N$ by volume-proportional splits lets the same adaptive strategy be used when $N$ is not a power of two, and taking two points per stratum gives a simple variance estimate.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The $\varepsilon$-oracle proof suggests the practical bottleneck is not tree depth but near-ties between split directions; averaging over several near-equal trees (a forest-style extension) would be a natural way to stabilise the rate constant in high dimension.
The same decision-tree partition could be reused as a control variate, since the tree leaves define a piecewise-constant function whose integral is exactly computable; this would combine stratification and centering without extra evaluations of $f$.
If monotonicity is the real barrier, applying a monotone transformation to each coordinate of the preliminary sample before growing the tree might extend the rates to non-monotone but coordinatewise monotone transforms of $f$; this is not tested in the paper.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Here's my take. The idea is genuinely worth knowing: train a CART-style regression tree on a pilot sample, use the leaves as strata, and get a stratified estimator with O(N log N) total cost. That is a real step beyond Haber-style cubic grids (which need N = k^s) and beyond kNN or Gaussian-process control variates with quadratic or cubic cost. The writing is clear, the sparsity-adaptivity point is nice, Theorem 1 is clean, and Example 1 is an honest demonstration that the monotonicity assumption is not a technical artifact. The authors also cite the relevant prior work and do not oversell the empirical results.

The load-bearing issue is in the proof of Theorem 2. The claimed flaw is real. The proof bounds the variance-reduction factor for any split direction j by h(B_j) = B_j/(1+B_j). It then defines k = argmax_j alpha_j^2 mu_j^2 and observes that B_j >= B_k. Since h is strictly increasing, B_j >= B_k implies h(B_j) >= h(B_k), so that lower bound cannot be used to conclude the factor is <= h(B_k). The uniform contraction factor in Eq. (9) is therefore not established, and Theorem 3, which inherits Eq. (9), is unsupported as written. This is not a missing technical detail; the direction of the inequality is wrong. A different argument might salvage the rate, for example by controlling B_{j*} directly, but that argument is not in the manuscript.

Everything else is minor by comparison: no code, no error bars on the RMSE plots, and the monotone/Lipschitz class is narrow, but the authors are explicit about that scope.

Who this is for: Monte Carlo researchers and practitioners who want a cheap adaptive variance-reduction method. The numerics suggest the method can help even in dimension 50, independent of the rate proof. I would not cite it for the rate, but I would mention it as related work on adaptive stratification.

Recommendation: this deserves a serious referee, because the idea and experiments are worth engaging with. The referee should be asked to check Theorem 2's proof carefully. In current form the headline convergence-rate claim is unproven; the paper would need a repaired proof before it is fully established.

Referee Report

2 major / 4 minor

Summary. The paper proposes AdaStrat, an adaptive stratified Monte Carlo estimator for integrals on [0,1]^s. A preliminary sample of size N is used to grow a regression tree with the CART split criterion; the leaves of the tree define strata of equal volume, and a second independent sample is drawn uniformly from each stratum. The paper claims that for certain integrands (linear functions, or functions that are strictly increasing with two-sided Lipschitz bounds), the oracle-tree version of the estimator converges at rate O(N^{-1/2-r(s)}) with explicit r(s)>0, and that the practical estimator retains this rate up to an arbitrarily small loss (Theorem 3). Numerical experiments compare the method with standard Monte Carlo and Haber's estimator in dimensions up to 50.

Significance. If the rate claims hold after repairing the proof, the paper makes a solid contribution: it provides an adaptive stratification method with near-linear computational cost and explicit convergence-rate gains, in a regime where classical stratified estimators such as Haber's are impractical. The oracle analysis in Theorems 1 and 2 is a clean way to obtain explicit rates, and the paper is unusually transparent about the scope of its assumptions, including a counterexample (Example 1) showing that the monotonicity condition is not a technical artifact. A notable strength is that the method has no data-fitted constants: the only tuning parameter is the tree depth, set to k, so the rate analysis does not rely on learned parameters. The central claim currently rests on a proof of Theorem 2 that contains a reversed inequality; the issue appears locally fixable, but until it is corrected the main theorem is not established.

major comments (2)

[Section 3.2, proof of Theorem 2] After Eq. (8), the variance-reduction factor for a split along direction j is bounded above by h(x_j)=x_j/(1+x_j), where x_j=(4/3)∑_{i=1}^s (β_i^2/α_i^2)(α_i^2 μ_i^2)/(α_j^2 μ_j^2). Since k=argmax_j α_j^2 μ_j^2, one has x_j ≥ x_k for every j, and because h is strictly increasing on [0,∞), h(x_j) ≥ h(x_k). The manuscript instead concludes the opposite inequality and uses it to obtain the contraction factor in Eq. (9); the direction is reversed and Theorem 2 is not proved as written. The argument can be repaired: the oracle split j*=argmin_j Δ(R,j) satisfies C2(j*)≤C2(k), and because Δ(R)=C1(j)+C2(j) is the same for all j, the factor satisfies C2(j*)/Δ(R) ≤ C2(k)/(C1(k)+C2(k)) ≤ h(x_k). This repair should be stated explicitly.
[Section 3.3, proof of Theorem 3] The proof bounds the probability that the estimated tree differs from an ε-oracle tree before depth L_k by O(N^{-2}) and separately bounds the variance of an ε-oracle tree, but it never assembles the unconditional RMSE. Since the estimator is unbiased conditional on the tree, the proof should explicitly write var(bI_AdaStrat)=E[var(bI_AdaStrat|tree)] ≤ P(good) q_ε^{L_k} var(bI_MC) + P(bad) C/N, with P(bad)=O(N^{-2}) and C a constant (for instance using var(bI_AdaStrat|tree)≤var(bI_MC) for any tree), and then conclude the stated rate. The current sentence 'we conclude the proof by remarking...' skips this decomposition. The final display also says 'var[bIMC]' where the variance of the AdaStrat estimator is meant; this is a typo, but the missing composition is a load-bearing gap in the proof of the main claim.

minor comments (4)

[Section 1.2] After the display of the control variate bound, the text says '∥f−bf_N∥∞ = O(k^{-1}) = O(N^{1/s})'; the exponent on N should be −1/s.
[Section 2.1] In the definition of R[j]_−, the displayed formula for (a'_i,b'_i) is self-referential and contains a typo; it should define the interval split at the midpoint of the j-th coordinate.
[Section 3.3, Lemma 2] The bound is written with 'eκ' in two places; this should be 'e^κ'.
[Section 3.3, proof of Theorem 3] The 'Remark that' condition for the choice of L_k is written in a confusing way: the role of δ is not defined (presumably δ=e^{-κ}), and the chain of implications should be rewritten for readability.

Circularity Check

0 steps flagged · score 1.0 of 10

No significant circularity: the rate proof derives the variance contraction from the stated Lipschitz assumptions and controls the estimated tree via concentration inequalities; self-citations are not load-bearing.

full rationale

The central claim, Theorem 3, is derived rather than assumed. The oracle-tree analysis in Theorem 2 starts from the exact variance decomposition (8) and the exact variance-reduction factor (7); the only inputs are the assumed lower and upper Lipschitz bounds on f, and the contraction factor (9) is a theorem consequence, not a fitted quantity. The estimated-tree step in Theorem 3 uses the Maurer-Pontil inequality and Lemma 2 to bound the probability of deviating from the epsilon-oracle tree; the constants there depend on the assumed alpha_i, beta_i bounds and the fixed depth k, not on estimated parameters. The paper also states in Section 3.1 that it established its results 'directly, from first principles,' and the proof does not invoke a prior stratification-rate theorem as a premise. Self-citations (Chopin and Gerber 2024; Oates et al. 2017; Chopin and Ridgway 2017) appear only as background, benchmarks, or experimental setup, and none is load-bearing for the convergence-rate theorems. The paper explicitly delimits its scope: strict monotonicity is required for the rate results, and the method does not reach optimal rates for C^r functions; these are stated limitations, not hidden inputs. The referee-identified inequality reversal in the proof of Theorem 2 is a potential correctness defect in the proof, but it is not circularity: it concerns the direction of a bound involving the increasing function x/(1+x), not an equivalence between the theorem's conclusion and its assumptions.

Assumptions & free parameters 1 free parameters · 4 assumptions · 0 invented entities

The method itself has no fitted free parameters: the only tuning choice is the tree depth, set to k (full depth). The theoretical rates rest on explicit function-class assumptions and standard concentration inequalities. No new entities are postulated.

free parameters (1)

maximum tree depth d = d = k (N = 2^k)
The convergence-rate results (Theorems 1-3) are proven for d = k, i.e., growing the tree to full depth so the number of strata equals N. The paper acknowledges that this would overfit in standard tree regression, but in stratification splitting always reduces variance. This is a design choice rather than a data-fitted constant, but the theory depends on it.

assumptions (4)

domain assumption The integrand f satisfies f(x) = lambda^T x (Theorem 1) or is strictly increasing with coordinate-wise lower and upper Lipschitz bounds alpha_i and beta_i (Theorem 2).
The convergence-rate theorems are stated for these function classes. The practical method does not require the user to know these bounds, but the theoretical rate claims only apply to these classes.
standard math Maurer and Pontil (2009) Theorem 10, an empirical-variance concentration inequality, holds as stated.
Used in Lemma 1 and Lemma 2 to bound the probability that the estimated tree makes an erroneous split decision.
domain assumption The decision tree is grown with mid-point splits along axes, using the CART criterion (4), to full depth k.
The method and its theory are defined for this specific splitting protocol. Section 2.3 sketches a generalization to arbitrary N, but the proofs focus on the power-of-two case.
standard math Splitting a stratum into two equal sub-strata always reduces the variance of the stratified estimator.
This follows from the law of total variance, and the paper relies on it in Section 3.1 and in the proof of Theorem 3.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Adaptive stratified Monte Carlo using decision trees." pith.science (2026). https://pith.science/paper/725BRMVD

@misc{pith2026250104842,
  author       = {Pith},
  title        = {Pith review of: Adaptive stratified Monte Carlo using decision trees},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/725BRMVD}},
  note         = {Machine review of arXiv:2501.04842}
}

abstract

It has been known for a long time that stratification is one possible strategy to obtain higher convergence rates for the Monte Carlo estimation of integrals over the hyper-cube $[0, 1]^s$ of dimension $s$. However, stratified estimators such as Haber's are not practical as $s$ grows, as they require $\mathcal{O}(k^s)$ evaluations for some $k\geq 2$. We propose an adaptive stratification strategy, where the strata are derived from a a decision tree applied to a preliminary sample. We show that this strategy leads to higher convergence rates, that is, the corresponding estimators converge at rate $\mathcal{O}(N^{-1/2-r})$ for some $r>0$ for certain classes of functions. Empirically, we show through numerical experiments that the method may improve on standard Monte Carlo even when $s$ is large.

Figures

Figures reproduced from arXiv: 2501.04842 by the authors.

**Figure 1.** Cubic stratification in dimension s = 5, with k = 5, N = 25. The black dots represent the centers cn. We mention one last possible control variate strategy: assume N = k s , for some integer k ≥ 2, and split the hyper-cube [0, 1]s into N hyper-cubes of edge length 1/k (and hence volume 1/N). Let Cn, n = 1, . . . , N, denote these N ‘sub-cubes’, let cn be the centre of Cn, and define fbN as fbN (x) = X N n=1 f(cn)1Cn… view at source ↗

**Figure 2.** Partition of [0, 1]2 into rectangles (left) that correspond to a decision tree (right). assuming |A| ≥ 2. In case |A| ≤ 1, we abuse notations and set EmpVar(A) = +∞. This non-standard convention will be convenient in the specific context of this paper. We denote by RMSE Ib the RMSE of an estimator Ib. Since we deal only with unbiased estimator, this quantity is always the square root of var(Ib), the variance of I… view at source ↗

**Figure 3.** Toy example, s = 5. RMSE vs N (estimated from 256 independent realisations) for the following unbiased estimators: standard Monte Carlo (MC, red), Haber of order one (Haber, black), and our adaptive decision-tree strategy (DT, yellow). Dashed lines represent the theoretical convergence rates of these three estimators [PITH_FULL_IMAGE:figures/full_fig_p011_3.png] view at source ↗

Figures from the paper (3 more)

**Figure 4.** Figure 4: Toy example, s = 15 (left) and s = 50 (right). Same plot as [PITH_FULL_IMAGE:figures/full_fig_p011_4.png]

**Figure 5.** Figure 5: Toy example, s = 5, N = 210. Box-plots of variance estimates for the Monte Carlo estimator and our adaptive stratified estimator. 4.2 Bayesian model choice We consider a Bayesian statistical model, with parameter β ∈ R d , prior distribution p(β), and likelihood L(y|β)…

**Figure 6.** Figure 6: Bayesian model choice example, s = 5, 10, 15, 20: RMSE (over 256 independent realisations) versus N for the following unbiased estimators: standard Monte Carlo (MC, red), and our adaptive stratified Monte Carlo estimators (DT, yellow). Specifically, we assume a Gaussia…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

17 extracted references · 16 canonical work pages

[1]

and Roberts, G

Andrieu, C. and Roberts, G. O. (2009). The pseudo-marginal approach for efficient M onte C arlo computations. Ann. Statist. , 37(2):697--725

work page 2009
[2]

Bahvalov, N. S. (1959). Approximate computation of multiple integrals. Vestnik Moskov. Univ. Ser. Mat. Meh. Astr. Fiz. Him. , 1959(4):3--18

work page 1959
[3]

and Scornet, E

Biau, G. and Scornet, E. (2016). A random forest guided tour. TEST , 25(2):197--227

work page 2016
[4]

and Gerber, M

Chopin, N. and Gerber, M. (2024). Higher-order M onte C arlo through cubic stratification. SIAM J. Numer. Anal. , 62(1):229--247

work page 2024
[5]

and Ridgway, J

Chopin, N. and Ridgway, J. (2017). Leave P ima I ndians alone: binary regression as a benchmark for B ayesian computation. Statist. Sci. , 32(1):64--87

work page 2017
[6]

Haber, S. (1966). A modified M onte- C arlo quadrature. Mathematics of Computation , 20(95):361--368

work page 1966
[7]

Haber, S. (1967). A modified M onte- C arlo quadrature. ii. Mathematics of Computation , 21(99):388--397

work page 1967
[8]

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The elements of statistical learning , volume 2 of Springer Series in Statistics . Springer, New York, second edition. Data mining, inference, and prediction

work page 2009

Show all 17 references

[9]

Leluc, R., Portier, F., Segers, J., and Zhuman, A. (2023). Speeding up Monte Carlo Integration : Control Neighbors for Optimal Convergence . arxiv 2305.06151

2023 arXiv
[10]

Lemieux, C. (2009). Monte Carlo and Quasi-Monte Carlo Sampling (Springer Series in Statistics) . Springer

2009
[11]

and Pontil, M

Maurer, A. and Pontil, M. (2009). Empirical B ernstein bounds and sample-variance penalization. In COLT 2009 - The 22nd Conference on Learning Theory, Montreal, Quebec, Canada, June 18-21, 2009

2009
[12]

Morgan, J. N. and Sonquist, J. A. (1963). Problems in the analysis of survey data, and a proposal. J. Am. Stat. Assoc. , 58:415--434

1963
[13]

Novak, E. (2016). Some results on the complexity of numerical integration. Monte Carlo and Quasi-Monte Carlo Methods , pages 161--183

2016
[14]

J., Girolami, M., and Chopin, N

Oates, C. J., Girolami, M., and Chopin, N. (2017). Control functionals for M onte C arlo integration. J. R. Stat. Soc. Ser. B. Stat. Methodol. , 79(3):695--718

2017
[15]

Owen, A. B. (2018). Monte Carlo theory, methods and examples (in progress, available on the author's web-page)

2018
[16]

and Monro, S

Robbins, H. and Monro, S. (1951). A stochastic approximation method. Ann. Math. Statistics , 22:400--407

1951
[17]

Scornet, E., Biau, G., and Vert, J.-P. (2015). Consistency of random forests. Ann. Statist. , 43(4):1716--1741

2015

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.

[1] [1]

and Roberts, G

Andrieu, C. and Roberts, G. O. (2009). The pseudo-marginal approach for efficient M onte C arlo computations. Ann. Statist. , 37(2):697--725

work page 2009

[2] [2]

Bahvalov, N. S. (1959). Approximate computation of multiple integrals. Vestnik Moskov. Univ. Ser. Mat. Meh. Astr. Fiz. Him. , 1959(4):3--18

work page 1959

[3] [3]

and Scornet, E

Biau, G. and Scornet, E. (2016). A random forest guided tour. TEST , 25(2):197--227

work page 2016

[4] [4]

and Gerber, M

Chopin, N. and Gerber, M. (2024). Higher-order M onte C arlo through cubic stratification. SIAM J. Numer. Anal. , 62(1):229--247

work page 2024

[5] [5]

and Ridgway, J

Chopin, N. and Ridgway, J. (2017). Leave P ima I ndians alone: binary regression as a benchmark for B ayesian computation. Statist. Sci. , 32(1):64--87

work page 2017

[6] [6]

Haber, S. (1966). A modified M onte- C arlo quadrature. Mathematics of Computation , 20(95):361--368

work page 1966

[7] [7]

Haber, S. (1967). A modified M onte- C arlo quadrature. ii. Mathematics of Computation , 21(99):388--397

work page 1967

[8] [8]

Hastie, T., Tibshirani, R., and Friedman, J. (2009). The elements of statistical learning , volume 2 of Springer Series in Statistics . Springer, New York, second edition. Data mining, inference, and prediction

work page 2009