Society for Industrial and Applied Mathematics: SIAM/ASA Journal on Uncertainty Quantification: Table of Contents

The Bayesian Finite Element Method in Inverse Problems: A Critical Comparison between Probabilistic Models for Discretization Error

Anne Poot — 2026-04-01T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 287-312, June 2026.
Abstract.When using the finite element method (FEM) in inverse problems, its discretization error can produce parameter estimates that are inaccurate and overconfident. The Bayesian finite element method (BFEM) provides a probabilistic model for the epistemic uncertainty due to discretization error. In this work, we apply BFEM to various inverse problems and compare its performance to the random mesh finite element method (RM-FEM) and the statistical finite element method (statFEM), which serve as a frequentist and inference-based counterpart to BFEM. We find that by propagating this uncertainty to the posterior, BFEM can produce more accurate parameter estimates and prevent overconfidence compared to FEM. Because the BFEM covariance operator is designed to leave uncertainty only in the appropriate space, orthogonal to the FEM basis, BFEM is able to outperform RM-FEM, which does not have such a structure to its covariance. Although it is also possible to use a model misspecification formulation such as statFEM to infer the discretization error downstream rather than model it at the source, the feasibility of such an approach is contingent on the availability of sufficient data. We find that the BFEM is the most robust way to consistently propagate uncertainty due to discretization error to the posterior of a Bayesian inverse problem.

Dirichlet–Neumann Averaging: The DNA of Efficient Gaussian Process Simulation

Robert Kutri — 2026-04-02T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 313-340, June 2026.
Abstract.Gaussian processes (GPs) and Gaussian random fields (GRFs) are essential for modeling spatially varying stochastic phenomena. Yet the efficient generation of corresponding realizations on high-resolution grids remains challenging, particularly when a large number of realizations are required. This paper presents two novel contributions. First, we propose a new methodology based on Dirichlet–Neumann averaging (DNA) to generate GPs and GRFs with isotropic covariance on regularly spaced grids. The combination of discrete cosine and sine transforms in the DNA sampling approach allows for rapid evaluations without the need for modification or padding of the desired covariance function. While this introduces an error in the covariance, our numerical experiments show that this error is negligible for most relevant applications, representing a trade-off between efficiency and precision. We provide explicit error estimates for Matérn covariances. The second contribution links our new methodology to the stochastic partial differential equation (SPDE) approach for sampling GRFs. We demonstrate that the concepts developed in our methodology can also guide the selection of boundary conditions in the SPDE framework. We prove that averaging specific GRFs sampled via the SPDE approach yields genuinely isotropic realizations without domain extension, with the error bounds established in the first part remaining valid.

Local Sensitivity Analysis for Bayesian Inverse Problems

Jürgen Dölz — 2026-04-02T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 341-368, June 2026.
Abstract.We present an extension of local sensitivity analysis, also referred to as the perturbation approach for uncertainty quantification, to Bayesian inverse problems. More precisely, we show how moments of random variables with respect to the posterior distribution can be approximated efficiently by asymptotic expansions. This is under the assumption that the measurement operators and prediction functions are sufficiently smooth and that their corresponding stochastic moments with respect to the prior distribution exist. Numerical experiments are presented to the illustrate the theoretical results.

Bayesian Inference for Non-synchronously Observed Diffusions

Ajay Jasra — 2026-04-07T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 369-393, June 2026.
Abstract.We consider the problem of Bayesian inference for bi-variate data observed in time but with observation times which occur non-synchronously. In particular, this occurs in a wide variety of applications in finance, such as high-frequency trading or crude oil futures trading. We adopt a diffusion model for the data and formulate a Bayesian model with priors on unknown parameters along with a latent representation for the so-called missing data. We then consider computational methodology to fit the model using Markov chain Monte Carlo (MCMC). We have to resort to time-discretization methods, as the complete data likelihood is intractable and this can cause considerable issues for MCMC when the data are observed in low frequencies. In the context of high frequency observations, we present a simple particle MCMC method based on an Euler–Maruyama time discretization, which can be enhanced using multilevel Monte Carlo (MLMC). In the low frequency observation regime, we introduce a novel bridging representation of the posterior in continuous time to deal with the issues of MCMC in this case. This representation is discretized and fitted using MCMC and MLMC. We apply our methodology to real and simulated data to establish the efficacy of our methodology.

Tensor-Variate Gaussian Process Regression for Efficient Emulation of Complex Systems: Comparing Regressor and Covariance Structures in Outer Product and Parallel Partial Emulators

D. Semochkina — 2026-04-20T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 394-414, June 2026.
Abstract.Multioutput Gaussian process regression has become an important tool in uncertainty quantification, for building emulators of computationally expensive simulators, and other areas such as multitask machine learning. We present a holistic development of tensor-variate Gaussian process (TvGP) regression, appropriate for arbitrary dimensional outputs where a Kronecker product structure is appropriate for the covariance. We show how two common approaches to problems with two-dimensional output, outer product emulators (OPE) and parallel partial emulators, are special cases of TvGP regression and hence can be extended to higher output dimensions. Focusing on the important special case of matrix output, we investigate the relative performance of these two approaches. The key distinction is the additional dependence structure assumed by the OPE, and we demonstrate when this is advantageous through two case studies, including application to a spatial-temporal influenza simulator.

Gradient-Free Sequential Bayesian Experimental Design via Interacting Particle Systems

Robert Gruhlke — 2026-04-23T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 415-455, June 2026.
Abstract.We introduce a gradient-free framework for Bayesian Optimal Experimental Design (BOED) in sequential settings, aimed at complex systems where gradient information is unavailable. Our method combines Ensemble Kalman Inversion (EKI) for design optimization with the Affine-Invariant Interacting Langevin Dynamics (ALDI) sampler for efficient posterior sampling—both of which are derivative-free and ensemble-based. To address the computational challenges posed by nested expectations in BOED, we propose variational Gaussian and parametrized Laplace approximations that provide tractable upper and lower bounds on the Expected Information Gain (EIG). These approximations enable scalable utility estimation in high-dimensional spaces and PDE-constrained inverse problems. We demonstrate the performance of our framework through numerical experiments ranging from linear Gaussian models to PDE-based inference tasks, highlighting the method’s robustness, accuracy, and efficiency in information-driven experimental design.

Cross-Validation Based Adaptive Sampling for Multilevel Gaussian Process Models

Louise M. Kimpton — 2026-04-27T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 456-481, June 2026.
Abstract.Complex computer codes or models can often be run in a hierarchy of different levels of complexity ranging from the very basic to the sophisticated. The top levels in this hierarchy are typically expensive to run, which limits the number of possible runs. To make use of simulations over all levels, and crucially improve predictions at the top level, we use multilevel Gaussian process emulators (GPs). The accuracy of the GP greatly depends on the design of the training points. In this paper, we present a multilevel adaptive sampling algorithm to sequentially increase the set of design points to optimally improve the fit of the GP. The normalized squared expected leave-one-out cross-validation error (ES-LOO) is calculated at all unobserved locations, and a new design point is chosen using expected improvement combined with a repulsion function to find the maximum ES-LOO. We use ES-LOO to obtain a model-free measure of prediction error at each level. This criterion is calculated for each model level weighted by an associated cost for the code at that level. Hence, at each iteration, our algorithm optimises for both the new point location and the model level. The algorithm is extended to batch selection as well as single point selection, where batches can be designed for single levels or optimally across all levels. We apply the new multilevel ES-LOO algorithm to a range of examples, and make comparisons to existing design methods. This includes comparisons with single-shot methods as well as more recent sequential design algorithms. A toolbox containing relevant code can be found at: github.com/EXA-UQ/EXAUQ-Toolbox.

High-Dimensional Stochastic Finite Volumes Using the Tensor Train Format

Juliette Dubois — 2026-04-29T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 482-507, June 2026.
Abstract.A method for the uncertainty quantification of nonlinear hyperbolic conservation laws with many uncertain parameters is presented. The method combines stochastic finite volume methods and tensor trains in a novel way: the dimensions of physical space and time are kept as full tensors, while all stochastic dimensions are compressed together into a tensor train. The resulting hybrid format has one tensor train for each spatial cell and each time step. The MUSCL scheme is adapted to the proposed hybrid format, and its feasibility is demonstrated through several test cases. For the scalar Burgers’ equation, we conduct a convergence study and compare the results with those obtained using the full tensor train (full-TT) format with three stochastic parameters. The equation is then solved for an increasing number of stochastic dimensions. For systems of conservation laws, we focus on the Euler equations. A parameter study and a comparison with the full-TT format are carried out for the Sod shock tube problem. As a more complex application, we investigate the Shu–Osher problem, which involves intricate wave interactions. The presented method opens new avenues for integrating uncertainty quantification with established numerical schemes for hyperbolic conservation laws.

A Kernel-Based Approach for Gaussian Process Modeling with Functional Information

D. Andrew Brown — 2026-05-05T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 508-533, June 2026.
Abstract.Gaussian processes (GPs) are ubiquitous tools for modeling and predicting continuous processes in physical and engineering sciences. This is partly due to the fact that one may employ a GP as an interpolator while facilitating straightforward uncertainty quantification at other locations. In addition to training data, it is sometimes the case that available information is not in the form of a finite collection of points. For example, boundary value problems contain information on the boundary of a domain, or underlying physics lead to known behavior on an entire uncountable subset of the domain of interest. While an approximation to such known information may be obtained via pseudo-training points in the known subset, such a procedure is ad hoc with little guidance on the number of points to use, nor the behavior as the number of pseudo-observations grows large. We propose and construct GPs that unify, via reproducing kernel Hilbert space, the typical finite training data case with the case of having uncountable information by exploiting the equivalence of conditional expectation and orthogonal projections in Hilbert space. We show existence of the proposed process and establish that it is the limit of a conventional GP conditioned on an increasing number of training points. We illustrate the flexibility and advantages of our proposed approach via numerical experiments.

Neural Field Equations with Random Data

Daniele Avitabile — 2026-05-13T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 534-567, June 2026.
Abstract.We study neural field equations, which are prototypical models of large-scale cortical activity, subject to random data. We view these spatially extended, nonlocal evolution equations as a Cauchy problem on abstract Banach spaces, with randomness in the synaptic kernel, firing rate function, external stimuli, and initial conditions. We determine conditions on the random data that guarantee existence, uniqueness, and measurability of the solution for uncertainty quantification (UQ) and examine the regularity of the solution in relation to the regularity of the inputs. We present results for linear and nonlinear neural fields and for the two most common functional setups in the numerical analysis of this problem. In addition to the continuous problem, we analyze in abstract form neural fields that have been spatially discretized, setting the foundations for analyzing UQ schemes.

Deterministic Kalman Filters for Dynamical Systems with Parametric Uncertainty

Karl Kunisch — 2026-05-27T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 568-596, June 2026.
Abstract.The Kalman(-Bucy) filter is the natural choice for the state reconstruction of disturbed, linear dynamical systems based on flawed and incomplete measurements. Taking a deterministic viewpoint, this work investigates possible extensions of the concept to systems with uncertain dynamics and noise covariances. In a theoretical analysis, error bounds in terms of the variance of the uncertainties are derived. The article concludes with a numerical implementation of two example systems, allowing for a comparison of the estimators.

Hierarchical Gaussian Random Field Sampling for Multilevel Markov Chain Monte Carlo: Coupling Stochastic Partial Differential Equation and the Karhunen–Loève Decomposition

Sohail Reddy — 2026-05-27T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 2, Page 597-620, June 2026.
Abstract.This work introduces structure preserving hierarchical decompositions for sampling Gaussian random fields (GRFs) within the context of multilevel Bayesian inference in high-dimensional space. Existing scalable hierarchical sampling methods, such as those based on stochastic partial differential equations (SPDEs), often reduce the dimensionality of the sample space at the cost of accuracy of inference. Other approaches, such that those based on Karhunen-Loève (KL) expansions, offer sample space dimensionality reduction but sacrifice GRF representation accuracy and ergodicity of the Markov chain Monte Carlo (MCMC) sampler and are computationally expensive for high-dimensional problems. The proposed method integrates the dimensionality reduction capabilities of KL expansions with the scalability of SPDE-based sampling, thereby providing a robust, unified framework for high-dimensional uncertainty quantification (UQ) that is scalable and accurate, preserves ergodicity, and offers dimensionality reduction of the sample space. The hierarchy in our multilevel algorithm is derived from the geometric multigrid hierarchy. By constructing a hierarchical decomposition that maintains the covariance structure across the levels in the hierarchy, the approach enables efficient coarse-to-fine sampling while ensuring that all samples are drawn from the desired distribution. The effectiveness of the proposed method is demonstrated on a benchmark subsurface flow problem, demonstrating its effectiveness in improving computational efficiency and statistical accuracy. Our proposed technique is more efficient and accurate and displays better convergence properties than existing methods for high-dimensional Bayesian inference problems.

Active Learning via Heteroskedastic Rational Kriging

Shangkun Wang — 2026-01-02T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 1-18, March 2026.
Abstract.Active learning methods for emulating complex computer models that rely on stationary Gaussian processes tend to produce design points that uniformly fill the entire experimental region, which can be wasteful for functions which vary only in small regions. In this article, we propose a new Gaussian process model that captures the heteroskedasticity of the function. This is achieved by modifying the recently proposed rational kriging model, which has enough flexibility to seamlessly incorporate heteroskedasticity. Active learning using this new model can place design points in the more interesting regions of the response surface, and thus obtain surrogate models with better accuracy. The proposed active learning method is compared with the state-of-the-art methods using simulations and two real datasets. It is found to have comparable or better performance relative to other nonstationary Gaussian process-based methods, but to be faster by orders of magnitude.

Goal-Oriented Bayesian Optimal Experimental Design for Nonlinear Models Using Markov Chain Monte Carlo

Shijie Zhong — 2026-01-02T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 19-47, March 2026.
Abstract.Optimal experimental design (OED) provides a systematic approach to quantify and maximize the value of experimental data. Under a Bayesian approach, conventional OED maximizes the expected information gain (EIG) on model parameters. However, we are often interested not in the parameters themselves but in predictive quantities of interest (QoIs) that depend on the parameters in a nonlinear manner. We present a computational framework of predictive goal-oriented OED (GO-OED) suitable for nonlinear observation and prediction models that seeks the experimental design providing the greatest EIG on the QoIs. In particular, we propose a nested Monte Carlo estimator for the QoI EIG, featuring Markov chain Monte Carlo for posterior sampling and kernel density estimation for evaluating the posterior-predictive density and its Kullback–Leibler divergence from the prior-predictive density. The GO-OED design is then found by maximizing the EIG over the design space using Bayesian optimization. We demonstrate the effectiveness of the overall nonlinear GO-OED method and illustrate its difference versus the conventional non-GO-OED through various test problems and an application of sensor placement for source inversion in a convection-diffusion field.

Random Fourier Features Based Gaussian Process Models for Stochastic Simulations

Ying Wu — 2026-01-06T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 48-76, March 2026.
Abstract.Stochastic simulations are increasingly used to describe complex systems with uncertainties. To better characterize the uncertainty of such a simulation, we propose a random Fourier features method to fit full Bayesian Gaussian process emulators for these simulations. The random Fourier features technique uses low-dimensional features of the correlation function of the Gaussian process to achieve a low-rank approximation of the correlation matrix. We prove the convergence of the proposed random Fourier features method. Simulation results show that the proposed method can significantly reduce computation time while maintaining high prediction accuracy. The advantages of the proposed method are also illustrated using a modern subway simulation.

Domain Uncertainty Quantification for the Lippmann–Schwinger Volume Integral Equation

Fernando Henríquez — 2026-01-08T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 77-110, March 2026.
Abstract.In this work, we consider the propagation of acoustic waves in unbounded domains characterized by a constant wavenumber, except possibly within a bounded region. The geometry of this inhomogeneity is assumed to be uncertain, and we are particularly interested in studying how this uncertainty propagates throughout the physical model under consideration. A key step in our analysis involves recasting the physical model—originally defined in an unbounded domain—into a computationally tractable formulation based on volume integral equations, specifically the Lippmann–Schwinger equation. We show that both the leading operator in this volume integral formulation and its solution depend holomorphically on shape variations of the support of the inhomogeneity. This property, known as shape holomorphy, is crucial for the analysis and implementation of various methods used in computational uncertainty quantification (UQ). We explore the implications of this result for both forward and inverse UQ and provide numerical experiments that illustrate and confirm the theoretical predictions.

Subspace Splitting Fast Sampling from Gaussian Posterior Distributions of Linear Inverse Problems

Daniela Calvetti — 2026-02-03T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 111-141, March 2026.
Abstract.It is well-known that the posterior density of linear inverse problems with Gaussian prior and Gaussian likelihood is also Gaussian, hence completely described by its covariance and expectation. Sampling from a Gaussian posterior may be important in the analysis of various non-Gaussian inverse problems in which estimates from a Gaussian posterior distribution constitute an intermediate stage in a Bayesian workflow. Sampling from a Gaussian distribution is straightforward if the Cholesky factorization of the covariance matrix or its inverse is available; however, when the unknown is high dimensional, the computation of the posterior covariance may be unfeasible. If the linear inverse problem is underdetermined, it is possible to exploit the orthogonality of the fundamental subspaces associated with the coefficient matrix together with the idea behind the randomize-then-optimize approach to design a low complexity posterior sampler that does not require the posterior covariance to be formed. The performance of the proposed sampler is illustrated with a few computed examples, including non-Gaussian problems with nonlinear forward model, and hierarchical models comprising a conditionally Gaussian submodel.

Iterative Methods for Full-Scale Gaussian Process Approximations for Large Spatial Data

Tim Gyger — 2026-02-23T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 142-167, March 2026.
Abstract.Gaussian processes are flexible probabilistic regression models which are widely used in statistics and machine learning. However, a drawback is their limited scalability to large data sets. To alleviate this, full-scale approximations (FSAs) combine predictive process methods and covariance tapering, thus approximating both global and local structures. We show how iterative methods can be used to reduce computational costs in calculating likelihoods, gradients, and predictive distributions with FSAs. In particular, we introduce a novel preconditioner and show theoretically and empirically that it accelerates the conjugate gradient method’s convergence speed and mitigates its sensitivity with respect to the FSA parameters and the eigenvalue structure of the original covariance matrix, and we demonstrate empirically that it outperforms a state-of-the-art pivoted Cholesky preconditioner. Furthermore, we introduce an accurate and fast way to calculate predictive variances using stochastic simulation and iterative methods. In addition, we show how our newly proposed fully independent training conditional (FITC) preconditioner can also be used in iterative methods for Vecchia approximations. In our experiments, it outperforms existing state-of-the-art preconditioners for Vecchia approximations. All methods are implemented in a free C++ software library with high-level Python and R packages (https://github.com/fabsig/GPBoost).

PCE-Net: High-Dimensional Surrogate Modeling for Learning Uncertainty

Paz Fink Shustin — 2026-02-26T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 168-196, March 2026.
Abstract.Learning data representations under uncertainty is an important task that emerges in numerous scientific computing and data analysis applications. However, uncertainty quantification techniques are computationally intensive and become prohibitively expensive for high-dimensional data. In this study, we introduce a dimensionality reduction surrogate modeling approach for representation learning and uncertainty quantification that aims to deal with data of moderate to high dimensions. The approach involves a two-stage learning process: (1) employing a variational autoencoder to learn a low-dimensional representation of the input data distribution, and (2) harnessing a polynomial chaos expansion formulation to map the low-dimensional distribution to the output target. The model enables us, without any prior statistical assumptions on the data, to (a) capture the system dynamics efficiently in the low-dimensional latent space, (b) learn under uncertainty, a representation of the data and a mapping between input and output distributions, (c) estimate this uncertainty in the high-dimensional data system, and (d) match high-order moments of the output distribution. Numerical results are presented to illustrate the performance of the proposed method.

Batch Bayesian Optimization via Particle Gradient Flows

Enrico Crovini — 2026-03-05T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 197-220, March 2026.
Abstract.Bayesian optimization (BO) methods seek to find global optima of objective functions which are only available as a black-box or are expensive to evaluate. Such methods construct a surrogate model for the objective function, quantifying the uncertainty in that surrogate through Bayesian inference. Objective evaluations are sequentially determined by maximizing an acquisition function at each step. However, this ancilliary optimization problem can be highly nontrivial to solve, due to the nonconcavity of the acquisition function, particularly in the case of batch Bayesian optimization, where multiple points are selected in every step. In this work we reformulate batch BO as an optimization problem over the space of probability measures. We construct a new acquisition function based on multipoint expected improvement, which is concave over the space of probability measures. Practical schemes for solving this “inner" optimization problem arise naturally as gradient flows of this objective function. We demonstrate the efficacy of this new method on different benchmark functions and compare with state-of-the-art batch BO methods.

Neural Network Approaches for Variance Reduction in Fluctuation Formulas

G. A. Pavliotis — 2026-03-05T08:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 221-255, March 2026.
Abstract.We propose a method utilizing physics-informed neural networks (PINNs) to solve Poisson equations that serve as control variates in the computation of transport coefficients via fluctuation formulas, such as the Green–Kubo and generalized Einstein-like formulas. By leveraging approximate solutions to the Poisson equation constructed through neural networks, our approach significantly reduces the variance of the estimator at hand. We provide an extensive numerical analysis of the estimators and detail a methodology for training neural networks to solve these Poisson equations. The approximate solutions are then incorporated into Monte Carlo simulations as effective control variates, demonstrating the suitability of the method for moderately high-dimensional problems where fully deterministic solutions are computationally infeasible.

Local Transfer Learning Gaussian Process Modeling, with Applications to Surrogate Modeling of Expensive Computer Simulators

Xinming Wang — 2026-03-10T07:00:00Z

SIAM/ASA Journal on Uncertainty Quantification, Volume 14, Issue 1, Page 256-286, March 2026.
Abstract.A critical bottleneck for scientific progress is the costly nature of computer simulations for complex systems. Surrogate models provide an appealing solution: such models are trained on simulator evaluations, then used to emulate and quantify uncertainty on the expensive simulator at unexplored inputs. In many applications, one often has available data on related systems. For example, in designing a new jet turbine, there may be existing studies on turbines with similar configurations. A key question is how information from such “source” systems can be transferred for effective surrogate training on the “target” system of interest. We thus propose a new LOcal transfer Learning Gaussian Process (LOL-GP) model, which leverages a carefully designed Gaussian process to transfer such information for surrogate modeling. The key novelty of the LOL-GP is a latent regularization model, which identifies regions where transfer should be performed and regions where it should be avoided. Such a “local transfer” property is present in many scientific systems: for certain parameters, systems may behave similarly and thus transfer is beneficial; for other parameters, they may behave differently and thus transfer is detrimental. By accounting for local transfer, the LOL-GP can temper the risk of “negative transfer,” i.e., the risk of worsening predictive performance from information transfer. We derive a Gibbs sampling algorithm for efficient posterior predictive sampling on the LOL-GP for both the multisource and multifidelity transfer settings. We then show, via a suite of numerical experiments and an application for jet turbine design, the improved surrogate performance of the LOL-GP over existing methods.