PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian treed Gaussian process models with an application to computer modeling

Robert B. Gramacy, Herbert K. H. Lee

arXiv:0710.4536v10stat.MEstat.APstat.CO

TL;DR

The paper develops a Bayesian Gaussian-process expansion for modeling complex simulator behavior with nonstationary structure. It combines stationary Gaussian processes with treed partitioning and reports effective Bayesian model averaging, smooth fitted functions, and a semiparametric, tractable, highly accurate model.

  • Problem

    Increasingly sophisticated simulators create modeling challenges, while standard methods have disadvantages and may be tractable only for relatively small datasets.

  • Method

    The approach expands Gaussian processes using Bayesian modeling, stationary process components, treed partitioning, and a tree prior favoring minimal depth.

  • Results

    Bayesian model averaging yields effective fits across several examples, with mean fitted functions typically quite smooth in practice and good behavior for large rocket-booster simulator regions.

  • Takeaways & Limitations

    The resulting Bayesian treed model is semiparametric, tractable, and highly accurate, and is demonstrated beyond the rocket-booster application.

  • Takeaways & Limitations

    Computational constraints include matrix operations such as multiplications involving K^-1, with factorization assumed available when computing K^-1.

Abstract

from arXiv · show

Motivated by a computer experiment for the design of a rocket booster, this paper explores nonstationary modeling methodologies that couple stationary Gaussian processes with treed partitioning. Partitioning is a simple but effective method for dealing with nonstationarity. The methodological developments and statistical computing details which make this approach efficient are described in detail. In addition to providing an analysis of the rocket booster simulator, our approach is demonstrated to be effective in other arenas.

1 Introduction

The paper develops treed Gaussian processes to address nonstationarity, heteroscedasticity, and computational scale in emulating an expensive rocket-booster simulator. Bayesian partitioning combines local stationary GP models with model averaging to provide flexible predictive uncertainty and tractable inference.

  • Motivation: NASA’s simulator typically required five to twenty hours per run, motivating a statistical emulator or surrogate model.The emulator was intended to model simulator behavior more efficiently than repeatedly running the simulator.
  • Limitations of standard GPs: The standard stationary smooth GP approach was inadequate for NASA data because of nonstationarity, heteroscedasticity, and dataset size.These limitations motivated an expansion beyond the standard GP formulation.
  • Limitations of standard GPs: Standard GP inference typically requires computing time in O(N3) to calculate inverses of N × N covariance matrices.The computational burden grows with the number of data points.
  • Limitations of standard GPs: Stationary GPs use one covariance structure across the input space, which may be too strong when regional behavior differs.The paper also notes that stationary predictive error does not directly depend on locally observed response values.
  • Treed Gaussian processes: Partitioning fits separate stationary GP models within regions, creating a straightforward nonstationary model while reducing computational demands through smaller datasets.Bayesian model averaging allows predictive uncertainty to vary beyond the constraints of a stationary model.
  • Results and implications: Bayesian model averaging typically produces smooth mean fitted functions, while retaining the ability to represent discontinuities when supported by the data.The paper combines stationary GPs and treed partitioning into a tractable model for nonparametric regression.

2 The Langley Glide-Back Booster Simulation

The rocket-booster simulator models lift across speed, angle of attack, and quantized sideslip levels, but costly and unreliable numerical convergence complicates surrogate modeling. Sharp regime changes and suspected simulator artifacts motivate smoothing, nonstationary modeling, and spatially varying uncertainty estimates.

  • Simulation design: The Langley Glide-Back Booster simulator models lift as a function of Mach speed, angle of attack, and six quantized sideslip levels.Speed spans Mach 0–6, angle of attack spans −5° to 30°, and the simulator was run at 3041 locations.
  • Observed response surface: The lift surface has a sharp ridge near Mach one, marking a difficult transition between subsonic and supersonic regimes.The discontinuity is especially visible at larger angles of attack, making continuity versus discontinuity a modeling choice.
  • Modeling requirements: A stationary Gaussian process is inadequate because the surface is nonsmooth and prediction uncertainty is not constant across the input space.The desired model must remain computationally feasible while providing local uncertainty estimates.
  • Observed response surface: False convergence produces isolated spikes, noisy depressions, and deviations that may be numerical rather than physical.The authors therefore want the surrogate to smooth suspect points while retaining uncertainty about whether some high-speed behavior is genuinely unstable.
  • Modeling requirements: Engineers need both the mean lift surface and associated predictive uncertainty for design diagnosis, trajectory planning, and contingency planning.These uncertainties matter because they vary across the surface.

3 Related work

The paper situates its approach among Gaussian-process and Bayesian partition methods for nonstationary spatial modeling. It combines stationary GPs with treed partitions to gain local flexibility while retaining interpretable structure and computational advantages.

  • Core approach: The proposed framework combines standard Gaussian processes with treed partitioning in a Bayesian hierarchical model with model averaging.The approach targets nonparametric and semiparametric nonstationary modeling.
  • Gaussian processes: Stationary GP models describe smooth spatial dependence through correlation functions, while nugget terms represent measurement error and improve numerical stability.The nugget is positive and helps prevent the correlation matrix from becoming numerically singular.
  • Partition models: Treed partitioning recursively divides the input space with binary, axis-aligned splits and fits an independent model in each leaf.This structure provides a flexible nonstationary process while allowing minimum data requirements per partition.
  • Treed Gaussian processes: The paper extends Bayesian treed models by fitting stationary Gaussian processes, rather than constant or linear models, within tree leaves.This retains local GP behavior while allowing region-specific modeling.
  • Treed Gaussian processes: Compared with Voronoi tessellations, treed GPs target fewer distinct partitions, producing a simpler and more interpretable overall model.Treed structure is less intricate, whereas Voronoi partitions can be complex and difficult to interpret.
  • Advantages and trade-offs: Unlike approaches that enforce continuity, treed models can represent discontinuities and may reduce computation by inverting smaller correlation matrices.The tree structure introduces flexibility while reducing computational burden relative to a single GP.

4 Treed Gaussian process models

The treed Gaussian process model fits local stationary GPs with linear trends inside recursively defined regions, then averages over tree structures to account for partition uncertainty.

  • Model construction: The model fits stationary Gaussian processes with linear trends independently within each leaf of a recursive partition.This extends earlier treed models that used constant or linear leaf models.
  • Model construction: Reversible-jump MCMC integrates over possible tree structures under a tree-generating prior.The tree is treated as uncertain rather than fixed.
  • Prediction: Prediction is conditioned on each tree and then averaged over the posterior to provide a full accounting of uncertainty.Posterior averaging incorporates uncertainty in the partition structure.

4.1 Hierarchical Model

The hierarchical model assigns region-specific GP parameters within a tree while allowing shared prior structure across regions. It is designed for regime boundaries where discontinuous response surfaces are plausible.

  • Regional structure: A tree recursively partitions the input space into R non-overlapping regions, each containing its own observations and GP component.For region ν, the data are represented as D_ν = {X_ν, Z_ν} with n_ν observations.
  • Regime boundaries: The model permits discontinuous response surfaces at boundaries between physical regimes such as subsonic and supersonic flight.This is treated as a feature for phase-transition problems rather than a defect requiring enforced continuity.
  • Regional structure: The model uses region-specific correlation parameters, including range and nugget terms, with priors that can encode a nonstationary global covariance structure.The correlation family may be isotropic or separable, with a fixed smoothness power and random range and nugget parameters.
  • Prior specification: A mixture-of-Gammas prior for the range parameter supports both wavy and smooth or approximately linear GP surfaces.The mixture places mass on small ranges and on smoother parameterizations.
  • Prior specification: The nugget prior is exponential, while alternative hierarchical specifications can share information across regional parameters.The paper also discusses hierarchical dependence for collections of correlation, variance, and nugget parameters.

4.2 Estimation

Estimation alternates posterior updates for treed GP parameters and tree structures, using Gibbs, Metropolis-Hastings, and reversible-jump MCMC steps. Specialized tree proposals improve exploration while preserving valid partitions and accounting for tree-prior changes.

  • Parameter estimation: Posterior sampling alternates draws of region-specific parameters conditional on the tree with updates of hierarchical prior parameters.Most parameters use Gibbs steps, while covariance parameters require Metropolis-Hastings draws.
  • Tree estimation: RJ-MCMC integrates over uncertainty in the tree structure while proposing changes that alter the partition and parameter-space dimension.Grow and prune operations add or remove regions and therefore require reversible-jump calculations.
  • Tree estimation: Change and swap proposals modify existing split configurations, while grow and prune proposals respectively add and remove partitions with corresponding parameter proposals or absorption.Swap exchanges adjacent parent-child splitting rules, whereas change moves a split point along one input dimension.
  • Tree estimation: Rotations provide more candidate nodes for pruning and help the Markov chain escape local minima in the marginal posterior of T.They address problematic same-variable internal-node configurations where a swap can create an empty child region.
  • Tree estimation: Rotations preserve leaf partitions, so their likelihood ratio is 1; acceptance is driven by the tree prior, which favors minimal depth.Depth changes in two subtrees make the acceptance-ratio calculation non-trivial.
  • Tree estimation: Grow proposals assign one child the parent correlation structure and draw the sibling’s correlation function from the prior to maintain an ergodic, reversible chain.The construction also handles new region parameters and the associated Jacobian term.

4.3 Treed GP Prediction

Prediction is computed conditionally on each tree and Gaussian-process parameter state, then aggregated over posterior samples. This averaging preserves flexibility for discontinuities while smoothing predictions near likely partition boundaries.

  • Conditional prediction: Conditional on a particular tree, the posterior predictive surface is discontinuous across that tree’s partition boundaries.The conditional predictive distribution is specified through its mean and variance under the GP model.
  • Posterior aggregation: Posterior averaging over tree and parameter samples tends to smooth the predictive mean near likely partition boundaries.Tree operations such as grow, prune, change, and swap contribute to this aggregate posterior prediction.
  • Model behavior: When data indicate a non-smooth process, the treed GP retains flexibility to model discontinuities.The model therefore does not require every posterior predictive fit to be continuous across partitions.
  • Model behavior: For data consistent with a continuous process, the treed GP fits can be almost indistinguishable from continuous fits.The motorcycle-data example is cited as an instance of this behavior.

4.4 Implementation

The implementation combines compiled numerical code, posterior-sampling strategies, and parallel prediction to make treed GP computation practical. It also provides an R interface and uses optimized linear-algebra libraries.

  • Software: The treed GP model is implemented in a mixture of C and C++ and exposes an R interface through the tgp package.The package was made available on CRAN.
  • Numerical computation: C code interfaces with platform-specific BLAS/Lapack libraries or ATLAS for the linear algebra required by GP computation.These options support fast execution across multiple platforms.
  • Sampling: Conditioned on the tree, Metropolis-Hastings parameters use uniform sliding-window proposals centered on the last accepted setting.The scheme is illustrated for a proposed nugget parameter.
  • Parallelization: Prediction after conditioning on tree and parameter values is parallelized with a producer-consumer model using PThreads.The implementation is designed to exploit multiple processors.
  • Parallelization: Parallel computation achieves speed-ups of at least a factor of two as multi-processor machines become commonplace.The reported gain concerns the prediction computation after conditioning on (T, θ).

4.5 Illustration

The motorcycle-data illustration shows why a stationary GP is inadequate for input-dependent uncertainty and changing smoothness. The treed GP captures regional uncertainty and smoother local behavior, while typically selecting about three regions and running faster than a GP-expert mixture.

  • Data and motivation: The motorcycle dataset measures head acceleration over time after impact and exhibits input-dependent noise, or heteroscedasticity.The example is used to illustrate a nonstationary covariance structure.
  • Stationary versus treed fits: The stationary GP cannot capture heteroscedasticity, and central-region variability makes both ends of its fit excessively wiggly.The transition from the flat initial region is consequently fitted with an upward curve before descending.
  • Partition structure: The treed GP almost always partitions the motorcycle data into three regions, occasionally four, and rarely two.The reported average is 3.11 partitions.
  • Computational comparison: The treed GP runs in about 3 minutes for 25,000 MCMC rounds, compared with roughly one hour for the GP-expert mixture using 11,000 rounds.The timings were reported on different processor speeds.
  • Posterior averaging: Although individual partition models are typically discontinuous, their averaged fitted function is essentially continuous.The mean fitted function in Figure 3 is described as essentially that of a continuous function.

4.6 Limiting linear models

The model can replace a full Gaussian process with a limiting linear model in partitions where the response is sufficiently simple, reducing computational cost while retaining GPs where curvature is needed.

  • The approach allows different partitions to use either Gaussian processes or simpler linear models.Replacing a full GP with a linear model greatly reduces the resources required to update that partition.
  • A standard linear model is a limiting case of the implemented Gaussian process and is more parsimonious and computationally efficient.A model-switching prior enables practical implementation of this choice.
  • The motorcycle accident example illustrates that the leftmost partition could use a line, the center requires a GP, and the rightmost is mostly linear.This demonstrates how partition-specific model choices match local response complexity.
  • In the motorcycle accident fits, the leftmost partition is nearly always flat, the rightmost is often flat, and the center is typically curved.The center can nevertheless be essentially piecewise linear when the estimated range parameter d is large.
  • The treed and non-treed Gaussian-process models with jumps to the limiting linear model are implemented in the tgp package on CRAN.The analyses use the full formulation described by the authors.

5 Rocket Booster Model Results

The rocket-booster analysis applies treed Gaussian processes through RJ-MCMC to model lift across operating conditions. The model captures smooth behavior and sharp regime changes, while predictive uncertainty reflects both data coverage and numerical instability.

  • Ten independent RJ-MCMC chains were used for the rocket-booster treed GP model, with 15,000 rounds per chain and the first 5,000 discarded as burn-in.Every tenth subsequent draw was treated as a posterior-distribution sample, and 10,000 samples were saved.
  • About 60 hours were required for the rocket-booster computation, whereas a stationary GP would have taken a minimum of 708 hours for the same number of samples.The comparison uses a 3041 × 3041 matrix inversion taking about 17 seconds and is described as a gross underestimate of stationary-GP cost.
  • The treed GP models lift surfaces across six sideslip levels using posterior predictive means and 5%-to-95% predictive-quantile widths.The mean surfaces cover speed and angle of attack, while variance surfaces show experimental-run locations and darker shades for higher values.
  • The treed GP captures mostly smooth lift behavior and the sharp transition at Mach one, while smoothing potential false convergences.Aggressive partitioning occurs near Mach one, with additional partitions addressing a low-speed singularity and numerical instabilities near Mach three.
  • Predictive variability increases where lift changes rapidly and where false numerical convergence is possible, but remains lower near Mach one at sideslip level zero because many samples cover that region.Higher sideslip levels show increased uncertainty near Mach one, and noisy depression data contribute uncertainty around higher angles of attack at sideslip one-half.
  • The fitted mean is essentially continuous even though individual MCMC treed-GP fits often contain partition-induced discontinuities.The smooth mean results from averaging individual components that may not be continuous.

6 Conclusion

The Bayesian treed GP is presented as a fully Bayesian, flexible nonstationary model that combines treed partitioning with GP components and includes linear models as a special case. Its region-specific predictive uncertainty supports active learning and adaptive exploration of computer experiments.

  • 6 Conclusion: The authors developed the treed Gaussian process model for the rocket booster computer experiment and characterize it as a simple and efficient method for nonstationary modeling.They also identify broader uses beyond the booster application.
  • 6 Conclusion: The Bayesian treed GP provides a fully Bayesian approach to nonstationary modeling with a hierarchical correlation-function parameterization that can be replaced by other correlation families.The model treats the correlation function as a modular component within the hierarchical specification.
  • 6 Conclusion: The model is uniquely nonstationary, semiparametric, tractable, and highly accurate, while containing the Bayesian treed linear model as a special case.Its limiting linear-model parameterization remains useful for Bayesian posterior estimation and prediction.
  • 6 Conclusion: Many computer experiments contain substantial linearity, including large regions of the rocket booster simulator's input space.This motivates retaining linear behavior within the broader nonstationary model.
  • 6 Conclusion: The Bayesian treed GP yields a full posterior predictive distribution with a nonstationary, region-specific estimate of predictive variance.This provides uncertainty estimates that vary across regions of the input domain.
  • 6 Conclusion: Region-specific predictive uncertainty can support active learning and an efficient framework for adaptively exploring computer experiment parameters.The authors identify sequential design as a major potential application of these characteristics.
Loading 0710.4536v10…