PaperlayerView exact revision on arXiv

Source-linked AI summary

A Residual Tree Gaussian Process Modeling Framework for High-Dimensional Data

Pulong Ma, Li Ma

arXiv:2610.02893v1stat.MEmath.STstat.COstat.MLstat.TH

TL;DR

Large spatial datasets with heterogeneous structure are challenging to model with Gaussian processes on high-dimensional domains. ResTGP uses residual processes along an adaptively learned tree to model multiscale dependence, with inference linear in sample size given the tree. The paper also establishes posterior consistency for estimating continuous functions.

  • Problem

    Existing GP methods have limitations in addressing large spatial datasets and heterogeneous structures on multidimensional domains beyond geostatistical applications.

    Source: §1 · p. 2

  • Method

    ResTGP recursively constructs predictive and residual processes along a random partition tree, using node residuals to represent finer-level dependencies.

    Source: §1 · p. 4

  • Results

    Given the tree, likelihood computation has cost linear in sample size, and the paper establishes posterior consistency for estimating continuous functions.

    Source: §1 · p. 5

  • Takeaways & Limitations

    The model supports Bayesian learning of tree partitions while representing large- and small-scale dependencies at different resolutions.

    Source: §1 · p. 4

    The framework uses axis-aligned partitions and a single tree structure; the authors note that adaptive non-axis-aligned partitions may be beneficial.

    Source: §7 · p. 23

Abstract

from arXiv · show

With the advance of measurement technologies and increasing computing power, large spatial data with heterogeneous structures are often collected over high-dimensional domains. Existing Gaussian process (GP) models and computational strategies are often inadequate for analyzing such datasets in multi-dimensional domains. To address these challenges, we develop a Bayesian residual tree GP methodology called ResTGP for large spatial data with potentially heterogeneous structures in multi-dimensional domains. The key idea is to decompose a Gaussian process at a cascade of resolutions along a dyadic tree through iteratively computing predictive and residual processes so that the residual process on each tree node, both interior and leaf, becomes sufficient for the finer-level dependency within that node. This allows characterization of the underlying covariance structure in a flexible, multi-scale manner while achieving divide-and-conquer on the data domain, which leads to computational efficiency. To allow efficient tree inference, we introduce a computational strategy for Bayesian inference based on recursive message passing, which scales linearly with the sample size given the tree. This paper also proves posterior consistency of the model for estimating continuous functions in a nonparametric regression framework. Extensive numerical examples and the storm surge application confirm the advantages of the proposed method.

1 Introduction

The paper addresses the difficulty of modeling large, heterogeneous spatial data on high-dimensional domains with Gaussian processes. ResTGP combines adaptive tree learning with multiresolution dependence modeling through residual processes, enabling efficient inference.

1 Introduction

  • Existing GP approaches often do not jointly address large datasets, heterogeneous structure, and domains with dimension beyond geostatistical applications.Some spatial GP methods are not designed for higher-dimensional domains, while other approaches to nonstationarity cannot handle large datasets in multidimensional domains.Source: §1 · p. 2
  • Treed GPs model dependence only within leaves, whereas multiresolution approaches relax that restriction but generally lack adaptive tree learning or divide-and-conquer computation.When the tree is correctly specified, multiresolution methods can capture long- and short-range dependence more effectively than treed GPs.Source: §1 · pp. 3–4
  • ResTGP combines stochastic tree learning with multiresolution GP modeling to retain adaptive partitions without imposing global blockwise independence.The model instead uses resolution-specific blockwise independence to factorize the likelihood over tree nodes.Source: §1 · pp. 4–5
  • ResTGP recursively constructs predictive and residual processes from root to leaves, with each node’s residual serving as a sufficient statistic for finer dependencies within that node.The residual is passed down the tree as the data for the next level.Source: §1 · p. 4
  • Bayesian tree inference favors splitting when a node’s residuals are not well captured by the baseline covariance model, up to a maximum tree depth.This lets the model represent large-scale dependence at lower resolutions and small-scale dependence at deeper resolutions.Source: §1 · p. 4
  • Given the tree, the likelihood factorization supports Bayesian computation linear in sample size, and the paper establishes posterior consistency for estimating continuous functions.The cost claim is under a per-resolution blockwise independence assumption.Source: §1 · p. 5

2 A Residual Tree GP Modeling Framework

ResTGP recursively decomposes a Gaussian process along a stochastic tree, using knot-based predictive processes and residuals to model dependence across scales. Its baseline models, adaptive partition prior, and knot selection support heterogeneous spatial behavior while retaining computationally efficient inference.

2.1 Overview of the model

  • Recursive predictive and residual processes make each residual sufficient for the next iteration, with computation focused on observations not yet selected as knots.At each node, the residual is formed by subtracting the knot-based predictor from the process.Source: §2.1 · pp. 5, 7
  • The decomposition yields mutually independent predictive and residual components, and factors the sampling density across knot sets and the final residual process.This exact decomposition alone does not provide computational benefits; the paper couples it with tree-based modeling strategies.Source: §2.1 · p. 7

2.2 The baseline model M0(A) on leaf nodes

  • Stochastic trees allow across-tree nonstationarity, while region-specific baseline models allow within-tree heterogeneity as node regions shrink.Leaf options include a residual-covariance predictive-process baseline and a heteroskedastic white-noise baseline; flexible parametric covariance models are left for future investigation because large-sample marginal likelihoods remain challenging.Source: §2.2 · pp. 7–8

2.3 The dyadic tree prior π(T )

  • The tree prior recursively creates binary partitions, with stopping, split dimensions, and cut locations governing the resulting dyadic partition.Cut locations can be adaptive, and the cut-point distribution may be continuous or discrete.Source: §2.3 · pp. 8–9

2.4 Knot selection

  • Knot sets should remain relatively small to keep the linear systems for multi-scale predictive and residual processes computationally efficient.Maximin distance design is suggested as a starting point; optimizing knot selection by risk functions is left for future work because the search may be computationally burdensome.Source: §2.4 · p. 10

3 Posterior Inference

Posterior inference uses recursive likelihood calculations and sequential Monte Carlo to integrate or sample over tree structures, while predictive distributions are computed recursively. The procedures support scalable inference and valid joint predictive uncertainty quantification.

3.1 Likelihood computation through recursion

  • After integrating out the partition-tree prior, the marginal likelihood is computed by a bottom-up recursion, with leaf-node likelihood q0(A) as its base case.The recursion offers an alternative to MCMC-based posterior computation.Source: §3.1 · p. 10
  • For fixed tree depth D, one likelihood evaluation costs O(nr^2D^2 + nrdD), scaling linearly in sample size; sequential data additions require updates only at affected nodes.The recursion integrates out the tree prior, and the update property supports online streaming data and sequential design.Source: §3.1 · pp. 10–11

3.2 Sequential Monte Carlo tree-sampling algorithms

  • A one-step look-ahead sequential Monte Carlo algorithm samples tree structures using sequential target distributions and particle-based proposals.Path sharing evaluates node information once across particles, with likelihood-evaluation cost O(nr^2D^2 + nrdD) and memory O(NnrD).Source: §3.2 · pp. 11–12

3.3 Predictive inference

  • ResTGP provides a valid joint predictive distribution over multiple locations, supporting legitimate predictive uncertainty quantification.For a fixed tree, prediction combines contributions along the location’s ancestor nodes with the residual process at its leaf; the residual predictive mean can also be marginalized over trees recursively.Source: §3.3 · p. 13

4 Theoretical results

The paper establishes posterior consistency for ResTGP under stated regularity conditions and examines consistency for continuous regression functions. The results cover fixed and random designs, with a stronger loss guarantee under a uniform-envelope-conditioned prior.

4.2 General consistency conditions

  • C1–C3 provide sufficient conditions for almost-sure posterior consistency, combining a posterior-denominator bound, exponential prior-mass control outside sieves, and testing conditions.The general criterion does not require independence between experiments.Source: §4.2 · p. 15

4.3 Assumptions

  • The consistency results assume a continuous true regression function, specified fixed or random designs, and a Bayesian CART tree prior with bounded, node-specific knot sets.The stated results also impose conditions on the covariance-parameter prior and measurement-error variance prior.Source: §4.3 · p. 15

4.4 Posterior consistency

  • Under Assumptions 1–4, posterior distributions for the regression function and variance parameter converge to the truth under their corresponding topologies as sample size grows.The paper establishes almost-sure posterior consistency for the white-noise leaf model.Source: §4 · p. 13, §4.4 · p. 16
  • Under the model-induced prior, functions converge in measure under the empirical input distribution for fixed design and under the population input distribution for random design.For each fixed error tolerance, posterior mass concentrates on functions exceeding that tolerance on a vanishing fraction of inputs or a set of vanishing population probability.Source: §4.4 · p. 16
  • Conditioning the prior on ∥g∥∞≤M yields a stronger loss guarantee when ∥g0∥∞<M; under random design, a density bounded below away from zero permits the same conclusion in L1(dx).The result is stated in Corollary 3.Source: §4.4 · p. 16

5 Numerical Illustrations

The numerical illustrations compare ResTGP variants with existing methods on likelihood approximation, input-dependent noise, and higher-dimensional prediction. ResTGP models show advantages in several settings, with performance depending on the leaf model and inference approach.

5.1 Likelihood approximation

  • The likelihood experiment compares methods across d = 2–20 and different sample sizes, using maximin designs and r = 30 knots or nearest neighbors.ResTGP-Full is fastest but can be sensitive to its prespecified tree; PP is slower because it uses SMC for tree inference, and timings across implementations do not necessarily reflect theoretical complexity.Source: §5.1 · p. 17
  • For d ≤ 5, all three ResTGP variants approximate the true likelihood well; at d = 10 or 20, ResTGP-PP and ResTGP-Full outperform sSGV with r = 30.At high dimensions, log-likelihood percentage differences range from 5% to 25% depending on n; increasing r from 30 to 100 reduces the reported difference from 20% to 12% for PP and Full.Source: §5.1 · p. 17, Figure 2 · p. 18

5.2 Adaptation to input dependent noise

  • In the 320-observation noise experiment, all methods give satisfactory predictive means, but their ability to capture heterogeneous noise differs.The simulation evaluates held-out predictions using CRPS, log-score, RMSE, and 95% interval coverage.Source: §5.2 · p. 18
  • With input-dependent noise, WN gives the best CRPS and log-score, while PP and sSGV do not adapt to the heterogeneous noise structure.The conclusion also holds for n = 3,200, where WN outperforms other methods, including BART, on CRPS and log-score.Source: §5.2 · pp. 18–19, Figure 3 · p. 19

5.3 Higher-dimensional examples

  • For n = 10,000 and d = 10 or 20, ResTGP-WN alone shows input-dependent predictive standard deviations among the compared models.BART has nearly flat predictive standard deviations, while sSGV, ResTGP-Full, and ResTGP-PP fail to show adaptation to the input-dependent noise.Source: §5.3 · pp. 19–20
  • For d = 10 and 20, ResTGP models predict better than BART and sSGV; PP has the smallest RMSE, while WN with PMMH leads on some distributional scores.PP gives the best CRPS for d = 10; for d = 20, WN with PMMH outperforms other methods in CRPS and log-score, and PMMH costs substantially more computing time than MAP.Source: Table 1 · p. 20, §5.3 · p. 20

6 Application to Storm Surge

The storm-surge application evaluates ResTGP variants against sSGV and BART on held-out ADCIRC simulations at one coastal location. ResTGP improves predictive performance over sSGV, with PP and WN leading on different predictive scores.

6 Application to Storm Surge

  • The comparison uses 44,801 training and 11,476 test simulations, with seven rescaled inputs and storm surges at a single coastal location.The full simulator design comprises 56,277 runs; generalization to massive spatial locations across the domain is investigated elsewhere.Source: §6 · p. 21
  • Exploratory analysis indicates varying spatial-dependence ranges across input parameters, modest variance variation, and response trends that may vary with landfall location.The response also includes very small values near zero across input settings, suggesting a regime change.Source: §6 · p. 21
  • On 11,476 held-out simulations, ResTGP improves predictive performance over sSGV; PP has the smallest RMSE, while WN has the smallest CRPS and log-score.ResTGP takes more computing time because it uses SMC for tree inference.Source: §6 · pp. 21–22, Table 2 · p. 22
  • WN shows more heterogeneous predictive uncertainty than PP, a feature relevant to using spatially heterogeneous extreme-event tail probabilities in risk decisions.PP has the smallest predictive standard deviations; the difference partly explains WN's lower log-score but similar CRPS to PP.Source: §6 · p. 22

7 Concluding Remarks

The paper presents ResTGP as a multiscale tree-based GP framework for large datasets over multidimensional domains. It reports scalable likelihood computation, a Bayesian tree-search strategy, posterior consistency, and future extensions.

7 Concluding Remarks

  • ResTGP decomposes GPs across tree resolutions and adapts local information along the tree, with PP and WN as two leaf-model choices.The likelihood is computed recursively with cost linear in the number of tree nodes given the tree; SMC supports Bayesian exploration of tree space.Source: §7 · p. 22
  • As a nonparametric prior, ResTGP has almost-sure posterior consistency for estimating continuous functions under several loss functions.The paper states this provides enough support for predicting the function class.Source: §7 · p. 22
  • Multiscale modeling across recursive partitions addresses long-range dependence between partition regions while supporting spatial adaptation and computational scalability.The framework integrates GP and tree-based methods for higher-dimensional spaces.Source: §7 · p. 22
  • Future directions include multivariate and space-time dependence, non-axis-aligned partitioning, and ensembles of additive trees.The current framework uses a single tree and is described as performing satisfactorily; ensembles could further improve predictive accuracy.Source: §7 · pp. 22–23

Supplementary Material

The supplementary material contains additional details and R code to reproduce the results.

Supplementary Material

Supplementary Material: A Residual Tree Gaussian Process Modeling Framework for High-Dimensional Data

The supplement collects posterior-inference details, additional numerical results, technical results and proofs, and validation results.

Supplementary Material

S.1 Posterior Inference

Posterior inference targets the tree and global covariance parameters; the supplement describes SMC-based empirical Bayes estimation and a PMMH sampler for fully Bayesian inference.

S.1 Posterior Inference

  • The posterior target is the joint distribution of the tree T and global covariance parameters θ given the data.Source: §S.1 · p. 28
  • The SMC algorithm estimates the marginal tree posterior p(T) when θ is given and can be used for empirical Bayes estimation.Source: §S.1 · p. 28
  • A particle marginal Metropolis–Hastings sampler performs fully Bayesian inference for both T and θ.Source: §S.1 · p. 28

S.1.1 Sequential Monte Carlo algorithm

The SMC implementation shares computations along common tree paths and uses staged tree expansion and resampling. Its computational costs and predictive-performance trade-offs depend on the cut rule and implementation settings.

S.1.1 Sequential Monte Carlo algorithm

  • The implementation identifies node information by its full ancestor path, computes it once, and lets particles reuse it while growing trees with an SMC look-ahead strategy.Particles that agree along a path share residuals and covariance factors there, avoiding duplicate numerical work while retaining the particle population and weights.Source: §S.1.1 · pp. 28–29
  • SMC has shared computation, particle propagation and weighting, and tempered resampling stages; the one-step posterior proposal gives every sampled decision the same incremental weight.Source: §S.1.1 · pp. 28, 30
  • For fixed N, r, d, M, and D, the displayed PP and midpoint-WN numerical costs are linear in n; balanced growth yields at most O(n log^2 n) when other factors are fixed.Fresh sorting adds an n log n term for adaptive WN; for a fixed tree, cached likelihood computation costs O(nrdD + nr^2D^2).Source: §S.1.1 · p. 31
  • Peak cache memory and computational advantage depend on visited ancestral states, realized depth, candidate construction, and cache policy.Releasing obsolete caches bounds active numerical cache memory by O(max_{1≤t≤D}{γ_t n r t + |P_t|r^2t}), excluding temporary candidate buffers and retained tree records.Source: §S.1.1 · p. 31
  • Adaptive cuts often give better predictive performance, but ResTGP-PP with adaptive cuts is significantly more expensive than with middle cuts; ResTGP-WN remains scalable.Source: §S.1.1 · p. 32

S.1.2 Particle marginal Metropolis-Hastings sampler

The sampler section contrasts faster MLE or MAP estimation with fully Bayesian PMMH and lists default settings used in the numerical examples.

S.1.2 Particle marginal Metropolis-Hastings sampler

  • MLE or MAP is computationally faster than PMMH, but the authors report that it may sacrifice statistical accuracy on some occasions.Source: §S.1.2 · p. 34
  • The nominal tree depth is D = ⌈log2(n/r)⌉ and the default particle size is N = 100D.The realized SMC depth may differ for PP and WN, while the full leaf model always has maximum depth D.Source: §S.1.2 · p. 34
  • The numerical examples use r = 30 knots, selected by a maximin design sequentially from root to leaf.Source: §S.1.2 · p. 34

S.2.1 Likelihood approximation

With 100 knots or nearest neighbors, both ResTGP variants consistently outperform sSGV in likelihood approximation, particularly at higher input dimensions; ResTGP-Full stands out in timing.

S.2.1 Likelihood approximation

  • ResTGP-Full sets maximum tree depth when leaf observations fall below 30, while ResTGP-PP defaults to 100⌈log2(n/r)⌉ particles for SMC tree inference.ResTGP-PP uses default hyperparameters.Source: §S.2.1 · p. 35
  • With r=100, ResTGP-Full and ResTGP-PP consistently outperform sSGV in likelihood approximation, especially at higher input dimensions.Figure S.1 compares percentage log-likelihood error against exact GP across input dimensions and sample sizes.Source: §S.2.1 · p. 35, Figure S.1 · p. 36
  • ResTGP-Full and ResTGP-PP have similar approximation accuracy despite using full versus predictive-process residual covariance at the leaves.The authors attribute this partly to tree-dependent knot selection during SMC inference and partly to near-equivalence of the residual structures.Source: §S.2.1 · p. 35
  • Timing results make ResTGP-Full stand out among the compared methods.The figure’s bottom panels show timing; the passage does not report specific timing values.Source: §S.2.1 · p. 35, Figure S.1 · p. 36

S.2.2 Adaptation to input-dependent noise

This section reports ResTGP-WN’s adaptation to input-dependent noise and predictive metrics, alongside high-dimensional timing, diagnostics, and supporting Gaussian-process decompositions. It also presents proof-side calculations and regularity checks, verifying conditions for posterior-consistency results and noting that the criterion does not require iid observations or independence across experiments.

S.2.2.1 Results under n = 320

  • With n = 320, ResTGP-WN has the best CRPS and log-score, ResTGP-PP the smallest RMSE, and tgp coverage closest to nominal.ResTGP-WN timing is competitive with SGV on this small dataset, while ResTGP-PP with adaptive cut is slowest.Source: §S.2.2.1 · p. 37, Table S.1 · p. 37

S.2.2.2 Results under n = 3200

  • With n = 3200, ResTGP models outperform SGV and BART on CRPS, log-score, and CVG; ResTGP-WN leads the models on those metrics and better captures heterogeneity.The comparison excludes tgp because of excessive computing time; Figure S.2 shows predictive means, intervals, and standard deviations across 100 test locations.Source: §S.2.2.2 · pp. 37–38, Table S.2 · p. 38, Figure S.2 · p. 39

S.3.1 Methods settings

  • For the high-dimensional timing comparison, models were fitted on a 10-core, 64GB machine, with matched knot and neighbor settings r = 60 for sSGV and ResTGP.ResTGP-PP used middle cut, ResTGP-WN used uniform cut, and ResTGP-WN PMMH used 3,000 burn-in iterations and 2,000 posterior draws.Source: §S.2.3 · p. 40, §S.3.1 · p. 44

S.2.3 High-dimensional examples

  • In high-dimensional examples, sSGV is fastest and ResTGP-WN with PMMH is slowest; optimization-based ResTGP-WN takes 8 to 23 minutes as d increases from 10 to 20.ResTGP-WN and ResTGP-PP become more expensive as d increases because of particle counts and tree-structure inference; ESS varies across tree levels, with resampling restoring ESS and log evidence stabilizing during expansion.Source: §S.2.3 · p. 40, Figure S.5 · p. 43

S.3.2 Exploratory analysis

  • The supplementary exploratory analysis describes the inputs used in the storm surge application and provides storm-characteristics parameters.The supplied passage does not state the parameter values.Source: §S.3.2 · p. 44, Table S.3 · p. 44

S.3.3 Additional results

  • Additional storm-surge results show predictive means against held-out observations and predictive standard deviations against each input under sSGV, ResTGP-PP, and ResTGP-WN.Within each panel, dots show standard deviations for one variable as other variables vary.Source: §S.3.3 · p. 44, Figure S.7 · p. 45, Figure S.8 · p. 46

S.4.1 Ancillary results

  • For a partitioned Gaussian vector, successive conditional residuals are independent, factorizing its density into lower-dimensional computations that can be extended to k subsets.The decomposition repeatedly regresses one subset on another and treats the resulting residuals as new data.Source: §S.4.1 · pp. 47–48

S.4.2 Proof of Theorem 1

  • Theorem 1’s proof extends the density factorization inductively by splitting a node into a residual region and a child node, whose processes are independent by the spatial BLUP property.The proof uses Stein (1999, Section 1.5).Source: §S.4.2 · p. 49

S.4.4 Proof of Theorem 3

  • The posterior predictive distribution of ξA(x0) is Student-t with 2ã degrees of freedom; adding hT(x0) yields the stated result.Source: §S.4.4 · p. 52

S.5.1 Model properties

  • The residual tree GP has an equivalent finite-dimensional representation with q≤(r+1)K, and its tree is finite almost surely with an exponential tail bound on node count.Proposition P2 gives P(N>u)≤e^(2m0−1−cT u), with cT>0.Source: §S.5.1 · pp. 53–54
  • For each tree, the realized function varies continuously in (θ,z) under the uniform norm, and the prior-to-function map is Borel.The covariance is differentiable in θ with continuous derivatives, supporting the stated continuity properties.Source: §S.5.1 · pp. 55–57

S.5.2 Verification of Condition C1

  • Condition C1 holds under either design for sufficiently small δ, with ε=C*δ when nε^2≥1.The prior assigns uniformly positive mass to δ-neighborhoods of the truth, and the proof establishes the required likelihood conditions.Source: §S.5.2 · pp. 58–59

S.5.3 Verification of Condition C2

  • Proposition P6 considers 0 < α ≤ αmax, where αmax = min{1, (4r̃ log 2)^−1}; θ and v range over Assumption 4 domains, with fixed covariance coordinates held at prescribed values, for sufficiently large n.Source: §S.5.3 · p. 59

S.5.4 Verification of Condition C3

  • For every fixed 0<ε≤1, Proposition P7 gives an αε≤αmax and deterministic Borel tests for Condition C3.The proof constructs tests using metric-entropy bounds for the sieve.Source: §S.5.4 · pp. 60–61
  • The Condition C3 covering bound is uniform in the input probability measure, and fixed covariance coordinates require only one grid point.The proof constructs tests from a finite cover, with local test errors bounded by e^(-nϵ^2/8).Source: §S.5.4 · p. 61

S.5.5 Proof of the criterion in Section 4.2

  • The criterion applies to general experiments without assuming iid observations or independence across experiments.Source: §S.5.5 · p. 62

S.5.6 Proof of Theorem 4

  • Applying the criterion with the prior and sampling densities from Section 4 verifies Theorem 4’s conclusion.Source: §S.5.6 · p. 62

S.5.8 Proof of Corollary 3

  • For Corollary 3, prior-independent moment bounds establish C1, the tests remain unchanged, and boundedness of g converts d_P control to L1(P) control.On the event ∥g∥∞ ≤ M, the proof obtains ∥g − g0∥L1(P) ≤ (1 + 2M)d_P(g, g0).Source: §S.5.8 · p. 63
Loading 2610.02893v1…