Source-linked AI summary
Efficient emulators of computer experiments using compactly supported correlation functions, with an application to cosmology
Cari G. Kaufman, Derek Bingham, Salman Habib, Katrin Heitmann, Joshua A. Frieman
TL;DR
Large designs make likelihood calculations for standard strictly positive-correlation GP emulators computationally infeasible. The paper combines compactly supported correlations, regression terms, and range restrictions to build a sparse emulator; in the photometric-redshift application, it reports higher predictive efficiency than the standard model alongside computational efficiency.
Problem
For large designs, the standard GP emulator's likelihood calculations become computationally infeasible because they scale as O(n^3).
Source: §1 · p. 5Method
The method uses anisotropic compactly supported correlations with regression functions in the GP mean and range restrictions that induce sparsity while allowing ranges to trade off.
Source: §2.1 · p. 6, §2.2 · p. 7Results
In the photometric-redshift analysis, the proposed methods were computationally efficient and had higher Nash–Sutcliffe efficiency than the traditional model.
Source: §5 · p. 16Takeaways & Limitations
The application demonstrates fitting the model to 20,000 training points and making predictions at 80,000 validation inputs.
Source: §5.2 · p. 19In the simulation study, prediction intervals were often more conservative than nominal, especially for two-dimensional inputs, and the authors caution that these results may not represent all simulators.
Source: §4 · pp. 14–15
Abstract
from arXiv · showhide
Statistical emulators of computer simulators have proven to be useful in a variety of applications. The widely adopted model for emulator building, using a Gaussian process model with strictly positive correlation function, is computationally intractable when the number of simulator evaluations is large. We propose a new model that uses a combination of low-order regression terms and compactly supported correlation functions to recreate the desired predictive behavior of the emulator at a fraction of the computational cost. Following the usual approach of taking the correlation to be a product of correlations in each input dimension, we show how to impose restrictions on the ranges of the correlations, giving sparsity, while also allowing the ranges to trade off against one another, thereby giving good predictive performance. We illustrate the method using data from a computer simulator of photometric redshift with 20,000 simulator evaluations and 80,000 predictions.
1. Gaussian process models for computer experiments.
The standard GP emulator models simulator output with a mean structure and a covariance determined by a correlation function. For strictly positive correlations, likelihood calculations scale as O(n^3), making the usual model infeasible at large sample sizes.
1. Gaussian process models for computer experiments.
- The model represents simulator output as a GP with fixed regression functions in its mean and a zero-mean GP residual.The covariance is σ^2 times a correlation function; σ^2 is the residual marginal variance and θ controls correlation.Source: §1 · p. 4
- Inference repeatedly evaluates a likelihood based on the observed-response vector, regression matrix, and correlation matrix, while prediction at a new input has a conditional normal distribution.The conditional prediction uses correlations between the new input and observed responses.Source: §1 · pp. 4–5
- The standard specification uses a constant mean and a product of one-dimensional power-exponential correlations, with dimension-specific parameters allowing anisotropy.The correlation parameters satisfy φ_k > 0 and 1 ≤ α_k ≤ 2.Source: §1 · p. 5
- Strictly positive correlations make likelihood operations grow as O(n^3), rendering them infeasible for moderate and impossible for large designs.The costly operations are computing the log determinant of the correlation matrix and solving a linear system involving it.Source: §1 · pp. 5–6
2. Building computationally efficient emulators.
The proposed emulator combines compactly supported correlations, regression terms in the GP mean, and parameter restrictions that preserve computational efficiency. Together, these choices target sparse computation while retaining flexible predictive behavior.
2. Building computationally efficient emulators.
- The method combines compactly supported correlations for small-scale variability, regression terms for large-scale variability, and efficiency-preserving parameter constraints.The paper presents these three innovations as working together to produce a flexible, fast emulator.Source: §2 · p. 6
2.1. Compactly supported correlation functions.
- Product-form compactly supported correlations create covariance-matrix zeros that enable sparse matrix algorithms while retaining dimension-specific ranges for anisotropy.Each one-dimensional correlation is zero at distances at or beyond its cutoff; the cutoffs also control sparsity and must be restricted for computational savings.Source: §2.1 · pp. 6–7
2.2. Regression functions.
- Regression terms model the simulator's large-scale mean structure, potentially shortening residual correlation ranges and allowing predictions similar to a full-correlation model at large-data scale.The paper reports that combining a small basis with compact support yields behavior similar to the true-correlation example, while that full-correlation model cannot be applied to large data sets.Source: §2.2 · pp. 7–8
2.3. Prior distributions.
- The prior on dimension-specific cutoffs constrains matrix sparsity while allowing ranges to trade off, but may not suit unimportant inputs with effectively infinite ranges.The approach assumes such inputs have already been identified and fixed through prior screening; integrating out β and σ^2 also leaves posterior sampling over τ.Source: §2.3 · pp. 8–9
3. Implementation and computational considerations.
The implementation samples range parameters using a sparse representation of the correlation matrix, then computes predictions from selected posterior iterations. In reference-data timing tests, sparse correlations reduced most likelihood calculations by one to three orders of magnitude.
3. Implementation and computational considerations.
- The two-stage workflow samples the marginal posterior of the ranges after integrating out β and σ^2, then estimates predictive means and variances from selected iterations.Conditional means and variances are combined using the laws of iterated expectation.Source: §3 · pp. 9–10
- After each proposed range vector, the algorithm builds a sparse correlation representation and evaluates the likelihood using only nonzero entries and a sparse Cholesky decomposition.The operations include identifying pairs within every dimension's cutoff, computing their product correlations, and decomposing the resulting matrix.Source: §3 · p. 10
- One to three orders of magnitude: sparse correlations reduced the time for most likelihood calculations in the reference-data tests.Figure 2 reports average timings across sample sizes and sparsity levels, with log-scale axes.Source: §3 · p. 10, Figure 2 · p. 11
4. Simulation study.
The simulation study compares the sparse emulator with the standard GP model across simulated process conditions and sample sizes. The sparse method captures similar explained variability in most cases by n = 1,100, while its prediction intervals can be conservative, particularly in two dimensions.
4. Simulation study.
- The study compares predictions under data generated from the standard model across dimensions, correlation settings, sample sizes from 100 to 1,100, and repeated random Latin hypercube designs.Predictions are evaluated at 512 orthogonal array-based Latin hypercube points; the evaluation design was chosen for integral approximation.Source: §4 · pp. 11–12
- The sparse model uses fifth-order Legendre polynomial regression terms, with two-dimensional interactions, and a truncated-power correlation cutoff chosen for 2% or 5% maximum off-diagonal proportions.The authors report that fourth- or fifth-order polynomials suffice for most applications, but do not adapt polynomial order to each dataset in this study.Source: §4 · p. 12
- Across the eight process conditions, standard-model NSE is higher for smoother processes and longer effective ranges, and the sparse method's NSE is higher when 5% rather than 2% of off-diagonal elements are allowed.The difference between the two sparsity settings is minor relative to differences among the processes.Source: §4 · pp. 13–14, Figure 3 · p. 13
- By n = 1,100, the sparse method captures about the same total variability as the standard method in most cases.The comparison uses Nash–Sutcliffe efficiency (NSE) for predictions from the posterior mean.Source: Figure 3 · p. 13, §4 · p. 14
- Proposed-model 95% prediction intervals are often more conservative than nominal, particularly in two dimensions, because sparsity-imposed shorter ranges keep interval widths from shrinking fast enough as n grows.The authors note that coverage in the photometric-redshift application was close to nominal in a held-out sample, with range-parameter samples well away from the sparsity boundary.Source: §4 · pp. 14–15, Figure 4 · p. 14
- The authors caution that the simulation findings should not be taken as representative of all simulators encountered in practice.In the photometric-redshift application, exploratory held-out analysis suggested coverage close to nominal.Source: §4 · p. 15
5. Application to photometric redshifts.
The application uses photometric measurements in four wavebands to predict spectroscopic redshifts, with a preliminary model comparison followed by a full analysis. The final fit predicts 80,000 validation cases from 20,000 training points.
5. Application to photometric redshifts.
- The application predicts spectroscopic redshift from four-band galaxy photometry, an alternative to high-resolution spectra that are difficult and time-consuming to obtain for large numbers of distant, dim galaxies.The training and validation data were generated to mimic Dark Energy Survey observations, and the analysis treats the simulator as noiseless.Source: §5 · pp. 15–16
5.1. Preliminary analysis and model comparison.
- For the OLS basis-function screen, p = 4 and m = 2 gave the highest NSE, 0.799, with 53 regression terms.The screen compared predictions on 500 validation points after fitting 2,000 training points.Source: §5.1 · pp. 16–17, Figure 5 · p. 17
- The largest NSE was for sparse correlation structure (8) with α = 1 and Legendre polynomials through degree 4, followed closely by the same model with α = 3/2; coverage was very close to nominal 95% with regression terms and sparse correlation.Source: §5.1 · pp. 17–18
- With Legendre regression terms, the posterior sum of correlation ranges moved away from the 0.185 cutoff; in the intercept-only model it remained near the boundary.The authors state that this cutoff is too small for the intercept-only GP to capture all data variability, but adequate for residual variability after adding the regression terms; the trace plots also show better mixing with regression terms.Source: §5.1 · p. 18, Figure 6 · p. 19
5.2. Full analysis.
- NSE was 0.831 and empirical coverage of 95% credible intervals was 94% for predictions on 80,000 validation points from 20,000 training points.The full analysis used 3,000 MCMC iterations, discarding 500 for burn-in.Source: §5.2 · p. 19
6. Discussion.
The discussion presents an anisotropic compactly supported GP with a regression mean and recommends selecting basis functions through a preliminary study. It also cautions that neither polynomial nor constant-mean structure is entirely satisfactory for extrapolation.
6. Discussion
- The proposed emulator combines an anisotropic compactly supported covariance with a regression mean, using a flat prior for regression parameters and a preliminary study to select basis functions.With the flat prior, all parameters except τ may be integrated out.Source: §6 · pp. 19–20
- Figure 7 plots photometric redshift predictions against spectroscopic redshift values, with transparency making densely over-plotted regions appear darker.Source: Figure 7 · p. 20
- Only a few Nash–Sutcliffe efficiency points were gained from shrinkage, and only for small data sets; the authors prefer preliminary-study basis selection for their target setting.They judge the added computational cost of shrinkage not worth the potential improvement.Source: §6 · p. 20
- Polynomial means can be problematic for extrapolation, while a constant-mean GP reverts to the global mean and ignores new inputs far outside the observed range.The authors say neither approach is entirely satisfactory and advise users to consider the mean structure when extrapolating.Source: §6 · p. 21
APPENDIX: POSTERIOR SAMPLING AND PREDICTION
The appendix samples the covariance-range parameters with an adaptive Metropolis algorithm, then uses conditional predictive means and variances to estimate the overall predictive distribution. These estimates support approximate pointwise credible intervals; global intervals require sampling the conditional multivariate distributions.
APPENDIX: POSTERIOR SAMPLING AND PREDICTION
- An adaptive Metropolis sampler rejects proposals outside the constrained parameter space and otherwise accepts them according to the integrated-likelihood ratio.Its proposal covariance adapts toward the posterior covariance and a target acceptance rate while retaining ergodicity under vanishing adaptation.Source: APPENDIX: POSTERIOR SAMPLING AND… · p. 21
- Conditional predictive means and variances at selected τ samples are combined by Monte Carlo approximation to estimate the overall predictive mean and variance.The number of τ samples can be chosen using an estimate of the smallest effective sample size among the τ elements.Source: APPENDIX: POSTERIOR SAMPLING AND… · pp. 21–22
- The resulting estimates support approximate pointwise credible intervals; global intervals instead require sampling from the conditional multivariate t distributions at each iteration.Source: APPENDIX: POSTERIOR SAMPLING AND… · p. 22