Bayesian Neural Network Surrogates

• Bayesian neural network surrogates are probabilistic emulators that represent uncertainty in expensive mappings by approximating posterior predictive distributions.
• They employ diverse inference methods, including variational, sampling-based, and ensemble techniques, to tackle intractable posteriors and enhance model accuracy.
• Their architectures adapt to various applications, from physics-constrained PDE solvers to gradient-informed Bayesian optimization, facilitating robust decision-making.

Bayesian neural network surrogates are neural-network emulators in which the surrogate mapping is endowed with a posterior over parameters, or with an explicitly approximated posterior predictive, so that prediction at a new input $x^*$ is represented by $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$. In the recent literature, this formulation is used to emulate expensive black-box objectives in Bayesian optimization, high-fidelity physical simulators, PDE solution operators, and structural-analysis pipelines, while exposing epistemic uncertainty for acquisition, exploration, calibration, and risk-aware triage (Fotias et al., 29 Jul 2025, Park et al., 30 Jan 2026, Zhang et al., 2021).

1. Probabilistic formulation and surrogate semantics

A standard finite-width BNN surrogate begins with a regression model $y_i=f(x_i;w)+\varepsilon_i$, $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$, together with a Gaussian prior over parameters, typically $p(w)=\mathcal N(w;0,\sigma_w^2 I)$. The posterior $p(w\mid D)\propto p(D\mid w)p(w)$ is intractable for nonlinear networks, and the surrogate is therefore defined operationally through an approximate posterior or a direct approximation to the posterior predictive (Fotias et al., 29 Jul 2025).

The likelihood need not be restricted to ordinary supervised residuals. In gradient-informed BO surrogates, noisy gradients are modeled through $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$, $\eta_i\sim\mathcal N(0,\sigma_{n,\nabla}^2 I_D)$, and the training objective augments the usual expected log-likelihood with a weighted gradient term $\lambda_\nabla$ (Makrygiorgos et al., 14 Apr 2025). In physics-constrained settings, the likelihood can instead be built from discretized PDE residuals or from residual vectors induced by weak-form discretization, so the surrogate is trained without labels in the usual sense (Zhang et al., 2021). A more general formulation replaces the ELBO-centric view entirely and trains an implicit posterior model by directly maximizing a Monte Carlo approximation to the posterior predictive; in that setting, “only a likelihood needs to be specified,” including loss functions and synthetic likelihoods (Dabrowski et al., 2022).

This breadth of likelihood construction is central to the notion of a BNN surrogate. The surrogate is not merely a stochastic regressor: it is a probabilistic emulator whose training target may be an objective function, a simulator response, a PDE residual, or a multi-modal observation model. This suggests that the defining feature of the class is posterior predictive uncertainty over expensive mappings rather than any single inference algorithm or architecture.

2. Posterior approximation mechanisms

The dominant inference mechanisms in the cited literature are variational, sampling-based, and ensemble-based. Mean-field variational inference uses a fully factorized Gaussian $q(w)$ and optimizes an ELBO of the form

$$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$0

with stochastic optimization and reparameterization; this appears in BO surrogates, bridge-assessment surrogates, and gradient-informed BNNs (Fotias et al., 29 Jul 2025, Kuhn et al., 29 Sep 2025, Makrygiorgos et al., 14 Apr 2025). TensorFlow Probability Flipout layers provide another variational route in PDE surrogates, where pseudo-independent perturbations reduce gradient variance while preserving an ELBO interpretation (Zhang et al., 2021).

Sampling-based methods remain the most faithful but most expensive approximations. In BO benchmarks, the literature considers HMC, NUTS, SGLD, and SGHMC; the 2023 comparative BO study concludes that HMC is the most successful approximate inference procedure for fully stochastic BNNs, while SGHMC yields less faithful uncertainty (Li et al., 2023). In high-dimensional stochastic PDE and inverse-problem settings, HMC is used directly on BNN parameters, with predictive moments obtained by Monte Carlo over posterior draws $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$1 (Jung et al., 2022).

A distinct line of work replaces explicit density approximation with an implicit posterior sampler $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$2 or $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$3, where $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$4 or $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$5 and $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$6 is drawn from a simple base distribution. The training objective is a Monte Carlo estimate of the posterior predictive, and the conditional version can produce input-dependent parameter sets $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$7, thereby increasing the apparent complexity of the primary model (Dabrowski et al., 2022).

Anchored ensembling provides a separate bridge between standard deep-learning practice and Bayesian posterior sampling. Pearce et al. define an “anchored” MAP problem in which each ensemble member is regularized around its own prior draw $$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$8, using a loss

$$p(y^*\mid x^*,D)=\int p(y^*\mid x^*,w)\,p(w\mid D)\,dw$$9

Their theorem states that for wide neural networks with Gaussian prior and Gaussian likelihood, as width $y_i=f(x_i;w)+\varepsilon_i$0 and $y_i=f(x_i;w)+\varepsilon_i$1 is scaled as $y_i=f(x_i;w)+\varepsilon_i$2, each solution to the anchored MAP problem is distributed as an exact posterior draw, so averaging over independent solves converges to the true Bayesian posterior predictive (Pearce et al., 2018). This result is notable because it reinterprets a scalable ensemble construction as an approximate BNN surrogate with principled posterior semantics.

3. Architectural patterns for surrogate construction

The architectural design space is heterogeneous because surrogate tasks differ sharply in fidelity structure, dimensionality, modality, and physical constraints.

B-INN replaces generic deep architectures with an interpolation-plus-tensor-decomposition construction. In one dimension, it represents $y_i=f(x_i;w)+\varepsilon_i$3, while in $y_i=f(x_i;w)+\varepsilon_i$4 dimensions it adopts a rank-$y_i=f(x_i;w)+\varepsilon_i$5 CP decomposition,

$y_i=f(x_i;w)+\varepsilon_i$6

with each $y_i=f(x_i;w)+\varepsilon_i$7 expanded in fixed basis functions. Bayesian inference proceeds blockwise through alternating-direction Bayesian linear regressions with closed-form Gaussian updates, and the function-space limit is a GP after suitable $y_i=f(x_i;w)+\varepsilon_i$8 normalization (Park et al., 30 Jan 2026).

Multi-fidelity GPBNN separates low- and high-fidelity uncertainty handling. The low-fidelity code is modeled by a GP, and the high-fidelity code is modeled by a one-hidden-layer BNN whose input augments $y_i=f(x_i;w)+\varepsilon_i$9 with low-fidelity emulator information. Three transfer schemes are described: Mean-Std, Quantiles, and a Gauss–Hermite method in which quadrature nodes $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$0 are propagated through the BNN. Posterior sampling for the high-fidelity BNN is performed with NUTS (Kerleguer et al., 2023).

Multi-modal surrogate architectures exploit auxiliary observations. One approach is a joint model that emits a concatenated mean vector for the primary quantity of interest and the auxiliary modalities; another is a layered model that first predicts each auxiliary modality and then feeds their posterior means, together with $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$1, into a second BNN for the main output. In both cases, the cited work uses a Matrix-Normal–Inverse-Wishart conjugate prior on the last layer and performs SVI only for the non-conjugate parameters, while analytically recomputing the conditional last-layer posterior at each step (Taylor et al., 26 Sep 2025).

Physics-constrained surrogates use architectures adapted to PDE structure. One study employs encoder–decoder CNNs with deterministic or probabilistic convolutional layers, taking Dirichlet and Neumann boundary conditions as image-like channels and defining the loss through a vectorized discretized residual of the weak form (Zhang et al., 2021). Another employs Fourier-feature networks in which $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$2 is followed by a small fully connected network, with the likelihood defined by PDE residuals and sensor data (Jung et al., 2022).

A further pattern injects prior information in function space rather than directly in weight space. In a mechanics surrogate, realizations from a functional prior—such as a low-fidelity model or GP—are first fit by small networks, producing a set of parameter vectors whose empirical covariance is factorized by thin SVD. The resulting low-rank correlated Gaussian prior is then used in an anchored-ensemble objective, allowing the weight-space prior to inherit mean trends and uncertainty profiles from function space (Ghorbanian et al., 2024).

4. Downstream workflows and operational uses

A major use of BNN surrogates is Bayesian optimization. In the CCS study, the GP surrogate is replaced by a BNN and refit after each BO iteration; single-objective acquisition uses MC-EI, $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$3, while the multi-objective setting uses EHVI approximated by averaging hypervolume increments over posterior samples (Fotias et al., 29 Jul 2025). The broader BO comparison likewise evaluates finite-width BNNs, infinite-width BNNs, linearized Laplace approximations, deep kernel learning, and deep ensembles as surrogate families, with acquisition functions including EI and GP-UCB and optimization via L-BFGS, Sobol candidate selection, or gradient-free search depending on differentiability and domain structure (Li et al., 2023).

Gradient information can be folded into the surrogate itself rather than only the acquisition. The gradient-informed BNN augments the ELBO with a gradient log-likelihood term and then reuses standard BO acquisition functions such as LCB and LogEI without further modification, because the predictive uncertainty already reflects gradient learning (Makrygiorgos et al., 14 Apr 2025).

Active learning is another recurrent use case. Anchored ensembles are explicitly proposed as a basis for BALD or mutual information calculations, and the same ensemble posterior predictive can be used to draw function samples for Bayesian optimization or to drive RL exploration through Q-network posterior variance (Pearce et al., 2018). B-INN turns this into a large-scale industrial active-learning loop by evaluating the average predictive standard deviation over spatial points, selecting the parameter $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$4 maximizing $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$5, and retraining by warm-start on the posterior means (Park et al., 30 Jan 2026).

Engineering surrogates often culminate in decision rules rather than just predictions. In bridge pre-assessment, independent BNN heads predict code-compliance factors and calibrated epistemic intervals, after which a triage rule labels a structure “red” if $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$6, “orange” if $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$7 but $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$8, and “green” if $\varepsilon_i\sim\mathcal N(0,\sigma_n^2)$9 (Kuhn et al., 29 Sep 2025). In multi-modal surrogates, the stated outer-loop applications include Bayesian optimization, inverse problems, and sensitivity analysis (Taylor et al., 26 Sep 2025). In PDE and SPDE settings, posterior samples are used directly for uncertainty quantification in forward and inverse problems, including posterior means, standard deviations, and covariance estimates over fields (Jung et al., 2022, Zhang et al., 2021).

5. Empirical behavior and comparative findings

Empirical results in BO are markedly problem dependent. The 2023 comparative study reports five central findings: the ranking of methods is highly problem dependent; HMC is the most successful approximate inference procedure for fully stochastic BNNs; full stochasticity may be unnecessary because deep kernel learning is relatively competitive; deep ensembles perform relatively poorly; and infinite-width BNNs are particularly promising, especially in high dimensions (Li et al., 2023). The same study further reports that i-NNGP consistently outperforms alternatives as dimensionality grows in its high-dimensional tests.

The CCS BO study presents a more guarded picture for finite-width BNN surrogates. In a single-objective sequestration-rate problem, SVI and MCMC converged fastest to the highest objective values, MC-dropout was slowest, GP plateaued after iteration 3, and the infinite-width BNN improved early and then stalled. In a multi-objective hypervolume setting, SVI achieved the largest hypervolume in one case study, but in a second case study with approximately 961 variables the GP outperformed all BNN surrogates in aggregate Pareto hypervolume and converged fastest; the authors therefore state that GPs often matched or outperformed BNNs even in high dimensions for these CCS tasks (Fotias et al., 29 Jul 2025).

B-INN is evaluated primarily against BNNs and GPs in large-scale physical-system settings. The paper states that B-INNs can be from 20 times to 10,000 times faster with a robust uncertainty estimation compared to Bayesian neural networks and Gaussian processes, and proves linear $p(w)=\mathcal N(w;0,\sigma_w^2 I)$0 inference complexity with respect to the number of training samples (Park et al., 30 Jan 2026). In one-dimensional synthetic regression, GP fits fail beyond $p(w)=\mathcal N(w;0,\sigma_w^2 I)$1, whereas B-INN with $p(w)=\mathcal N(w;0,\sigma_w^2 I)$2 trains in $p(w)=\mathcal N(w;0,\sigma_w^2 I)$3 on $p(w)=\mathcal N(w;0,\sigma_w^2 I)$4 points and achieves test RMSE matching GP within $p(w)=\mathcal N(w;0,\sigma_w^2 I)$5. In active learning on parametric PDEs, final RMSE is reported as approximately $p(w)=\mathcal N(w;0,\sigma_w^2 I)$6 for Poisson and $p(w)=\mathcal N(w;0,\sigma_w^2 I)$7 for heat conduction, with training times substantially below BNN-HMC and BNN-VI (Park et al., 30 Jan 2026).

Anchored ensembles are reported to match or outperform Deep Ensembles on UCI regression in high-epistemic tasks and to recover GP-like credible bands on toy one-dimensional regression, outperforming mean-field VI and MC Dropout in uncertainty calibration (Pearce et al., 2018). This is consistent with their theoretical positioning as a wide-network posterior sampler rather than merely a heuristic diversity device.

Bridge pre-assessment gives a concrete example of calibration-oriented surrogate evaluation. On a held-out test slice in the safety-critical range $p(w)=\mathcal N(w;0,\sigma_w^2 I)$8, the three independent BNN heads achieve RMSE $p(w)=\mathcal N(w;0,\sigma_w^2 I)$9, MAPE $p(w\mid D)\propto p(D\mid w)p(w)$0, and $p(w\mid D)\propto p(D\mid w)p(w)$1 for steel bending; RMSE $p(w\mid D)\propto p(D\mid w)p(w)$2, MAPE $p(w\mid D)\propto p(D\mid w)p(w)$3, and $p(w\mid D)\propto p(D\mid w)p(w)$4 for concrete bending; and RMSE $p(w\mid D)\propto p(D\mid w)p(w)$5, MAPE $p(w\mid D)\propto p(D\mid w)p(w)$6, and $p(w\mid D)\propto p(D\mid w)p(w)$7 for shear. Post-hoc scaling $p(w\mid D)\propto p(D\mid w)p(w)$8 is used because the initial BNN is described as mildly overconfident, and calibration curves after scaling are reported to align nominal and empirical coverage within $p(w\mid D)\propto p(D\mid w)p(w)$9 (Kuhn et al., 29 Sep 2025).

6. Limitations, controversies, and active directions

Several limitations recur across the literature. One is prior specification. In mechanics surrogate modeling, the cited work argues that the high dimensionality and non-physicality of the NN parameter space, together with the complex relationship between weights and outputs, makes both prior design and posterior inference challenging; it further argues that neglecting weight–weight correlations impairs the transfer of functional prior information into parameter space (Ghorbanian et al., 2024). This concern is closely related to the earlier objection that plain neural-network ensembles depart from Bayesian methodology even when they are used as uncertainty estimators, motivating anchored alternatives (Pearce et al., 2018).

A second limitation is optimization and model-selection sensitivity. B-INN notes that the choice of $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$0, length scales $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$1, and mode count $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$2 can affect convergence of alternating-direction updates, that block coordinate updates can be prone to local optima, and that basis functions are currently fixed (Park et al., 30 Jan 2026). Multi-modal conjugate-last-layer surrogates rely on Gaussian residuals and do not directly handle non-Gaussian or discrete outputs; very high-dimensional outputs may require PCA or other dimensionality reduction before Wishart-based covariance modeling, and a multi-modal acquisition function for BO remains to be developed (Taylor et al., 26 Sep 2025).

Sampling-based BNN surrogates face their own costs. In the high-dimensional PDE framework, architecture choice, Fourier-feature scales, and HMC tuning parameters $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$3 still require empirical adjustment, and fitting the Gaussian mixture for the joint data density may itself become challenging in very high data dimensions (Jung et al., 2022). In GPBNN, low-fidelity GP fitting retains an $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$4 component, while BNN posterior sampling is feasible because $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$5 and $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$6 in the reported setting (Kerleguer et al., 2023).

Another controversy concerns how often BNN surrogates actually dominate GP baselines. The CCS study explicitly recommends always comparing against a GP baseline in low dimensions and reserving costly MCMC for small-scale problems or cases where posterior fidelity is critical (Fotias et al., 29 Jul 2025). The broader BO study reaches a compatible conclusion in a different form: no single surrogate dominates across tasks, and the choice of surrogate should reflect the problem’s inductive biases rather than any universal preference for BNNs (Li et al., 2023).

Research directions in the cited material are correspondingly pragmatic. Gradient-informed BNNs extend BO surrogates to use local derivative information through automatic differentiation and show faster BO convergence as the number of decision variables increases (Makrygiorgos et al., 14 Apr 2025). Surrogate-assisted parallel tempering replaces some exact likelihood evaluations by a learned likelihood surrogate in order to reduce the cost of Bayesian neural learning; Chandra et al. report substantial time savings but also note that the surrogate tends to under-estimate $\nabla y_i=\nabla_x f(x_i;w)+\eta_i$7 systematically, especially on high-dimensional tasks, and propose local surrogate models, time-series surrogates, GP or deep-kernel surrogates for the likelihood model, and adaptive scheduling of surrogate usage as future directions (Chandra et al., 2018). Taken together, these directions suggest continued pressure toward hybrid constructions that preserve uncertainty quality while reducing the computational cost of posterior approximation.