REVIEW 2 major objections 4 minor 36 references
Hierarchical Bayesian Quadrature
T0 review · 2 major / 4 minor · reviewed 2026-07-14 · grok-4.5
Pith's one-line read Hierarchical Bayesian Quadrature adapts stationary GP models to nonstationary integrands by growing a tree of local surrogates and recombining their integrals with hierarchical conditioning.
desk verdict Solid adaptive BQ for nonstationary integrands: real empirical gains, usable algorithm, partial theory and soft hierarchical novelty. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
Tree conditioning: each branch node retains the stationary GP that was fitted before the split; the children’s integral estimates are treated as noisy observations of the corresponding sub-integrals under that parent prior, yielding a joint Gaussian update that reintroduces cross-subdomain correlations and propagates upward to the global integral.
What would settle it
Run the method on a nonstationary integrand whose iso-complexity contours are strongly diagonal or curved (so that axis-aligned rectangles cannot isolate them) and check whether the adaptive tree still reduces integral error relative to a single stationary GP at matched evaluation budgets; if it does not, the partition geometry is insufficient.
Extended reading notes
Core claim
An adaptively grown tree of local stationary Gaussian processes, whose leaf integral estimates are recombined by hierarchical GP conditioning and whose growth is gated by a BIC split test, recovers the accuracy of ordinary Bayesian quadrature on stationary integrands and substantially improves it on nonstationary ones, all without MCMC and while preserving closed-form kernel mean embeddings.
Load-bearing premise
That axis-aligned rectangular splits and stationary product kernels with closed-form mean embeddings are flexible enough to capture the nonstationarity that matters, and that the partial error bounds that assume a finished tree still justify the adaptive construction.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes Hierarchical Bayesian Quadrature (HBQ): an adaptive tree of axis-aligned hyperrectangular subdomains, each carrying a local stationary Matérn (or Wendland) GP with its own lengthscale and smoothness, whose leaf integral estimates are recombined bottom-up by hierarchical GP conditioning that reintroduces parent-kernel cross-subdomain correlations. Tree growth is gated by a BIC approximation to the split Bayes factor (with a 5d minimum points-per-leaf rule), and refinement is driven by a batched integral-variance-reduction score. Theorem 1 / Corollary 2 bound the parent integral error by the worst-case child error under quasi-uniform local designs. Experiments on Genz families (d=2 and d=3), an SIR model-evidence ridge, an on-model piecewise Matérn recovery task with calibration, and a multi-scale reaction–diffusion PDE show substantial gains over standard BQ on nonstationary integrands and parity on stationary ones, with a public Julia implementation.
Significance. If the empirical claims hold, HBQ is a practical Bayesian counterpart to Genz–Malik adaptive cubature that preserves closed-form kernel mean embeddings and calibrated integral uncertainty while automatically allocating budget to local complexity. The combination of BIC-gated splitting, hierarchical recombination of integral functionals, and the partial but clean recombination theory is a genuine contribution to probabilistic numerics for expensive low-to-moderate-dimensional nonstationary integrands (model evidences, multi-scale PDE/simulator responses). Strengths include reproducible code, multiple controlled baselines (MC/QMC/BQ), an on-model structure-recovery experiment with z-score calibration, and an explicit limitations section that correctly flags the partial theory and axis-aligned restriction.
major comments (2)
- [§3.2 / Supp. C.3] §3.2 Eqs. (15)–(17) and Supp. C.3: tree conditioning is presented as the primary recombination mechanism that reintroduces cross-subdomain correlations, yet on the Genz families the independent-sum ablation matches both median absolute error and posterior standard deviation to within 0.05% at every budget. The hierarchical novelty claim therefore rests almost entirely on adaptive partitioning rather than on the recombination hierarchy. Either demonstrate a setting (e.g., the piecewise Matérn or PDE) where tree conditioning materially improves calibration or accuracy, or reframe the contribution so that independent sum is the default and tree conditioning is optional.
- [§3.3 / §7.1] §3.3 Theorem 1 / Corollary 2 and §7.1: the convergence statements are deliberately conditional on a fixed final tree and fixed local kernels; they do not cover the adaptive construction, BIC decisions, or hyperparameter re-fitting that constitute the algorithm. The central claim that HBQ “achieves substantial gains … while matching … on stationary ones” is therefore supported only empirically. A short discussion of what would be required for a full adaptive rate (or an explicit statement that the theory is only a recombination sanity check) would make the theoretical contribution proportionate to what is proved.
minor comments (4)
- [§4.1] §4.1 / Eq. (23): the BIC penalty uses Δk = 3 (shared lengthscale) or d+2 (ARD) and a hard 5d minimum; Supp. C.4 shows robustness on the piecewise task, but a one-sentence justification for the 5d rule (beyond “reliable hyperparameter estimation”) would help readers who wish to retune.
- [Fig. 1 / §6.2] Fig. 1(c) and §6.2: SIR results use independent-sum recombination while the PDE and Genz main figures use tree conditioning; a brief note in the caption or text would avoid confusion about which recombination is active.
- [§5] §5 Related work: Nested kriging and BART-Int are cited appropriately; a short clause distinguishing HBQ’s closed-form integral functionals from Nested kriging’s second-level GP aggregation would sharpen the positioning.
- [Abstract / §5] Typographical: arXiv header and abstract both use “model evidence computations”; ensure consistent hyphenation (“model-evidence”) and that the Genz–Malik citation appears with the same spelling throughout.
Circularity Check
No significant circularity: algorithmic method with external ground-truth benchmarks and standard GP fitting; theory is conditional and openly partial.
full rationale
The paper's central claims are algorithmic (adaptive axis-aligned tree of local stationary Matérn/Wendland GPs, BIC-gated splits, hierarchical recombination of closed-form integral functionals) and empirical (gains on nonstationary Genz Corner Peak / reaction-diffusion / SIR ridge; parity on stationary families). Integral estimates follow standard BQ formulas (Eqs. 5-6, 16-17) applied to local kernels whose hyperparameters are profile-ML fitted per leaf (standard practice, not a derivation of the integral). Ground truths are external (analytic Genz, fine-grid trapezoid, Genz-Malik reference). Theorem 1 / Corollary 2 are explicitly conditional on a fixed final tree and kernels and do not claim to cover adaptive construction or hyperparameter adaptation (§3.3, §7.1). No prediction reduces by construction to a fitted target; no uniqueness theorem or ansatz is imported via self-citation as a load-bearing premise; related-work citations (nested kriging, Genz-Malik, treed GPs) are properly scoped as analogues. The method is self-contained against independent benchmarks.
Assumptions & free parameters
free parameters (6)
- BIC split penalty Δk (and optional multiplier in sensitivity)
- Minimum points per child leaf (5d)
- Child lengthscale log-normal prior variance σ²_prior and parent-anchored mean
- Root PC prior calibration P(ℓ < |Ω|/10) = 0.025
- Candidate Matérn smoothness set {1/2, 3/2, 5/2}
- Refinement batch size / init grid density (e.g. Sobol +20, initdiv)
assumptions (4)
- domain assumption Integrand is well modeled locally by zero-mean stationary Matérn/Wendland GPs with tensor-product structure on axis-aligned hyperrectangles so that kernel mean embeddings remain closed-form.
- ad hoc to paper BIC difference approximates the Bayes factor for split vs no-split well enough at small per-leaf N to control false splits.
- standard math Standard GP posterior and integral formulas (Eqs. 2–6) and Neumann-series expansion for tree conditioning (Theorem 1 proof).
- domain assumption Quasi-uniform leaf designs yield the classical Matérn BQ rates used in Corollary 2.
invented entities (2)
-
Tree-conditioning recombination of leaf integral estimates
-
HBQ adaptive loop (BIC split vs IVR refine on variance-reduction scores)
Cite this review
Pith. "Pith review of Hierarchical Bayesian Quadrature." pith.science (2026). https://pith.science/paper/6QMYJY7K
@misc{pith2026260710793,
author = {Pith},
title = {Pith review of: Hierarchical Bayesian Quadrature},
year = {2026},
howpublished = {\url{https://pith.science/paper/6QMYJY7K}},
note = {Machine review of arXiv:2607.10793}
}
read the original abstract
Numerical integration is a cornerstone of various scientific computing applications, such as engineering simulations and model evidence computations in probabilistic machine learning. Bayesian Quadrature uses Gaussian process surrogates that explicitly encode structural assumptions about the integrand to obtain integral estimates with quantified uncertainty. These surrogates are predominantly based on stationary covariance functions, which results in model misspecification for integrands exhibiting nonstationary behavior. We tackle this issue through an adaptively growing, tree-based partition of the integration domain into local stationary models. Our method recombines the local integral estimates through a hierarchy of GP conditioning that reintroduces cross-subdomain correlations, while model selection criteria control the tree growth to avoid unnecessary partitioning. The resulting algorithm is simple, requires no MCMC, and adapts its evaluation budget to local integrand complexity. On benchmark integration problems and a model evidence computation for an epidemiological model, Hierarchical Bayesian Quadrature achieves substantial gains over standard Bayesian Quadrature on nonstationary integrands while matching its performance on stationary ones.
Figures
Figures from the paper (7 more)
Reference graph
Works this paper leans on
-
[1]
Adachi, S
M. Adachi, S. Hayakawa, M. J rgensen, H. Oberhauser, and M. A. Osborne. Fast B ayesian inference with batch B ayesian quadrature via kernel recombination. In Advances in Neural Information Processing Systems, volume 35, pages 16533--16547, 2022
2022
-
[2]
Berntsen, T
J. Berntsen, T. O. Espelid, and A. Genz. Algorithm 698: DCUHRE : An adaptive multidemensional integration routine for a vector of integrals. ACM Transactions on Mathematical Software, 17 0 (4): 0 452--456, 1991
1991
-
[3]
Briol, C
F.-X. Briol, C. Oates, M. Girolami, and M. A. Osborne. F rank- W olfe B ayesian quadrature: Probabilistic integration with theoretical guarantees. In Advances in Neural Information Processing Systems , volume 28, pages 1162--1170, 2015
2015
-
[4]
Briol, C
F.-X. Briol, C. J. Oates, M. Girolami, M. A. Osborne, and D. Sejdinovic. Probabilistic integration: A role in statistical computation? (with discussion and rejoinder). Statistical Science, 34 0 (1): 0 1--22, 2019
2019
-
[5]
Briol, A
F.-X. Briol, A. Gessner, T. Karvonen, and M. Mahsereci. A dictionary of closed-form kernel mean embeddings. In Proceedings of the 1st International Conference on Probabilistic Numerics, volume 271 of PMLR, pages 84--95, 2025
2025
-
[6]
P. J. Davis and P. Rabinowitz. Methods of Numerical Integration. Computer Science and Applied Mathematics. Academic Press, 2nd edition, 1984
1984
-
[7]
Deisenroth and J
M. Deisenroth and J. W. Ng. Distributed G aussian processes. In Proceedings of the 32nd International Conference on Machine Learning , volume 37 of PMLR, pages 1481--1490, 2015
2015
-
[8]
J. Dick, F. Y. Kuo, and I. H. Sloan. High-dimensional integration: The quasi- M onte C arlo way. Acta Numerica, 22: 0 133--288, 2013
2013
Show all 36 references
-
[9]
Scalable global optimization via local Bayesian optimization
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, pages 5496--5507, 2019
2019
-
[10]
Evans and T
M. Evans and T. Swartz. Approximating Integrals via Monte Carlo and Deterministic Methods . Oxford University Press, 2000
2000
-
[11]
Fuglstad, D
G.-A. Fuglstad, D. Simpson, F. Lindgren, and H. Rue. Constructing priors that penalize the complexity of G aussian random fields. Journal of the American Statistical Association, 114 0 (525): 0 445--452, 2019
2019
-
[12]
A. Genz. Testing multidimensional integration routines. In Proceedings of International Conference on Tools , Methods and Languages for Scientific and Engineering Computation , pages 81--94. Elsevier North-Holland, 1984
1984
-
[13]
A. C. Genz and A. A. Malik. Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N -dimensional rectangular region. Journal of Computational and Applied Mathematics, 6 0 (4): 0 295--302, 1980
1980
-
[14]
M. N. Gibbs. B ayesian G aussian Processes for Regression and Classification . PhD thesis, University of Cambridge, 1997
1997
-
[15]
R. B. Gramacy and H. K. H. Lee. B ayesian treed G aussian process models with an application to computer modeling. Journal of the American Statistical Association, 103 0 (483): 0 1119--1130, 2008
2008
-
[16]
A simple method for comparing complex models: Bayesian model comparison for hierarchical multinomial processing tree models using warp-iii bridge sampling
Quentin Gronau, Eric-Jan Wagenmakers, Daniel Heck, and Dora Matzke. A simple method for comparing complex models: Bayesian model comparison for hierarchical multinomial processing tree models using warp-iii bridge sampling. Psychometrika, 84 0 (1): 0 261--284, 2019. doi:10.100...
2019 doi
-
[17]
Gunter, M
T. Gunter, M. A. Osborne, R. Garnett, P. Hennig, and S. J. Roberts. Sampling for inference in probabilistic models with fast B ayesian quadrature. In Advances in Neural Information Processing Systems , volume 27, pages 2789--2797, 2014
2014
-
[18]
Hennig, M
P. Hennig, M. A. Osborne, and M. Girolami. Probabilistic numerics and uncertainty in computations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471 0 (2179): 0 20150142, 2015
2015
-
[19]
Hennig, M
P. Hennig, M. A. Osborne, and H. P. Kersting. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022
2022
-
[20]
Max Hinne, David Leeftink, Marcel A. J. van Gerven, and Luca Ambrogioni. Bayesian model averaging for nonparametric discontinuity design. PLOS ONE, 17 0 (6): 0 e0270310, 2022. doi:10.1371/journal.pone.0270310
2022 doi
-
[21]
Kanagawa and P
M. Kanagawa and P. Hennig. Convergence guarantees for adaptive B ayesian quadrature methods. In Advances in Neural Information Processing Systems, volume 32, pages 6237--6248, 2019
2019
-
[22]
Kanagawa, B
M. Kanagawa, B. K. Sriperumbudur, and K. Fukumizu. Convergence analysis of deterministic kernel-based quadrature rules in misspecified settings. Foundations of Computational Mathematics, 20 0 (1): 0 155--194, 2020
2020
-
[23]
K. Li, D. Giles, T. Karvonen, S. Guillas, and F.-X. Briol. Multilevel B ayesian quadrature. In Proceedings of the 26th International Conference on Artificial Intelligence and Statistics, volume 206 of PMLR, pages 1845--1868, 2023
2023
-
[24]
Mahsereci and T
M. Mahsereci and T. Karvonen. B ayesian quadrature: G aussian processes for integration. arXiv:2602.16218v1, 2026
2026
-
[25]
Irons, Pierre Latouche, and Adrian E
Martin Metodiev, Marie Perrot-Dock \`e s, Sarah Ouadah, Nicholas J. Irons, Pierre Latouche, and Adrian E. Raftery. Easily Computed Marginal Likelihoods from Posterior Simulation Using the THAMES Estimator . Bayesian Analysis, 20 0 (3): 0 1003 -- 1030, 2025. doi:10.1214/24-BA1422
2025 doi
-
[26]
A. O'Hagan. Bayes-- Hermite quadrature. Journal of Statistical Planning and Inference, 29 0 (3): 0 245--260, 1991
1991
-
[27]
M. A. Osborne, D. Duvenaud, R. Garnett, C. E. Rasmussen, S. J. Roberts, and Z. Ghahramani. Active learning of model evidence using B ayesian quadrature. In Advances in Neural Information Processing Systems, volume 25, pages 46--54, 2012
2012
-
[28]
Paciorek and M
C. Paciorek and M. Schervish. Nonstationary covariance functions for G aussian process regression. In Advances in Neural Information Processing Systems , volume 16, 2003
2003
-
[29]
Park and D
C. Park and D. Apley. Patchwork kriging for large-scale G aussian process regression. Journal of Machine Learning Research, 19 0 (7): 0 1--43, 2018
2018
-
[30]
Pf \"o rtner, I
M. Pf \"o rtner, I. Steinwart, P. Hennig, and J. Wenger. Physics-informed G aussian process regression generalizes linear PDE solvers. arXiv:2212.12474v6, 2024
2024 arXiv
-
[31]
C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006
2006
-
[32]
Remes, M
S. Remes, M. Heinonen, and S. Kaski. Non-stationary spectral kernels. In Advances in Neural Information Processing Systems , volume 30, 2017
2017
-
[33]
Rulli \`e re, N
D. Rulli \`e re, N. Durrande, F. Bachoc, and C. Chevalier. Nested Kriging predictions for datasets with a large number of observations. Statistics and Computing, 28 0 (4): 0 849--867, 2018
2018
-
[34]
G. Schwarz. Estimating the dimension of a model. The Annals of Statistics, 6 0 (2): 0 461--464, 1978
1978
-
[35]
V. Tresp. A B ayesian committee machine. Neural Computation, 12 0 (11): 0 2719--2741, 2000
2000
-
[36]
H. Zhu, X. Liu, R. Kang, Z. Shen, S. Flaxman, and F.-X. Briol. Bayesian probabilistic numerical integration with tree-based models. In Advances in Neural Information Processing Systems , volume 33, pages 5837--5849, 2020
2020
Reviewed July 14, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.