PaperlayerView exact revision on arXiv

Source-linked AI summary

Practical heteroskedastic Gaussian process modeling for large simulation experiments

Mickael Binois, Robert B. Gramacy, Michael Ludkovski

arXiv:1611.05902v2stat.MEstat.CO

TL;DR

Large simulation efforts face intrinsic stochasticity and signal-to-noise ratios that can vary across experiments. The paper introduces computational shortcuts and a framework for dependent noise, while retaining complexity comparable to constant-noise Gaussian-process methods.

  • Problem

    Modeling methodology for large simulation efforts with intrinsic stochasticity is lagging, while signal-to-noise ratios can vary dramatically across experiments.

  • Method

    The approach introduces the Woodbury trick for computational shortcuts without approximation and addresses dependent noise within the same framework.

  • Results

    The dependent-noise approach retains computational complexity on par with the best implementations of constant-noise Gaussian-process methods.

  • Takeaways & Limitations

    The framework is illustrated on simulation experiments involving industrial design and epidemic management.

  • Takeaways & Limitations

    The paper identifies shortcomings in SK, including inference for unknowns that is not based completely on the likelihood and requires an additional crutch.

Abstract

from arXiv · show

We present a unified view of likelihood based Gaussian progress regression for simulation experiments exhibiting input-dependent noise. Replication plays an important role in that context, however previous methods leveraging replicates have either ignored the computational savings that come from such design, or have short-cut full likelihood-based inference to remain tractable. Starting with homoskedastic processes, we show how multiple applications of a well-known Woodbury identity facilitate inference for all parameters under the likelihood (without approximation), bypassing the typical full-data sized calculations. We then borrow a latent-variable idea from machine learning to address heteroskedasticity, adapting it to work within the same thrifty inferential framework, thereby simultaneously leveraging the computational and statistical efficiency of designs with replication. The result is an inferential scheme that can be characterized as single objective function, complete with closed form derivatives, for rapid library-based optimization. Illustrations are provided, including real-world simulation experiments from manufacturing and the management of epidemics.

1 Introduction

The paper develops likelihood-based Gaussian process methods for replicated simulation experiments with input-dependent noise, addressing computational and inferential limitations of earlier approaches. It combines Woodbury-based reductions with a jointly smoothed latent noise-field model and full likelihood inference.

  • Motivation: Replication reveals pure simulation variance, helps detect weak signals in noisy settings, and can reduce computation through pre-averaging.
  • Limitations of prior methods: Existing stochastic kriging methods use incomplete likelihood inference, require replication at each design site, and treat input-dependent noise as secondary.
  • Computational approach: The method applies Woodbury identities to reduce the cost of likelihood, derivative, inference, prediction, and classical kriging calculations under replication.
  • Heteroskedastic modeling: It extends this framework to heteroskedastic modeling by jointly inferring a spatial noise-field and the mean-field Gaussian process.
  • Heteroskedastic modeling: A Gaussian process prior smooths the latent log-variance process while preserving full derivatives for optimizing coupled GP hyperparameters with library routines.
  • Optimization: The smoothing mechanism does not bias the stationary-point predictor, instead annealing the likelihood and easing initialization to accelerate convergence.

2 GPs under replication and heteroskedasticity

Gaussian-process models represent noisy simulation outputs through a covariance structure and can accommodate input-dependent noise. Replication reduces the effective design size, but prior approaches often separate variance estimation from mean modeling or impose restrictive replication requirements.

  • The observation model decomposes each output into a latent function and Gaussian noise with variance r(x), which may be constant or input-dependent.
  • Standard GP prediction and likelihood inference require covariance-matrix inversions, determinants, and repeated O(N^3) decompositions.This cubic cost limits the experiment sizes that can be modeled directly.
  • Replication: Replication replaces N observations with n unique input locations, where n ≪ N, while allowing site-specific replicate counts and unreplicated sites.Replicate averages and the associated design structure provide the basis for reduced computations.
  • Prior approaches: Earlier stochastic-kriging approaches estimated replicate noise separately from other kernel parameters and evaluated a pre-averaged likelihood at O(n^3) cost.This treatment can require more than one replicate at every input and does not jointly model the mean and variance processes.
  • Proposed framework: The paper combines a GP mean with GP latent log variances and replication-aware Woodbury calculations to target joint, likelihood-based inference with O(n^3) demands.The stated goal is to retain smoothing and full likelihood inference while exploiting replication.
  • Prior approaches: Latent-variable GP methods jointly model the mean and noise processes, but previous replication-aware variants approximated full-data quantities or maximized latent log variances rather than averaging over them.The latter produces a solution that is not a smooth GP realization.

3 Fast GP inference and prediction under replication

Replication lets GP prediction and likelihood calculations be rewritten using matrices indexed by unique design sites. Woodbury identities preserve full-data results while reducing the dominant computational work to the unique-n scale.

  • Woodbury reduction: Woodbury identities map full-N GP prediction equations to unique-n matrices under replication.The construction uses replicated design blocks and diagonal replicate-count matrices.
  • Computational savings: The full-N matrices need not be constructed because their inverse and determinant calculations exploit the low-dimensional replication structure.This avoids direct operations on the full observation dimension.
  • Prediction: The unique-n predictive mean is identical to the full-N predictive mean, and the same identity holds for predictive variance.These results retain the full-data predictor's BLUP properties after conditioning on hyperparameters.
  • Computational savings: O(n^3) decomposition cost plus O(N+n) or O(n) extra terms yields essentially O(n^3) inference after the unique-n matrices are decomposed.The additional costs are negligible relative to the cubic unique-site decomposition.
  • Likelihood inference: The full-data estimate of the common variance factor corrects the variance contribution from replicate-level deviations.As replicate counts increase, the estimator converges to a replicate-based variance expression.
  • Optimization: Closed-form derivatives support numerical optimization, while the same calculations apply to homoskedastic and heteroskedastic models.The heteroskedastic scheme estimates latent variance parameters at unique design locations.

4 Practical heteroskedastic modeling

The method models input-dependent noise through a latent log-variance process smoothed over unique design locations, while replication enables likelihood inference with fewer latent variables and lower computational cost. Smoothing regularizes variance estimates, supports optimization, and remains compatible with full maximum-likelihood inference.

  • Heteroskedastic GP modeling learns a heterogeneous diagonal noise matrix rather than imposing a common variance across inputs.
  • Replication reduces the latent variance representation from the full data to unique design locations through Woodbury identities.The reduction uses the relationship U⊤Λ_NU = A_nΛ_n.
  • A regularizing GP maps latent variables to smooth, positive variance estimates across the input space.The smoothing construction guarantees smooth and positive variances for every latent-variable setting.
  • The smoothing parameter g acts as an annealing device that improves numerical stability without changing the global optimum.The objective is maximized at g = 0 with Δ_n = log Λ_n, so smoothing is redundant at the global maximum.
  • In the motorcycle illustration, smoothed variances look more sensible than unsmoothed estimates, while the available replication provides only minor Woodbury savings.The data have at most six replicates per input and often only one.
  • The resulting likelihood optimization has complexity O(N + n^3), essentially O(n^3), and requires no minimum replication level.The computational benefit is greatest when N is much larger than n.

5 Empirical comparison

Empirical comparisons show that the proposed methods match existing motorcycle-data accuracy while remaining fast, and that heteroskedastic models perform best in the ATO and epidemics experiments. Replication-aware fitting also yields substantial computational savings and more accurate variance surfaces than direct empirical-variance interpolation.

  • 5.1 Benchmark data: The WGP and WHGP methods achieve accuracy comparable to published motorcycle-data comparators, while WHGP completes 100 optimization iterations in under one second.
  • 5.2 Assemble to order: On ATO data, heteroskedastic comparators are the clear winners across the evaluated predictive-score comparisons.Scores combine predictions of the mean and variance across testing, training-site, combined, and in-sample settings.
  • 5.2 Assemble to order: The Woodbury-based heteroskedastic comparator is about 40× faster than the corresponding full-data homoskedastic comparator.Here the unique-site count is approximately one fifth of the full data size.
  • 5.2 Assemble to order: Running heteroskedastic optimization to convergence takes over 50× longer but produces minimal score improvement.The comparison is between truncated and fully optimized heteroskedastic comparators.
  • 5.3 Epidemics management: In the epidemics example, the proposed method recovers key variance-surface features with substantially fewer replicates than the reference design.Direct interpolation of empirical variances is reported to be much less accurate.
  • 5.3 Epidemics management: The proposed method outperforms empirical variance estimates even with many replicates, while empirical estimates are notably high-biased when a_i < 50.Stochastic Kriging mean predictions remain off even at inputs with a_i = 50 or a_i = 100.

6 Discussion

The framework combines computational shortcuts from Woodbury-based GP inference with a smoothed latent variance process for practical heteroskedastic modeling. It preserves input-dependent noise modeling while retaining computational complexity comparable to efficient constant-noise GP implementations.

  • The framework combines Woodbury computational shortcuts with a smoothed latent process for heteroskedastic GP inference.
  • Its latent-process construction replaces more expensive EM and MCMC simulation-based inference alternatives.
  • The resulting framework improves modeling fidelity by representing input-dependent noise while retaining complexity comparable to efficient constant-noise GP methods.
  • The approach also performs well under other forms of non-stationarity despite model mismatch.
  • Compared with free pseudo-input selection, fixing latent locations to replicated design sites restricts applicability when replication is infeasible or undesired.Free pseudo-input methods can provide additional computational savings, but rely on approximation and efficient pseudo-input selection.

A.1 Derivations supporting lemmas in Section 3

Under replication, matrix products involving the replication indicator and covariance matrices reduce full-data calculations to unique-site quantities. Woodbury’s identity then supplies the corresponding inverse and determinant expressions needed for likelihood inference and prediction.

  • Replication induces identities linking full-data matrices to unique-site matrices through U and the replication-count matrix A_n.
  • Woodbury’s identity provides the reduced expression for the covariance calculation under replication.
  • The same reduced calculations apply when the averaged response is replaced by the kernel vector for prediction.
  • Using the full response vector on both sides yields the corresponding full-data quadratic-form expression.
  • The determinant term admits an analogous replication-aware reduction.

A.2 Empirical demonstration of Woodbury speedups

The experiment compares homoskedastic GP inference with and without exploiting replication through the Woodbury trick. Replication-aware computation produces essentially identical estimates while substantially reducing runtime.

  • Runtime comparison: More than 4000 times faster, the Woodbury-trick version reduces computation when replicated inputs are represented by unique locations and multiplicities.The comparison bypasses the Woodbury trick in the slower implementation by treating the data as having no replicates.
  • Inference agreement: 1.099367 and 1.789508 versus 1.099308 and 1.789700 are the paired estimated lengthscales from the replication-aware and full-data implementations.The code checks that both methods perform the same inference by printing estimated lengthscales.
  • Design representation: The package accepts unique locations, responses, and replication counts, or detects replicates internally from a supplied design matrix.The internal find reps function builds the representation used for further calculations.
  • Latent-variable extension: The accompanying derivations and examples examine latent-variable optimization and empirically demonstrate consistency of hetGP optimization.The section proceeds from supporting derivations to a heteroskedastic example and an empirical consistency demonstration.

B.1 Proof of Lemma 4.1

The proof analyzes the likelihood at the proposed latent-variable maximizer by separating terms that remain invariant from terms that depend on the latent variance.

  • Proof of Lemma 4.1: The top part of Eq. (18) is invariant whenever Λ_n = Λ*ₙ, so the proof focuses on the remaining likelihood terms.The argument identifies the invariant portion before analyzing the latent-variance-dependent terms.
  • Proof of Lemma 4.1: log(ν(g)) and log |Υ(g)| are the remaining terms analyzed for their dependence on the latent variable g.The proof introduces these terms explicitly as the focus of the derivative argument.
  • Proof of Lemma 4.1: The derivative of log |Υ(g)| is positive because it is the trace of a positive definite matrix, so log |Υ(g)| increases with g.This term does not depend on Δ_n.

B.2 Illustration in hetGP

The hetGP illustration fits a heteroskedastic GP to motorcycle acceleration data and uses replication-aware prediction and visualization.

  • Examples: The example reuses the motorcycle data, with code for additional ATO and SIR Monte Carlo experiments available upon request.The motorcycle results correspond to the bottom row of Figure 1.
  • Motorcycle example: About one second is required to optimize the unknown parameters, including latent variances, with R’s default linear algebra libraries.The fit uses mleHetGP with the motorcycle times and acceleration data.
  • Replication handling: Raw times and acceleration responses are supplied directly, while find reps identifies the design’s replication for Woodbury calculations.The package navigates the small amount of replication through its internal representation.
  • Prediction visualization: The prediction code returns means, uncertainty bounds, and predictive noise variance across a 301-point time grid.The plots show acceleration observations with the fitted mean and quantile curves, alongside estimated noise variance.

B.3 Robustness in numerical optimization

The paper evaluates initialization for latent-noise optimization on motorcycle and two-dimensional Branin examples. Thoughtful initialization consistently outperforms random initialization across the reported diagnostics.

  • Robustness in numerical optimization: Figure 5 compares one thousand random starting points with the proposed starting point for motorcycle and 2d Branin data.The Branin score uses a 51 × 51 grid with 10 replicates per testing location.
  • Figure 5 comparison: The thoughtful-initialization hetGP fit is shown in solid red, while the homoskedastic fit is shown as a blue dotted reference.These encodings apply across all four panels.
  • Robustness in numerical optimization: More than 99% of the time, thoughtful initialization yields better estimates than random initialization.Figure 5 compares random starts with the proposed start using likelihoods, latent variables, and proper scores.
  • Figure 5 comparison: Random initialization sometimes produces egregiously over-estimated latent variables in the motorcycle example.The top-right panel displays the latent-variable estimates in log space.

C Non-stationarity

Input-dependent noise is difficult to distinguish from mean non-stationarity when data are sparse, but a heteroskedastic-only model can still be useful for some tasks. In the examples, estimated noise increases where mean dynamics change most rapidly, supporting improved accuracy and active-learning heuristics, while denser design or more replication can return a homogeneous model.

  • Input-dependent noise is difficult to distinguish from mean non-stationarity when little data are available.
  • A mis-specified heteroskedastic-only model may be better than a purely stationary model for many tasks.
  • In both examples, estimated noise increases in regions where mean dynamics change most rapidly.
  • This estimated-noise behavior improves accuracy and can facilitate active-learning heuristics in sequential design.
  • Focusing design effort in high-noise regions can be useful even when rapidly changing signal, rather than high noise, is the true explanation.
  • With further exploration or more replicates in the interesting region, the behavior disappears and a homogeneous model is returned.
Loading 1611.05902v2…