Gaussian-Process Emulators

A Gaussian process can interpolate expensive relativistic calculations while returning a covariance conditioned on its kernel, noise model, hyperparameters, and training data. That covariance is useful when calibrated against held-out high-fidelity results. It is not automatically the uncertainty of general relativity, the numerical simulation, or the astrophysical model. This page separates interpolation uncertainty from model discrepancy and treats confident error outside the training domain as a decisive failure.

Required background. Probability and Inference for Gravity supplies Gaussian conditioning, likelihoods, coverage, and calibration.

Helpful background. Formulations of General Relativity supplies the physical variables and constraints behind emulator outputs; Differential Equations for Relativity supplies high-fidelity solution and residual checks; Computational Relativity supplies versioned benchmarks, convergence, and reproducibility.

Conditional interpolation

For scalar training outputs $y=(y_1,\ldots,y_n)^T$ at parameters $X=(\lambda_1,\ldots,\lambda_n)$, assume

$$f(\lambda)\sim\operatorname{GP}(m(\lambda),k(\lambda,\lambda')), \qquad y=f(X)+\epsilon,$$

where $\epsilon$ is zero-mean Gaussian noise independent of $f$, with covariance $\Sigma_\epsilon$. Conditional on fixed kernel hyperparameters, $K_y=K(X,X)+\Sigma_\epsilon$. At a new point $\lambda_*$, the latent function $f_*$ has mean and variance

$$\mu_* = m_*+k_{*X}K_y^{-1}(y-m_X),$$ $$\sigma_*^2=k_{**}-k_{*X}K_y^{-1}k_{X*}.$$

The variance shown is for the latent function, not a new noisy output. For an independent new observation with variance $s_{\epsilon,*}^2$, add that variance to $\sigma_*^2$. Correlated numerical errors require the appropriate cross-covariances instead. Integrating over uncertain hyperparameters also adds uncertainty in the conditional mean; a single optimized kernel does not include it.

Implement these expressions through a factorization, never an explicit inverse. Record parameter scaling, mean, kernel family, amplitude and length scales, nugget or simulation-noise covariance, hyperpriors, optimization or sampling method, and conditioning. A tiny diagonal jitter added only to make a factorization succeed changes the probabilistic model and belongs in the record.

For a waveform, first reduce the time or frequency dependence and emulate coefficients $c_j(\lambda)$, amplitude and phase pieces, or values at empirical nodes. Joint coefficients may be correlated; fitting independent scalar processes is an approximation whose effect needs a held-out test. Doctor and collaborators lay out a Gaussian-process/reduced-order construction in Doctor et al. 2017, § II, Eqs. (3)–(10).

The posterior variance has a precise conditional meaning only if the kernel and likelihood describe the unknown function and the hyperparameter treatment is adequate. It usually shrinks near training points even when all simulations share the same physical bias. It can also be overconfident at a sharp feature that a smooth stationary kernel cannot express.

One observation, one prediction

Consider a zero-mean unit-variance GP with one noiseless observation $y_1=2$. If the covariance between the observed and requested point is $\rho$, conditioning gives

$$\mu_*=2\rho,\qquad \sigma_*^2=1-\rho^2.$$

For $\rho=0.8$, the prediction is $1.6$ with standard deviation $0.6$. A nominal pointwise 95% Gaussian interval is approximately $1.6\pm1.96(0.6)$. At the training point, $\rho=1$ gives exact interpolation and zero variance. Far away, a decaying kernel gives $\rho\to0$: the prediction returns to its prior mean and prior variance.

That last variance is the chosen prior scale, not a measured extrapolation error. An omitted physical effect larger than that scale can therefore make the far-away prediction confidently wrong.

Interpolation is not model discrepancy

Suppose observations or a more faithful calculation satisfy

$$y_{\rm ref}(\lambda) =f_{\rm sim}(\lambda)+\delta(\lambda)+\varepsilon,$$

where the emulator approximates $f_{\rm sim}$ and $\delta$ is simulator or model discrepancy. The GP interpolation variance for $f_{\rm sim}$ does not contain $\delta$ unless a discrepancy process is explicitly introduced and constrained by separate evidence. Calibration can then confound a physical parameter change with discrepancy: different combinations explain the same residual.

Kennedy and O’Hagan formulate probabilistic calibration with discrepancy in Kennedy and O’Hagan 2001, §§ 2–4. Brynjarsdóttir and O’Hagan show why weak or inappropriate discrepancy assumptions can bias physical calibration in Brynjarsdóttir and O’Hagan 2014, §§ 2–4. The lesson is not that discrepancy modeling is impossible; it is that interpolation covariance, numerical error, model-form error, and observational noise must be named separately.

For compact-binary waveforms, relevant omissions can include finite-resolution and extraction error, an approximate waveform family, missing modes, eccentricity or matter outside the model, and detector calibration. Adding these variances as if they were independent Gaussians is itself an assumption.

Comparison of methods and results

The evidence below was checked through 2026-08-18.

Position	Evidence in its tested regime	Boundary of the claim
established baseline	direct interpolation, reduced bases, and empirical interpolation with held-out high-fidelity comparisons	no predictive probability unless a calibrated stochastic model is added
useful controlled gain	Doctor et al. used 120 IMRPhenomD training waveforms over $q\in[1,3]$ and equal aligned spins $\chi_1=\chi_2\in[-0.5,0.5]$, demonstrating interpolation against that waveform model	the result is conditional on that waveform model, parameter box, preprocessing, kernel, and tests
adaptive promise	GP variance can guide acquisition toward under-resolved regions	it succeeds only when the variance tracks actual error; kernel misspecification can hide a feature
negative or limiting evidence	discrepancy and parameter calibration can be nonidentifiable and prior-sensitive	a narrow posterior can reflect assumptions rather than physical knowledge
not established	universal coverage under parameter shift, inclusion of simulator bias in ordinary GP variance, or reliable extrapolation	each requires new domain-specific evidence

The positive benchmark and uncertainty-guided sampling appear in Doctor et al. 2017, § III, Figs. 4–9, and § IV, Algorithm 1 and Figs. 10–12; kernel and cost limitations are discussed in Doctor et al. 2017, § V, Figs. 13–15. The inspected arXiv text has a numerical inconsistency: §III.C states a maximum mismatch $4.3\times10^{-5}$, while the Fig. 9 caption states $4.3\times10^{-3}$. We therefore do not use either number as an unambiguous acceptance target. The source supports useful interpolation inside the demonstrated domain. The competing position is that calibration and coverage can fail at boundaries, omitted regions, and model discrepancy, preventing a broader physical uncertainty claim.

Calibration benchmark

Freeze a versioned compact-binary waveform family, its reduced representation, parameter coordinates, alignment, norm, and software environment. Before training, partition it into:

• interpolation points scattered through the declared domain;
• complete faces and corners of $q\in[1,3]$, $\chi\in[-0.5,0.5]$;
• one contiguous withheld interior region;
• a narrow or rapidly varying feature region if the family contains one; and
• out-of-domain points used only to test warning or refusal.

For every held-out scalar coefficient, compute

$$z_i=\frac{y_i-\mu_i}{\sigma_i},$$

using the predictive variance appropriate to the reference (latent value or noisy observation). If $\sigma_i=0$, compare the absolute residual and report a violated deterministic prediction rather than silently dividing by zero. Report empirical 68% and 95% interval coverage, negative log predictive density, residual-versus-variance plots, and the worst standardized residual. For reconstructed waveforms, also report median, 95th-percentile, and maximum detector-weighted mismatch and the parameter location of each failure. Coverage is evaluated across genuinely held-out physical points, not time samples from a waveform used in training.

Compare at least a Matérn and squared-exponential kernel, reasonable hyperpriors or optimization restarts, and a simple distance-to-training baseline. Repeat the train/split procedure over declared random seeds. Increase training density, reduced-basis size, optimizer accuracy, and arithmetic precision separately. A finite acceptance rule might require empirical 95% coverage between 90% and 98% on the predeclared in-domain holdout and maximum mismatch below $10^{-4}$; both the thresholds and actual discrepancies must be reported. The contiguous-hole and boundary results carry more weight than an IID random split.

The independent check reconstructs held-out high-fidelity waveforms directly and recomputes overlaps with a separate quadrature. It also tests a deliberately adverse boundary or omitted-feature case. A low mean mismatch with undercoverage is a failed uncertainty claim.

The environment record pins data revisions and checksums, split manifest, preprocessing, basis, kernel and priors, optimizer, seeds, precision, linear-algebra and GP libraries, operating system or container, compiler or interpreter, accelerator, and inference code. Timings state hardware, warm-up, batching, and whether training cost is amortized.

Reproducing the calculation and recognizing failure

The principal failure modes are leakage through preprocessing, kernel smoothness that suppresses a sharp feature, ill-conditioned covariance, hyperparameter optimization trapped in a seed-dependent solution, uncalibrated variance, and extrapolation labeled interpolation. Another is adding emulator variance to physical uncertainties without tracing dependence or simulator discrepancy.

A clean reconstruction verifies inputs; rebuilds the reduced representation and fixed splits; trains every kernel and seed; sweeps training coverage, capacity, optimization tolerance, and precision; evaluates standardized residuals, coverage, likelihood scores, mismatch, conditioning, and cost; runs the adverse boundary and hole tests; and emits all thresholds and discrepancies in machine-readable form. Failure to cover or warn on the adverse case rejects the uncertainty claim even if mean in-domain error is small.

Review is triggered by a new simulation catalog or parameter domain, a changed kernel or discrepancy prior, an independent coverage failure, a software or arithmetic change that alters conditioning, or a reproducible benchmark that materially expands the validated region. Until then, the claim remains conditional on the dated domain and model.

Check your understanding

A GP trained on a biased simulator has posterior standard deviation $10^{-5}$ near a dense cluster of training points. A new high-fidelity calculation differs from the GP mean by $10^{-2}$ there. Does the small GP variance make the high-fidelity result implausible?

Solution

Not by itself. The variance measures uncertainty about interpolation of the simulator under the chosen GP model. A shared simulator discrepancy can remain large where training is dense. Treat the $10^{-2}$ residual as evidence against the combined interpolation-plus-discrepancy model, separate numerical and model-form contributions, and recalibrate coverage on held-out high-fidelity cases. Simply inflating the GP nugget would conceal rather than explain the failure.

Continue with Neural Surrogates for Relativistic Systems.

References

• Brynjarsdóttir, Jenný, and Anthony O’Hagan. “Learning about Physical Parameters: The Importance of Model Discrepancy.” Inverse Problems 30, no. 11 (2014): 114007. DOI.
• Doctor, Zoheyr, Ben Farr, Daniel E. Holz, and Michael Pürrer. “Statistical Gravitational Waveform Models: What to Simulate Next?” Physical Review D 96 (2017): 123011. arXiv:1706.05408. DOI.
• Kennedy, Marc C., and Anthony O’Hagan. “Bayesian Calibration of Computer Models.” Journal of the Royal Statistical Society: Series B 63, no. 3 (2001): 425–464. DOI.