Efficient Bayesian Physics Informed Neural Networks for Inverse Problems
  via Ensemble Kalman Inversion
Bayesian Physics Informed Neural Networks (B-PINNs) have gained significant
attention for inferring physical parameters and learning the forward solutions
for problems based on partial differential equations. However, the
overparameterized nature of neural networks poses a computational challenge for
high-dimensional posterior inference. Existing inference approaches, such as
particle-based or variance inference methods, are either computationally
expensive for high-dimensional posterior inference or provide unsatisfactory
uncertainty estimates. In this paper, we present a new efficient inference
algorithm for B-PINNs that uses Ensemble Kalman Inversion (EKI) for
high-dimensional inference tasks. We find that our proposed method can achieve
inference results with informative uncertainty estimates comparable to
Hamiltonian Monte Carlo (HMC)-based B-PINNs with a much reduced computational
cost. These findings suggest that our proposed approach has great potential for
uncertainty quantification in physics-informed machine learning for practical
applications.
Andrew Pensoneault, Xueyu Zhu
Efficient Bayesian Physics Informed Neural Networks for Inverse Problems via Ensemble Kalman Inversion

###### Abstract

Bayesian Physics Informed Neural Networks (B-PINNs) have gained significant attention for inferring physical parameters and learning the forward solutions for problems based on partial differential equations. However, the overparameterized nature of neural networks poses a computational challenge for high-dimensional posterior inference. Existing inference approaches, such as particle-based or variance inference methods, are either computationally expensive for high-dimensional posterior inference or provide unsatisfactory uncertainty estimates. In this paper, we present a new efficient inference algorithm for B-PINNs that uses Ensemble Kalman Inversion (EKI) for high-dimensional inference tasks. By reframing the setup of B-PINNs as a traditional Bayesian inverse problem, we can take advantage of EKI's key features: (1) gradient-free, (2) computational complexity scales linearly with the dimension of the parameter spaces, and (3) rapid convergence with typically \(\mathcal{O}(100)\) iterations. We demonstrate the applicability and performance of the proposed method through various types of numerical examples. We find that our proposed method can achieve inference results with informative uncertainty estimates comparable to Hamiltonian Monte Carlo (HMC)-based B-PINNs with a much reduced computational cost. These findings suggest that our proposed approach has great potential for uncertainty quantification in physics-informed machine learning for practical applications.

B -

Bayesian Physically Informed Neural Networks, Inverse Problems, Ensemble Kalman Inversion, Gradient-free

## 1 Introduction

Many applications in science and engineering can be accurately modeled by partial differential equations (PDEs). These equations often contain parameters corresponding to the physical properties of the system. In practice, these properties are often challenging to measure directly. Inverse problems arise when indirect measurements of the system are used to infer these model parameters. These problems are often ill-posed, with the available measurements being limited and noisy. Traditional approaches for solving these problems can be sensitive to noise and data scarcity. In addition, they often require many runs of sophisticated forward numerical solvers [46], which can be computationally expensive. Furthermore, it is not uncommon that the initial or boundary conditions are missing for real-world applications. In such a case, the traditional numerical forward solver might not even be able to run. To address these issues, there is growing interest in developing more efficient, robust, and flexible alternatives.

Recently, Scientific Machine Learning (Sci-ML) [4, 42, 37, 41], a set of approaches that combine domain-specific scientific and engineering knowledge with powerful machine learning tools, has received much attention. These approaches have proven effective in solving PDE-based inverse problems efficiently. One particularly promising approach is Physics Informed Neural Networks (PINNs) [42, 31]. PINNs construct a neural network approximation of the forward solution to the underlying problem while inferring model parameters simultaneously by minimizing data misfit regularized by the underlying governing PDE residual loss. In addition to obtaining estimates of these quantities, the presence of noise and lack of data make it essential to quantify the impact of uncertainty on the surrogate and parameters, especially in high-stakes
applications [50, 41]. However, standard PINN approaches provide only deterministic estimates. Non-Bayesian approaches to uncertainty quantification (UQ) in PINN, such as neural network dropout [48], have been explored. Despite their efficiency, the estimates provided by these methods tend to be less satisfactory, as pointed out in [48, 30].

Alternatively, several attempts to develop Bayesian PINNs (B-PINNs) have been explored, enabling uncertainty quantification in the neural network approximations and the corresponding physical parameter estimates. These methods utilize Bayesian Neural Networks (BNNs) as surrogates by treating weights of neural networks and physical parameters as random variables. In general, Markov Chain Monte Carlo (MCMC) methods are one of the most popular approaches for Bayesian inference tasks. These methods approximate the posterior distribution of parameters with a finite set of samples. While these approaches provide inference with asymptotic convergence guarantees, the use of standard MCMC methods in deep learning applications is limited by their poor scalability for large neural network architectures and large data sets [40]. In practice, Hamiltonian Monte Carlo (HMC) is the gold standard for inference for Bayesian neural networks [17]. HMC is an MCMC method that constructs and solves a Hamiltonian system from the posterior distribution. This approach enables higher acceptance rates than traditional MCMC methods. In the context of B-PINNs [48], an HMC-based B-PINN has been investigated for inference tasks. Despite this, the computational cost of inference remains high.

Alternatively, Variational Inference (VI) approaches to B-PINNs [48] posit that the posterior distribution of the physical and neural network parameters lives within a family of parameterized distributions and solves a deterministic optimization problem to find an optimal distribution within that family [7]. In general, VI approaches are more efficient than HMC approaches and scale well for large parameter spaces and data sizes. Variational inference methods, however, do not share the same theoretical guarantees as MCMC approaches and only provide estimates within the function space of the family of parameterized distributions. Additionally, the inference quality depends on the space of parameterized densities chosen. For example, in [48], KL divergence-based B-PINNs tend to provide less satisfactory estimates than the corresponding HMC B-PINNs and are less robust to measurement noise.

The particle-based VI approach bridges the gap between VI and MCMC methods, combining the strengths of both approaches to provide more efficient inference than MCMC methods while also providing greater flexibility through non-parametric estimates in contrast with the parametric approximations utilized in VI [47]. In the context of B-PINNs, Stein Variational Gradient Descent (SVGD) has been proposed for inference tasks to reconstruct idealized vascular flows with sparse and noisy velocity data [48]. However, SVGD tends to underestimate the uncertainty of the distribution for high dimensional problems, collapsing to several modes [3, 45].

Recently, Ensemble Kalman inversion (EKI) [29, 32, 28] has been introduced as a particle-based VI approach that uses the Ensemble Kalman filter algorithm to solve traditional Bayesian inverse problems. EKI methods have many appealing properties: they are gradient-free, easily parallelizable, robust to noise, and computationally scale linearly with ensemble size [36]. Additionally, these methods use low-rank approximate Hessian (and gradient) information, allowing for rapid convergence of the methods [49]. In the case of a Gaussian prior and linear Gaussian likelihood [33], Ensemble Kalman methods have asymptotic convergence to the correct posterior distribution in a Bayesian sense. Ensemble Kalman methods are asymptotically biased
when these assumptions are violated, yet they are computationally efficient compared to asymptotically unbiased methods and empirically provide reasonable estimates [19, 8]. Recently, these methods have also been used to train neural networks efficiently [32, 23, 14, 22, 49]. In [32], EKI was first proposed as a derivative-free approach to train the neural networks but primarily for traditional purely data-driven machine learning problems. Recently, a one-shot variant of EKI [22] has been used to learn maximum a posteriori (MAP) estimates of a NN surrogate and model parameters for several inverse problems involving PDEs; however, this approach still requires traditional numerical solvers or discretizations to enforce the underlying physical laws, which can be computationally expensive for large scale complex applications. Additionally, while EKI has been traditionally used to obtain MAP estimates of unknown parameters, recent works have begun to investigate the use of the EKI methods to efficiently draw approximate samples for Bayesian inference for traditional Bayesian inverse problems [8, 19, 25, 26, 15].

Motivated by the recent advances in EKI, we present a novel efficient inference method for B-PINNs based on EKI. Specifically, we first recast the setup of B-PINNs as a traditional Bayesian inverse problem so that EKI can be applied to obtain approximate posterior estimates. Furthermore, based on the variant of EKI in [32], we present an efficient sampling-based inference approach to B-PINNs. Because our approach inherits the properties of EKI, it provides efficient gradient-free inference and is well suited for large-scale neural networks thanks to the linear computational complexity of the dimension of unknown parameters [36]. Further, unlike the traditional setting for EKI [29], this approach replaces the expensive numerical forward solver with a NN surrogate trained on the measurements and the underlying physics laws to jointly estimate the physical parameters and forward solution with the corresponding uncertainty estimation. Because the trained neural network surrogate is cheap to evaluate, our proposed approach can draw a large number of samples efficiently and is thus expected to reduce sampling errors significantly. Through various classes of numerical examples, we show that the EKI method can efficiently infer physical parameters and learn valuable uncertainty estimates in the context of B-PINNs. Furthermore, the empirical results show that EKI B-PINNs can deliver comparable inference results to HMC B-PINNs while requiring significantly less computational time. To our best knowledge, this is the _first attempt_ to use EKI for efficient inference under the context of B-PINNs.

The rest of the paper is organized as follows: in section 2, we first introduce the setup of the PDE-based inverse problem and briefly review B-PINNs. In section 3, we first reframe the problem in terms of a Bayesian Inverse problem and introduce EKI under the context of Bayesian inverse problems. We then introduce the building blocks of our proposed EKI B-PINNs framework. We demonstrate the performance of this approach via several numerical examples in section 4. Finally, we conclude in section 5.

## 2 Problem Setup and Background

We consider the following partial differential equation (PDE):

\[\mathcal{N}_{x}(u(\mathbf{x});\mathbf{\lambda}) =f(\mathbf{x}) \mathbf{x}\in\Omega, \tag{1}\] \[\mathcal{B}_{x}(u(\mathbf{x});\mathbf{\lambda}) =b(\mathbf{x}) \mathbf{x}\in\partial\Omega, \tag{2}\]

where \(\mathcal{N}_{x}\) and \(\mathcal{B}_{x}\) denote the differential and boundary operators, respectively. The spatial domain \(D\subseteq\mathbb{R}^{d}\) has boundary \(\Gamma\), and \(\mathbf{\lambda}\in\mathbb{R}^{N_{\lambda}}\) represents a vector of unknown physical parameters. The forcing function \(f(\mathbf{x})\) and boundary function \(b(\mathbf{x})\) are given,
and \(\mathbf{u}(\mathbf{x})\) is the solution of the PDE. For time-dependent problems, we consider time \(t\) as a component of \(\mathbf{x}\) and consider domain \(\Omega\) and boundary \(\partial\Omega\) to additionally contain the temporal domain and initial boundary, respectively.

In this setting, we additionally have access to \(N_{u}\) measurements \(\mathcal{D}_{u}=\{(\mathbf{x}_{u}^{i},u^{i})\}_{i=1}^{N_{u}}\) of the forward solution at various locations. Given the available information, we aim to infer the physical parameters \(\mathbf{\lambda}\) with uncertainty estimation.

### Bayesian Physics Informed Neural Networks (B-PINNs)

Over the last several years, SciML approaches for solving inverse problems have received much attention [4, 31]. One of the most promising approaches is the Physics Informed Neural Networks (PINNs), which approximate the forward solution \(u(\mathbf{x})\) with a fully connected neural network surrogate \(\tilde{u}(\mathbf{x};\mathbf{\theta})\), parameterized by neural network's weight parameter \(\mathbf{\theta}\in\mathbb{R}^{N_{\theta}}\). Denote \(\mathbf{\xi}=[\mathbf{\theta},\mathbf{\lambda}]\) as the concatenation of the neural network and physical parameters. The neural networks are trained by minimizing a weighted sum of the data misfit regularized by initial/boundary data misfit and underlying PDE residual loss (1) over a set of discrete points - "residual points" and "boundary points," respectively, in the domain as follows:

\[\mathcal{D}_{f}=\{(\mathbf{x}_{f}^{i},f(\mathbf{x}_{f}^{i}))\}_{ i=1}^{N_{f}} =\{(\mathbf{x}_{f}^{i},f^{i})\}_{i=1}^{N_{f}} \tag{3}\] \[\mathcal{D}_{b}=\{(\mathbf{x}_{b}^{i},b(\mathbf{x}_{b}^{i}))\}_{ i=1}^{N_{b}} =\{(\mathbf{x}_{b}^{i},b^{i})\}_{i=1}^{N_{b}}, \tag{4}\]

with residual locations \(\mathbf{x}_{f}^{i}\in\Omega\) and boundary locations \(\mathbf{x}_{b}^{i}\in\partial\Omega\). With these notions, the corresponding PINN loss function is defined as follows:

\[\mathcal{L}(\mathbf{\xi})=\omega_{u}\mathcal{L}_{u}(\mathbf{\xi})+\omega_{f}\mathcal{ L}_{f}(\mathbf{\xi})+\omega_{b}\mathcal{L}_{b}(\mathbf{\xi}), \tag{5}\]

where

\[\mathcal{L}_{u}(\mathbf{\xi}) =\frac{1}{N_{u}}\sum_{i=1}^{N_{u}}|u^{i}-\tilde{u}(\mathbf{x}_{u}^ {i};\mathbf{\theta})|^{2}, \tag{6}\] \[\mathcal{L}_{f}(\mathbf{\xi}) =\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}|f^{i}-\mathcal{N}_{x}(\tilde{ u}(\mathbf{x}_{f}^{i};\mathbf{\theta});\mathbf{\lambda})|^{2},\] (7) \[\mathcal{L}_{b}(\mathbf{\xi}) =\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}|b^{i}-\mathcal{B}_{x}(\tilde{ u}(\mathbf{x}_{b}^{i};\mathbf{\theta});\mathbf{\lambda})|^{2}, \tag{8}\]

and \(\omega_{u}\), \(\omega_{f}\), and \(\omega_{b}\) are the weights for each term. In practice, this loss is often minimized with ADAM or L-BFGS optimizers. Standard PINNs typically provide only a deterministic estimate of the target parameters \(\mathbf{\xi}\)[42, 37, 31]. These estimates may be inaccurate and unreliable for problems with small or noisy datasets. Therefore, qualifying the uncertainty in the estimate would be desirable.

To account for the uncertainty, Bayesian Physics-Informed Neural Networks (B-PINNs) (e.g., [43, 48, 35, 2]) have been proposed. B-PINNs are built on Bayesian Neural Networks (BNNs) by treating neural network weights and biases \(\mathbf{\theta}\) and physical parameters \(\mathbf{\lambda}\) as random variables. By Bayes' theorem, the posterior distribution of the parameters \(\xi\) conditions on the forward measurements \(\mathcal{D}_{u}\), the residual points \(\mathcal{D}_{f}\), and the boundary points \(\mathcal{D}_{b}\) can be obtained as follows:

\[p(\mathbf{\xi}|\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b})\propto p(\mathbf{\xi} )p(\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b}|\mathbf{\xi}). \tag{9}\]
The choice of prior distribution \(p(\mathbf{\xi})\) and likelihood \(p(\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b}|\mathbf{\xi})\) will greatly affect the properties of the posterior distribution \(p(\mathbf{\xi}|\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b})\). A typical choice for the prior is to assume independence between the physical parameters \(\mathbf{\lambda}\) and neural network parameters \(\mathbf{\theta}\), i.e. \(p(\mathbf{\xi})=p(\mathbf{\theta})p(\mathbf{\lambda})\). Additionally, the neural network parameters \(\mathbf{\theta}=\{\theta^{i}\}_{i=1}^{N_{\theta}}\) are often assumed to follow independent zero-mean Gaussian distributions, i.e.

\[p(\mathbf{\theta})=\prod_{i=1}^{N_{\theta}}p(\theta^{i}),\quad p(\theta^{i})\sim \mathcal{N}\left(0,\sigma_{\theta}^{i}\right), \tag{10}\]

where \(\sigma_{\theta}^{i}\) is the standard deviation of the corresponding neural network parameter \(\theta^{i}\). For the likelihood, independence between the forward measurements \(\mathcal{D}_{u}\), residual points \(\mathcal{D}_{f}\), and boundary points \(\mathcal{D}_{b}\) is often assumed as follows:

\[p(\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b}|\mathbf{\xi})=p(\mathcal{D}_{u} |\mathbf{\xi})p(\mathcal{D}_{f}|\mathbf{\xi})p(\mathcal{D}_{b}|\mathbf{\xi}). \tag{11}\]

Each term within \(\mathcal{D}_{u}\), \(\mathcal{D}_{f}\), and \(\mathcal{D}_{b}\) is often assumed to follow a Gaussian distribution of the form

\[p(\mathcal{D}_{u}|\mathbf{\xi}) =\prod_{i=1}^{N_{u}}p(u^{i}|\mathbf{\xi}),\quad p(\mathcal{D}_{f}|\bm {\xi})=\prod_{i=1}^{N_{f}}p(f^{i}|\mathbf{\xi}),\quad p(\mathcal{D}_{b}|\mathbf{\xi})= \prod_{i=1}^{N_{b}}p(b^{i}|\mathbf{\xi}), \tag{12}\] \[p(u^{i}|\mathbf{\xi}) =\frac{1}{\sqrt{2\pi\sigma_{\eta_{u}}^{2}}}\exp\left(-\frac{ \left(u^{i}-\tilde{u}(\mathbf{x}_{u}^{i};\mathbf{\theta})\right)^{2}}{2\sigma_{\eta _{u}}^{2}}\right),\] (13) \[p(f^{i}|\mathbf{\xi}) =\frac{1}{\sqrt{2\pi\sigma_{\eta_{f}}^{2}}}\exp\left(-\frac{ \left(f^{i}-\mathcal{N}_{x}(\tilde{u}(\mathbf{x}_{u}^{i};\mathbf{\theta});\mathbf{ \lambda})\right)^{2}}{2\sigma_{\eta_{f}}^{2}}\right),\] (14) \[p(b^{i}|\mathbf{\xi}) =\frac{1}{\sqrt{2\pi\sigma_{\eta_{b}}^{2}}}\exp\left(-\frac{ \left(b^{i}-\mathcal{B}_{x}(\tilde{u}(\mathbf{x}_{b}^{i};\mathbf{\theta});\mathbf{ \lambda})\right)^{2}}{2\sigma_{\eta_{b}}^{2}}\right). \tag{15}\]

Here, \(\sigma_{\eta_{u}}\), \(\sigma_{\eta_{f}}\), and \(\sigma_{\eta_{b}}\) are the standard deviations for the forward measurement, residual point, and boundary point, respectively. Choice of physical parameter prior \(p(\mathbf{\lambda})\) is often problem-dependent, as this distribution represents domain knowledge of the corresponding physical property. Given these choices of prior and likelihood functions, one can construct the corresponding posterior distribution of the BNN. For most BNNs, the closed-form expressions for the posterior distribution are unavailable, and approximate inference methods must be employed. Due to the overparameterized nature of neural networks, the resulting Bayesian inference problem is often high-dimensional for even moderately sized BNNs.

### Hamiltonian Monte Carlo (HMC)

Next, we briefly review Hamiltonian Monte Carlo (HMC), a popular inference algorithm for B-PINNs [48] that serves as a baseline method for our proposed method. Hamiltonian Monte Carlo (HMC) is a powerful method for sampling based posterior inference [18] and has been utilized in the context of inference in Bayesian Neural Networks (BNNs) [39]. HMC employs Hamiltonian dynamics to propose states in parameter space with high acceptance in the Metropolis-Hastings acceptance step. Given the posterior distribution
\(p(\mathbf{\xi}|\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b})=e^{-U(\mathbf{\xi})}\), where \(U\) is the negative log-density of the posterior, we define the Hamiltonian dynamics as follows

\[H(\mathbf{\xi},\mathbf{r})=U(\mathbf{\xi})+\frac{1}{2}\mathbf{r}^{T}M^{-1}\mathbf{r}, \tag{16}\]

where \(\mathbf{r}\in\mathbb{R}^{N_{\xi}}\) is an auxiliary momentum vector, and \(M\in\mathbb{R}^{N_{\xi}\times N_{\xi}}\) is the corresponding mass matrix, often set to the identity \(I_{N_{\xi}}\). Starting from an initial sample of \(\mathbf{\xi}\), the HMC generates proposal samples by resampling momentum \(\mathbf{r}\sim\mathcal{N}(0,M)\) and advancing \((\mathbf{\xi},\mathbf{r})\) through Hamiltonian dynamics

\[\frac{d\mathbf{\xi}}{dt} =-M\mathbf{r}, \tag{17}\] \[\frac{d\mathbf{r}}{dt} =-\nabla U(\mathbf{\xi}). \tag{18}\]

This is often done via Leapfrog integration [18] for \(L\) steps given a step size \(\delta t\). Following this, a Metropolis-Hastings acceptance step is applied to determine if the given sample will be accepted. The details of the HMC are shown in Algorithm 1. The variant HMC B-PINN used in this paper is based on the version in [50], to which we refer interested readers for more details.

```
1:Input: \(\mathbf{\xi}_{0}\) (initial sample), \(\delta t\) (step size), \(L\) (leapfrog steps)
2:for\(i=1,...,J\)do
3:\(\mathbf{\xi}_{i}\leftarrow\mathbf{\xi}_{i-1}\)
4: Sample \(\mathbf{r}_{i}\sim\mathcal{N}(0,M)\)
5:\(\hat{\mathbf{\xi}}_{i}\leftarrow\mathbf{\xi}_{i}\)
6:\(\hat{\mathbf{r}}_{i}\leftarrow\mathbf{r}_{i}\)
7:for\(j=1,...,L\)do
8:\(\hat{\mathbf{r}}_{i}\leftarrow\hat{\mathbf{r}}_{i}-\frac{\delta t}{2}\nabla U(\hat{ \mathbf{\xi}}_{i})\)
9:\(\hat{\mathbf{\xi}}_{i}\leftarrow\hat{\mathbf{\xi}}_{i}+\delta tM^{-1}\hat{\mathbf{r}}_{i}\)
10:\(\hat{\mathbf{r}}_{i}\leftarrow\hat{\mathbf{r}}_{i}-\frac{\delta t}{2}\nabla U(\hat{ \mathbf{\xi}}_{i})\)
11:endfor
12: Sample \(p\sim\mathcal{U}(0,1)\)
13:\(\alpha\leftarrow\min[1,\exp(H(\hat{\mathbf{\xi}}_{i},\hat{\mathbf{r}}_{i})-H(\mathbf{\xi} _{i},\mathbf{r}_{i}))]\)
14:if\(p<\alpha\)then
15:\(\mathbf{\xi}_{i}\leftarrow\hat{\mathbf{\xi}}_{i}\)
16:endif
17:endfor
18:Return: \(\mathbf{\xi}_{1},...,\mathbf{\xi}_{J}\) (Posterior samples)
```

**Algorithm 1** Hamiltonian Monte Carlo (HMC)

From the HMC algorithm, we obtain a set of approximate samples from the B-PINNs posterior distribution \(p(\mathbf{\xi}|\mathcal{D}_{u},\mathcal{D}_{f},\mathcal{D}_{b})\). We shall use these samples to obtain uncertainty estimates of the approximate forward solution and physical parameters. Denote \(\bar{\lambda}\) and \(\bar{u}\) to be the sample mean of physical parameter \(\lambda\) and forward surrogate \(\tilde{u}\), respectively. Additionally, we denote the corresponding sample standard deviations \(s_{\lambda}\) and \(s_{\tilde{u}}\). We compute the sample statistics over the \(J\) samples \(\{\lambda_{j}\}_{j=1}^{J}\) and
\(\{\tilde{u}(\mathbf{x};\mathbf{\theta}_{j})\}_{j=1}^{J}\) obtain from the HMC as follows:

\[\bar{\lambda} =\frac{1}{J}\sum_{j=1}^{J}\lambda_{j},\quad\bar{u}(\mathbf{x})= \frac{1}{J}\sum_{j=1}^{J}\tilde{u}(\mathbf{x};\mathbf{\theta}_{j}), \tag{19}\] \[s_{\lambda} =\sqrt{\frac{\sum_{j=1}^{J}\left(\lambda_{j}-\bar{\lambda} \right)^{2}}{J-1}},\quad s_{\tilde{u}}(\mathbf{x})=\sqrt{\frac{\sum_{j=1}^{J} \left(\tilde{u}(\mathbf{x};\mathbf{\theta}_{j})-\bar{u}(\mathbf{x})\right)^{2}}{J- 1}}. \tag{20}\]

## 3 Ensemble Kalman Inversion-based B-PINNs

In this section, we first briefly review Ensemble Kalman Inversion (EKI) as an efficient method for solving Bayesian inverse problems. Following that, we present our proposed method, denoted _EKI B-PINNs_. Specifically, we first recast the setup of B-PINNs in the traditional Bayesian inverse problem setting and then employ EKI for efficient sampling-based inference.

### Ensemble Kalman Inversion (EKI)

Ensemble Kalman Inversion (EKI) [29, 27, 28, 32] is a popular class of methods that utilize the Ensemble Kalman Filter (EnKF) [20] in the context of traditional inverse problems. These methods are derivative-free, easily parallelizable, and scale well in high-dimension inverse problems with ensemble sizes much smaller than the total number of parameters [32]. Assume that the unknown parameters \(\mathbf{\xi}\in\mathbb{R}^{N_{\xi}}\) have a prior distribution \(p(\mathbf{\xi})\) and the observations \(\mathbf{y}\in\mathbb{R}^{N_{y}}\) are related to the parameters through the observation operator \(\mathcal{G}\):

\[\mathbf{y}=\mathcal{G}(\mathbf{\xi})+\mathbf{\eta}, \tag{21}\]

where \(\mathbf{\eta}\in\mathbb{R}^{N_{y}}\) is a zero-mean Gaussian random vector with observation covariance matrix \(R\in\mathbb{R}^{N_{y}\times N_{y}}\), i.e., \(\mathbf{\eta}\sim\mathcal{N}(0,R)\). In the Bayesian context, this problem corresponds to the following posterior distribution:

\[p(\mathbf{\xi}|\mathbf{y})\propto p(\mathbf{\xi})\exp\left(-\frac{\left\|R^{-1/2}( \mathbf{y}-\mathcal{G}(\mathbf{\xi}))\right\|_{2}^{2}}{2}\right). \tag{22}\]

Given a prior of the form \(p(\mathbf{\xi})\sim\mathcal{N}(\mathbf{\xi}_{0},C_{0})\), the posterior becomes

\[p(\mathbf{\xi}|\mathbf{y})\propto\exp\left(-\frac{\left\|C_{0}^{-1/2}(\mathbf{\xi}_{0} -\mathbf{\xi})\right\|_{2}^{2}+\left\|R^{-1/2}(\mathbf{y}-\mathcal{G}(\mathbf{\xi})) \right\|_{2}^{2}}{2}\right), \tag{23}\]

For weakly nonlinear systems, approximate samples from (23) can be obtained by minimizing an ensemble of loss functions of the form

\[f(\mathbf{\xi}|\mathbf{\xi}_{j},\mathbf{y}_{j})=\frac{1}{2}\left\|C_{0}^{-1/2}(\mathbf{\xi }_{j}-\mathbf{\xi})\right\|_{2}^{2}+\frac{1}{2}\Big{\|}R^{-1/2}(\mathbf{y}_{j}- \mathcal{G}(\mathbf{\xi}))\Big{\|}_{2}^{2}, \tag{24}\]

where \(\mathbf{\xi}_{j}\sim\mathcal{N}(\mathbf{\xi}_{0},C_{0})\) and \(\mathbf{y}_{j}\sim\mathcal{N}(\mathbf{y},R)\)[15, 21]. EnKF-based methods such as EKI can be derived as an approximation to the minimizer of an ensemble of loss functions of the form (24) [21].

In practice, EKI and its variants consider the following artificial dynamics state-space model formulation based on the original Bayesian inverse problem (21) so that
the EnKF update equations can be applied:

\[\mathbf{\xi}_{i} =\mathbf{\xi}_{i-1}+\mathbf{\epsilon}_{i},\quad\mathbf{\epsilon}_{i}\sim\mathcal{ N}(0,Q), \tag{11}\] \[\mathbf{y}_{i} =\mathcal{G}(\mathbf{\xi}_{i})+\mathbf{\eta}_{i},\quad\mathbf{\eta}_{i}\sim \mathcal{N}(0,R), \tag{12}\]

where \(\mathbf{\epsilon}_{i}\) is an artificial parameter noise term with the artificial parameter covariance \(Q\in\mathbb{R}^{N_{\xi}\times N_{\xi}}\) and \(\mathbf{\eta}_{i}\) represents the observation error with the observation covariance \(R\in\mathbb{R}^{N_{y}\times N_{y}}\). Given an initial ensemble of \(J\) ensemble members \(\{\mathbf{\xi}_{0}^{(j)}\}_{j=1}^{J}\), the iterative EKI methods iteratively correct the ensemble \(\{\mathbf{\xi}_{i}^{(j)}\}_{j=1}^{J}\) based on Kalman update equations similar to [32]:

\[\mathbf{\hat{\xi}}_{i}^{(j)} =\mathbf{\xi}_{i-1}^{(j)}+\mathbf{\epsilon}_{i}^{(j)},\quad\mathbf{\epsilon} _{i}^{(j)}\sim\mathcal{N}(0,Q), \tag{13}\] \[\hat{\mathbf{y}}_{i}^{(j)} =\mathcal{G}(\mathbf{\xi}_{i}^{(j)}),\] (14) \[\mathbf{\xi}_{i}^{(j)} =\mathbf{\hat{\xi}}_{i}^{(j)}+C_{i}^{\hat{\xi}\hat{y}}(C_{i}^{\hat{y }\hat{y}}+R)^{-1}(\mathbf{y}-\hat{\mathbf{y}}_{i}^{(j)}+\mathbf{\eta}_{i}^{(j)}), \quad\mathbf{\eta}_{i}^{(j)}\sim\mathcal{N}(0,R), \tag{15}\]

where \(C_{i}^{\hat{y}\hat{y}}\) and \(C_{i}^{\hat{\xi}\hat{y}}\) are the sample covariance matrices defined as follows:

\[C_{i}^{\hat{y}\hat{y}} =\frac{1}{J-1}\sum_{j=1}^{J}(\hat{\mathbf{y}}_{i}^{(j)}-\bar{ \mathbf{y}}_{i})(\hat{\mathbf{y}}_{i}^{(j)}-\bar{\mathbf{y}}_{i})^{T}, \tag{16}\] \[C_{i}^{\hat{\xi}\hat{y}} =\frac{1}{J-1}\sum_{j=1}^{J}(\mathbf{\hat{\xi}}_{i}^{(j)}-\bar{\mathbf{ \xi}}_{i})(\hat{\mathbf{y}}_{i}^{(j)}-\bar{\mathbf{y}}_{i})^{T}. \tag{17}\]

Here, \(\bar{\mathbf{\xi}}_{i}\) and \(\bar{\mathbf{y}}_{i}\) are the corresponding sample average of prior ensembles \(\{\mathbf{\hat{\xi}}_{i}^{(j)}\}\) and \(\{\hat{\mathbf{y}}_{i}^{(j)}\}\). We remark that for many EKI variants, \(Q=0\), which may lead to the ensemble collapsing [13]. This ensemble collapse is not observed for positive definite observation error covariance \(Q\), and desirable convergence properties with reasonable uncertainty estimates have been shown for traditional Bayesian inverse problems [26, 25]. In the case of Gaussian prior and linear measurement operator, convergence to the correct Bayesian posterior corresponding to the artificial dynamics (11)-(12) can be shown [33]. However, these assumptions are typically not satisfied, and thus the corresponding estimates will be biased. Nonetheless, empirical evidence suggested that EKI can still provide reasonable posterior estimates even when these assumptions are violated [8].

### Eki B-Pinn

While EKI-based methods have often been applied to estimate physical model parameters (e.g., [29, 11, 12]) under the context of traditional Bayesian inverse problems, these approaches often rely on existing numerical forward solvers or the corresponding discrete operators. As the EKI requires multiple evaluations of the numerical forward solvers over EKI iterations, this can be computationally expensive for large-scale complex applications. In contrast, our approach (EKI B-PINN) learns a neural network surrogate and infers the model parameters simultaneously without the need for a traditional numerical solver. Combining the inexpensive forward surrogate with EKI allows us to efficiently explore the high-dimensional posterior with larger ensemble sizes.

To employ EKI under the context of B-PINNs, we first recast the setup of B-PINNs in section 2.1 by interpreting the corresponding notation in the context of EKI. Recall the notation introduced in B-PINNs section 2.1: the forward measurement
vector \(\mathbf{u}=\{u^{i}\}_{i=1}^{N_{u}}\), the residual vector \(\mathbf{f}=\{f^{i}\}_{i=1}^{N_{f}}\) and boundary vector \(\mathbf{b}=\{b^{i}\}_{i=1}^{N_{b}}\) from datasets \(\mathcal{D}_{u}=\{(\mathbf{x}_{u}^{i},u^{i})\}_{i=1}^{N_{u}}\), \(\mathcal{D}_{f}=\{(\mathbf{x}_{f}^{i},f^{i})\}_{i=1}^{N_{f}}\) and \(\mathcal{D}_{b}=\{(\mathbf{x}_{b}^{i},b^{i})\}_{i=1}^{N_{b}}\). We utilize the concatenated physical and neural network parameters \(\boldsymbol{\xi}=[\boldsymbol{\lambda},\boldsymbol{\theta}]\). With the approximate solution \(\tilde{u}(\mathbf{x};\boldsymbol{\theta})\) parameterized by \(\boldsymbol{\theta}\), the PDE operator \(\mathcal{N}_{x}\), and the boundary operator \(\mathcal{B}_{x}\), we can define the corresponding observation operator \(\mathcal{G}_{u}\), \(\mathcal{G}_{f}\), and \(\mathcal{G}_{b}\):

\[\mathcal{G}_{u}(\boldsymbol{\xi}) =[\tilde{u}(\mathbf{x}_{u}^{1};\boldsymbol{\theta}),...,\tilde{u }(\mathbf{x}_{u}^{N_{u}};\boldsymbol{\theta})], \tag{3.12}\] \[\mathcal{G}_{f}(\boldsymbol{\xi}) =[\mathcal{N}_{x}(\tilde{u}(\mathbf{x}_{f}^{1};\boldsymbol{\theta });\boldsymbol{\lambda}),...,\mathcal{N}_{x}(\tilde{u}(\mathbf{x}_{f}^{N_{f}} ;\boldsymbol{\theta});\boldsymbol{\lambda})],\] (3.13) \[\mathcal{G}_{b}(\boldsymbol{\xi}) =[\mathcal{B}_{x}(\tilde{u}(\mathbf{x}_{u}^{1};\boldsymbol{ \theta});\boldsymbol{\lambda}),...,\mathcal{B}_{x}(\tilde{u}(\mathbf{x}_{b}^{ N_{b}};\boldsymbol{\theta});\boldsymbol{\lambda})]. \tag{3.14}\]

Given these notations, we now define our measurement vector \(\mathbf{y}\) and the corresponding measurement operator \(G(\boldsymbol{\xi})\) under the context of EKI as follows:

\[\mathbf{y} =[\mathbf{u},\mathbf{f},\mathbf{b}], \tag{3.15}\] \[\mathcal{G}(\boldsymbol{\xi}) =[\mathcal{G}_{u}(\boldsymbol{\xi}),\mathcal{G}_{f}(\boldsymbol{ \xi}),\mathcal{G}_{b}(\boldsymbol{\xi})]. \tag{3.16}\]
After identifying each component in the EKI setting, we employ the version of EKI in (3.7)-(3.9) to infer the parameters of B-PINNs. Following that, we can then compute the sample means and standard deviations of ensemble \(\{\lambda^{(j)}\}_{i=1}^{J}\) of the physical parameter \(\lambda\) and ensemble of approximate solutions \(\{\tilde{u}(\mathbf{x};\boldsymbol{\theta^{(j)}})\}_{j=1}^{J}\) via (2.19)-(2.20) as described in the section 2.2.

We remark that in EKI B-PINN, the underlying physics laws are enforced as soft constraints on the residual points. This approach contrasts with the standard EKI measurement operator, which enforces the physics law as hard constraints imposed by the traditional numerical solvers.

**Choice of Covariance matrices \(Q\) and \(R\).** The choice of observation covariance \(R\) and parameter evolution \(Q\) are often important in the Ensemble Kalman type algorithms, including the EKI algorithms. In the Ensemble Kalman methods literature, several attempts to automatically estimate these covariance matrices have also been suggested (e.g., [6, 34, 1]), however, there is no standard approach to estimating these matrices. In practice, \(R\) is often assumed to be known or estimated empirically from the instrument error and representation error between the states and observations [44]. The matrix \(Q\) is more difficult to estimate, as its dimension is often much larger than the number of available measurements.

In this paper, we assume \(R\) corresponds to the covariance matrix chosen for the Gaussian likelihood \(p(\mathbf{y}|\boldsymbol{\xi})\) in Section 2.1, i.e.,

\[R=\begin{bmatrix}\sigma_{\eta_{u}}^{2}I_{N_{u}}&0&0\\ 0&\sigma_{\eta_{f}}^{2}I_{N_{f}}&0\\ 0&0&\sigma_{\eta_{b}}^{2}I_{N_{b}}\end{bmatrix}. \tag{3.17}\]

The choice of \(Q\) is important in this setting as it prevents the collapse of the ensemble and improves the uncertainty quantification in the EKI. In this study, we assume \(Q\) takes the form

\[Q=\begin{bmatrix}\sigma_{\theta}^{2}I_{N_{\theta}}&0\\ 0&\sigma_{\lambda}^{2}I_{N_{\lambda}}\end{bmatrix}, \tag{3.18}\]

where \(\sigma_{\theta}\) and \(\sigma_{\lambda}\) are standard deviations associated with the neural weight parameters \(\boldsymbol{\theta}\) and the physical parameters \(\boldsymbol{\lambda}\), respectively.

**Stopping Criterion.** We consider a stopping criterion based on the Discrepancy Principle [27, 29, 24, 38], originally utilized as a stopping criterion in iterative regularized Gauss-Newton methods, which has been used as a common choice in many EKI formulations. The discrepancy principle suggests an acceptable choice of \(\boldsymbol{\xi}\) for the inverse problem (3.1) can be made when

\[||R^{-1/2}(y-\mathcal{G}(\boldsymbol{\xi}))||\leq||R^{-1/2}(y- \mathcal{G}(\boldsymbol{\xi}^{\dagger}))||, \tag{3.19}\]

where \(\boldsymbol{\xi}^{\dagger}\) is the true solution of (3.1). This choice avoids instabilities in the solution and provides a criterion for stopping the EKI iteration. As the \(\boldsymbol{\xi}^{\dagger}\) is not known, often it is stated as

\[||R^{-1/2}(y-G(\boldsymbol{\xi}))||\leq\eta \tag{3.20}\]

where \(\eta>0\) is some stopping threshold. In the context of EKI, a sample mean-based discrepancy principle stopping criterion has been employed [27]:

\[\left\|R^{-1/2}\left(y-\frac{1}{J}\sum_{j=1}^{J}\mathcal{G}( \boldsymbol{\xi_{i}^{(j)}})\right)\right\|\leq\eta. \tag{3.21}\]
The choice of \(\eta\) is dependent on \(R\), which can be an issue if a reasonable choice of \(R\) is not clear, such as in the case with residual points \(\mathcal{D}_{f}\) and boundary points \(\mathcal{D}_{b}\) in the context of B-PINNs. Instead, we consider the relative change in the discrepancy over several iterations. By defining \(D_{i}\) as the \(i\)th discrepancy metric as follows:

\[D_{i}=\left\|R^{-1/2}\left(y-\frac{1}{J}\sum_{j=1}^{J}\mathcal{G}(\mathbf{\xi_{i}^{ (j)}})\right)\right\|, \tag{3.22}\]

We shall stop at the iteration when the relative improvement of \(D_{i}\) over a fixed iteration window of length \(W\) does not improve by more than \(\tau\), i.e.,

\[\max_{j\in\{i-W,\dots,i\}}\frac{|D_{j}-D_{i}|}{D_{i}}<\tau. \tag{3.23}\]

**Remark 3.1**: _If an EKI iteration results in a failure state which occurs in the first several iterations, we can reinitialize the ensemble from the prior \(\mathbf{\xi}_{i}^{(j)}\sim p(\mathbf{\xi})\)._

**Complexity Analysis.** At each EKI iteration in Algorithm 2 has following computational complexity of \(\mathcal{O}(JN_{y}N_{\xi}+N_{y}^{3}+JN_{y}^{2})\) detailed as follows:

* Construction of matrix \(C^{\hat{y}\hat{y}}\) in (3.10)
* Construction of matrix \(C^{\hat{\xi}\hat{y}}\) in (3.11)
* Evaluation of \(C^{\hat{\xi}\hat{y}}(C^{\hat{y}\hat{y}}+R)^{-1}(\hat{\mathbf{y}}_{i}^{(j)}+\mathbf{\eta }_{i}^{(j)})\) in (3.9).

As a result, the computational complexity of EKI grows linearly with both parameter dimension and ensemble size and thus can easily scale to high-dimensional inverse problems efficiently. Nevertheless, the method does scale cubically with data dimension and hence is computationally prohibitive when naively applied to large data sets. In this scenario, mini-batching, as in standard neural network optimization problems, can be used to reduce the computational cost [32]. Furthermore, it is important to note that the EKI requires storing an ensemble of \(J\) neural network parameter sets, which can be memory-demanding for large ensemble sizes and large networks. Dimension reduction techniques might be helpful to address this issue [16].

## 4 Numerical Examples

In this section, we shall demonstrate the applicability and performance of EKI B-PINNs via various numerical examples. We also compare our approach with a variant of the HMC B-PINN to assess the inference accuracy and computational efficiency of our proposed method.

For each example, we generated a synthetic data set, \(\mathcal{D}_{u}\), by solving the corresponding problem and corrupting the solution with i.i.d zero-mean Gaussian noise, \(\mathcal{N}(0,\sigma_{u})\). To demonstrate the robustness of the EKI method, we considered two noise levels: \(\sigma_{u}=0.1\) and \(\sigma_{u}=0.01\). We generate residual points \(\mathcal{D}_{f}\) by evaluating \(f\) at locations generated using Latin hypercube sampling over the problem domain for 2D problems and equally spaced over the domain for 1D problems. For 2D problems, boundary points \(\mathcal{D}_{b}\) are placed equally spaced over the boundary. Boundary and residual points are assumed to be noise-free in this paper.

For all examples, the neural network architecture used for both B-PINNs consists of 2 hidden layers with 50 neurons in each layer and the tanh activation function. Neural network parameter dimension \(N_{\theta}\) for each type of example can be seen in Table 1. For each problem, we assume the standard deviations for the boundary and residual likelihood in (3.17): \(\sigma_{\eta_{b}}=0.01\) and \(\sigma_{\eta_{f}}=0.01\) for the B-PINNs unless otherwise specified. Furthermore, we assume that the measurement noise level \(\sigma_{u}\) is
```
1:Input: \(\mathbf{y}\) (observations), \(Q\) (parameter covariance), \(R\) (observation covariance), \(W\) (stopping window), \(\tau\) (stopping threshold)
2:Initialize prior samples for \(j=1,...,J\): \[\boldsymbol{\xi}_{0}^{(j)}\sim p(\boldsymbol{\xi}).\]
3:for\(i=1,...,I\)do
4: Obtain \(i\) th prior parameter and measurement ensembles for \(j=1,...,J\): \[\boldsymbol{\epsilon}_{i}^{(j)} \sim\mathcal{N}(0,Q).\] \[\boldsymbol{\hat{\xi}}_{i}^{(j)} =\boldsymbol{\xi}_{i}^{(j)}+\boldsymbol{\epsilon}_{i}^{(j)}.\] \[\boldsymbol{\hat{y}}_{i}^{(j)} =\mathcal{G}_{i}(\boldsymbol{\hat{\xi}}_{i}^{(j)}).\]
5: Evaluate sample mean and covariance terms: \[\boldsymbol{\bar{\xi}}_{i} =\frac{1}{J}\sum_{j=1}^{J}\boldsymbol{\hat{\xi}}_{i}^{(j)}.\] \[\bar{\mathbf{y}}_{i} =\frac{1}{J}\sum_{j=1}^{J}\boldsymbol{\hat{y}}_{i}^{(j)}.\] \[C_{i}^{\hat{y}\hat{y}} =\frac{1}{J-1}\sum_{j=1}^{J}(\boldsymbol{\hat{y}}_{i}^{(j)}- \bar{\mathbf{y}}_{i})(\boldsymbol{\hat{y}}_{i}^{(j)}-\bar{\mathbf{y}}_{i})^{ T}.\] \[C_{i}^{\hat{\xi}\hat{y}} =\frac{1}{J-1}\sum_{j=1}^{J}(\boldsymbol{\hat{\xi}}_{i}^{(j)}- \boldsymbol{\bar{\xi}}_{i})(\boldsymbol{\hat{y}}_{i}^{(j)}-\bar{\mathbf{y}}_{ i})^{T}.\]
6: Update posterior ensemble for \(j=1,...,J\): \[\boldsymbol{\eta}_{i}^{(j)} \sim\mathcal{N}(0,R),\] \[\boldsymbol{\xi}_{i}^{(j)} =\boldsymbol{\hat{\xi}}_{i}^{(j)}+C_{i}^{\hat{\xi}y}(C_{i}^{yy} +R)^{-1}(\mathbf{y}-\hat{\mathbf{y}}_{i}^{j}+\boldsymbol{\eta}_{i}^{(j)}).\]
7: Check discrepancy \[D_{i}=\left\|R^{1/2}\left(y-\frac{1}{J}\sum_{j=1}^{J}\mathcal{G}(\boldsymbol{ \xi}_{i}^{(j)})\right)\right\|,\]
8:if\(\max\limits_{j\in\{i-W,...,i\}}|D_{j}-D_{i}|/D_{i}<\tau\)then
9:\(I:=i\)
10: Break
11:endif
12:endfor
13:Return: \(\boldsymbol{\xi}_{I}^{(1)},...,\boldsymbol{\xi}_{I}^{(J)}\)
```

**Algorithm 3** Ensemble Kalman Inversion (EKI) B-PINNs
known for each example and set \(\sigma_{\eta_{u}}=\sigma_{u}\). Finally, unless otherwise specified, we use physical parameter prior \(\mathbf{\lambda}\sim\mathcal{N}(0,I_{N_{\lambda}})\) for both B-PINNs.

For the EKI B-PINNs, we choose an ensemble size of \(J=1000\) and stopping criterion parameters from (3.23): \(W=25\) and \(\tau=0.05\). Additionally, the artificial dynamics standard deviations for the parameters from (3.18) are chosen to be \(\sigma_{\lambda}=0.1\) and \(\sigma_{\theta}=0.002\). For the HMC B-PINNs, the leapfrog step \(L=50\) and the initial step size \(\delta t=0.1\) adaptively tuned to reach an acceptance rate of \(60\%\) during burn-in steps as in [50]. We draw a total of \(1000\) samples following \(1000\) burn-in steps.

**Metrics.** We examine the accuracy of the forward solution approximation and physical parameter estimation from the B-PINNs in the following metrics over an independent test set \(\{\mathbf{x}_{t}^{i}\}_{i=1}^{N_{t}}\):

\[e_{u}=\sqrt{\frac{\sum_{i=1}^{N_{t}}\left|u(\mathbf{x}_{t}^{i})-\bar{u}( \mathbf{x}_{t}^{i})\right|^{2}}{\sum_{i=1}^{N_{t}}\left|u(\mathbf{x}_{t}^{i}) \right|^{2}}},\quad e_{\lambda}=\frac{\left|\lambda-\bar{\lambda}\right|}{ \left|\lambda\right|}, \tag{4.1}\]

where \(u(x)\) and \(\lambda\) are the reference solution and the reference physical parameters. The sample means of the forward solution approximation \(\bar{u}\) and the physical parameter \(\bar{\lambda}\) are computed from the B-PINNs as defined in (2.19). The mean predictions are obtained in the final iteration for the EKI B-PINNs and over the \(1000\) samples generated for the HMC B-PINN. Furthermore, we assess the quality of our uncertainty estimates by examining the sample standard deviation of the estimated forward solutions and physical parameters \(\sigma_{\lambda}\) and \(\sigma_{\tilde{u}}\) defined in (2.20).

To demonstrate the efficiency of our proposed method, we compare the walltime of the EKI B-PINN and HMC B-PINN (including both burn-in and training time) experiments averaged over \(10\) trials. The average iterations utilized over \(10\) trials are also presented for the EKI B-PINN tests. All experiments use the JAX [9] library in single-precision floating-point on a single GPU (Tesla T4 GPU with \(16\) GB of memory).

### 1D linear Poisson equation

We first consider the 1D linear Poisson equation motivated by [10] as follows:

\[u_{xx}-k\cos(x)=0,\quad x\in[0,8], \tag{4.2}\] \[u(0)=0,\quad u(8)=k\cos(8), \tag{4.3}\]

where the exact solution is \(u(x)=k\cos(x)\). Assuming the unknown constant \(k=1.0\), we generated \(N_{u}=8\) equally spaced measurements over \([0,8]\) for measurements \(\mathbf{u}\), excluding the boundary points, and \(N_{b}=2\) boundary points at \(x=0\) and \(x=8\). Additionally, \(N_{f}=100\) equally spaced residual points are utilized.

Due to the linearity of the solution with respect to \(k\), one can derive a Gaussian "reference" posterior distribution for parameter \(k\) conditioned on knowledge of

\begin{table}
\begin{tabular}{|l|l|} \hline  & \(\mathbf{\theta}\) size \\ \hline
1D PDE - Examples 4.1, 4.2 & 5251 \\ \hline
2D PDE - Examples 4.3, 4.5 4.6 & 5301 \\ \hline System of ODEs - Example 4.4 & 5353 \\ \hline \end{tabular}
\end{table}
Table 4.1: The number of neural network parameters \(\mathbf{\theta}\) for each example.
the correct solution parameterization \(\tilde{u}(x;k)=k\cos(x)\)[10]. In the case where the forward solution and boundary standard deviations are equal (i.e., \(\sigma_{\eta_{u}}=\sigma_{\eta_{b}}\)), if we denote for simplicity \(\mathcal{D}_{u}\) as containing both forward solution and boundary data, the distribution takes the form:

\[p(k|\mathcal{D}_{u})\propto\exp\left(-\left(\frac{\sum_{i=1}^{N_{u}}(u^{i}-k \cos(x_{u}^{i}))^{2}}{2\sigma_{\eta_{u}}^{2}}+\frac{(k-k_{0})^{2}}{2\sigma_{k} ^{2}}\right)\right), \tag{4.4}\]

where \(k_{0}\), \(\sigma_{k}\) are the prior mean and standard deviation of \(k\).

In this example, we choose \(k_{0}=0,\sigma_{k}=1\) and consider the noise level \(\sigma_{u}=0.01\). We present the "reference" density and density estimates of the estimated posterior of \(k\) in Fig 4. Both approaches deliver comparable inference results for the posterior distribution of \(k\) compared to the "reference" density.

Table 4 provides the mean and one standard deviation of the posterior \(k\) approximations. The uncertainty estimates of \(k\) for both B-PINNs show that the true value of \(k\) is within one standard deviation of the posterior mean. Table 4 shows the mean relative error of the solution and parameters, with corresponding walltime and the number of EKI iterations utilized. From the table, it is clear that the EKI B-PINN can achieve comparable approximation quality to HMC B-PINN but with a 8-fold speed-up.

### 1D nonlinear Poisson equation

We now consider a 1D nonlinear Poisson equation as presented in [48] as follows:

\[\lambda u_{xx}+k\tanh(u)=f,\quad x\in[-0.7,0.7], \tag{4.5}\]

where \(\lambda=0.01\). Here, \(k=0.7\) is the unknown physical parameter to be inferred. By assuming the true solution \(u(x)=\sin^{3}(6x)\), the right-hand side \(f\) and boundary conditions can be analytically constructed. The objective of this problem is to infer the physical parameter \(k\) and the forward solution \(u\), along with corresponding uncertainty estimates, using \(N_{u}=6\) measurements \(\mathbf{u}\) equally spaced over the spatial domain, \(N_{b}=2\) boundary points \(\mathbf{b}\), and \(N_{f}=32\) equally spaced PDE residual points \(\mathbf{f}\).

Table 4 compares the posterior sample mean and one standard deviation of \(k\) obtained from the EKI B-PINN and HMC B-PINN. Both methods accurately capture

\begin{table}
\begin{tabular}{c|c|c c|c} \hline  & & \(e_{u}\) & \(e_{k}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & 0.63\% & 0.18\% & 5.81 seconds & 219 \\  & HMC & 1.43\% & 0.08\% & 46.66 seconds & - \\ \hline \end{tabular}
\end{table}
Table 4.3: Example 4.1: relative errors \(e_{u}\) of the forward solution \(u\) and \(e_{k}\) of parameter \(k\) for the noise level \(\sigma_{u}=0.01\), as well as average walltime. Average EKI iterations for the EKI B-PINN are also reported.
the true value of \(k=0.7\) within one standard deviation of the mean for both noise levels. Figure 4.2 compares the sample mean and standard deviation of the surrogate based on both approaches. The mean of the EKI B-PINN is in good agreement with the reference solution. The one standard deviation bound captures most of the error in the mean prediction.

Table 4.5 reports the relative errors of the forward solution and parameter estimates using EKI B-PINN and HMC B-PINN for two noise levels, along with the corresponding walltime and number of EKI iterations. Both methods are reasonably accurate, with estimates of \(k\) within \(1\%\) error for \(\sigma_{u}=0.01\) and \(5\%\) error for \(\sigma_{u}=0.1\). Additionally, the relative errors of the mean forward solution \(e_{u}\) for both B-PINNs reach similar levels of accuracy, approximately \(2\%\) and \(10\%\) for \(\sigma_{u}=0.01\) and \(\sigma_{u}=0.1\) respectively. Nevertheless, the EKI B-PINN achieves comparable inference results to the HMC B-PINN, but approximately \(9\) times faster.

\begin{table}
\begin{tabular}{c|c|c} \hline  & & \(k\) (mean \(\pm\) std) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & \(0.701\pm 0.006\) \\  & HMC & \(0.701\pm 0.006\) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & \(0.697\pm 0.012\) \\  & HMC & \(0.688\pm 0.022\) \\ \hline \end{tabular}
\end{table}
Table 4.4: Example 4.2: sample mean and standard deviation of parameter \(k\) for EKI B-PINN and HMC B-PINN for \(\sigma_{u}=0.01,0.1\) noise levels. The true value of \(k\) is \(0.7\).

Figure 4.1: Example 4.1: posterior distributions of \(k\) for EKI B-PINN and HMC B-PINN at \(\sigma_{u}=0.01\) with the corresponding reference posterior distribution.
\begin{table}
\begin{tabular}{c|c|c|c} \hline  & & \(k\) (mean\(\pm\) std) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & \(0.999\pm 0.006\) \\  & HMC & \(0.996\pm 0.006\) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & \(0.988\pm 0.017\) \\  & HMC & \(1.023\pm 0.032\) \\ \hline \end{tabular}
\end{table}
Table 4.6: Example 4.3: sample mean and standard deviation of parameter \(k\) for EKI B-PINN and HMC B-PINN for \(\sigma_{u}=0.01,0.1\) noise levels. The true value of \(k\) is 1.

Figure 4.2: Example 4.2: sample mean and standard deviation of \(u\) for EKI B-PINN and HMC B-PINN for \(\sigma_{u}=0.01,0.1\) noise levels.

\begin{table}
\begin{tabular}{c|c c|c|c|c} \hline  & & \(e_{u}\) & \(e_{k}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & \(1.19\%\) & \(0.19\%\) & \(6.43\) seconds & \(282\) \\  & HMC & \(1.12\%\) & \(0.10\%\) & \(55.55\) seconds & - \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & \(9.32\%\) & \(0.38\%\) & \(6.47\) seconds & \(289\) \\  & HMC & \(8.76\%\) & \(1.73\%\) & \(55.77\) seconds & - \\ \hline \end{tabular}
\end{table}
Table 4.5: Example 4.2: relative errors \(e_{u}\) of the forward solution \(u\) and \(e_{k}\) of parameter \(k\) for the noise levels \(\sigma_{u}=0.01,0.1\), as well as average walltime. Average EKI iterations for the EKI B-PINN are also reported.
### 2D nonlinear diffusion-reaction equation

Next, we examine the 2D nonlinear diffusion-reaction equation in [48] as follows:

\[\lambda\Delta u+ku^{2} =f,\quad(x,y)\in[-1,1]^{2}, \tag{4.6}\] \[u(x,-1) =u(x,1)=0,\] (4.7) \[u(-1,y) =u(1,y)=0, \tag{4.8}\]

where \(\lambda=0.01\) is known. Here, \(k=1\) is an unknown physical parameter, and the ground truth solution \(u(x,y)=\sin(\pi x)\sin(\pi y)\). We construct the source term \(f\) to satisfy the PDE with the given solution. For this problem, we have \(N_{u}=100\) measurements \(\mathbf{u}\) and \(N_{f}=100\) residual points \(\mathbf{f}\) both sampled via Latin Hypercube sampling over the spatial domain. Additionally, we have \(N_{b}=100\) boundary points \(\mathbf{b}\), which are equally spaced over the boundary. As in the previous example, we aim to estimate \(k\) and \(u\) with uncertainty estimates. The solution \(u\) and measurements \(\mathbf{u}\) can be seen in Figure 4.3.

Table 4.6 shows the one standard deviation confidence interval of the B-PINN estimates for parameter \(k\). Notably, the ground truth \(k=1.0\) falls within one standard deviation of the mean estimates for both noise levels. Additionally, Figure 4.4 shows the sample mean, one standard deviation, and the error of the forward surrogate based on both approaches. Both B-PINNs provide reasonably good mean estimates of the true solution for two measurement noise levels. Moreover, the standard deviation for both B-PINNs increases as the measurement noise level increases as expected, indicating that the uncertainty estimates provided are plausible. Although the standard deviations by the EKI B-PINN somewhat differ from those of the HMC B-PINN, both

\begin{table}
\begin{tabular}{c|c|c c|c|c} \hline  & & \(e_{u}\) & \(e_{k}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & 1.12\% & 0.04\% & 2.47 seconds & 53 \\  & HMC & 1.06\% & 0.03\% & 63.26 seconds & - \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & 3.64\% & 1.46\% & 3.41 seconds & 66 \\  & HMC & 2.53\% & 3.73\% & 64.71 seconds & - \\ \hline \end{tabular}
\end{table}
Table 4.7: Example 4.3: relative errors \(e_{u}\) of the forward solution \(u\) and \(e_{k}\) of parameter \(k\) for the noise levels \(\sigma_{u}=0.01,0.1\), as well as average walltime. Average EKI iterations for the EKI B-PINN are also reported.

Figure 4.3: Example 4.3: measurements of forward solution \(\mathbf{u}\) and boundary measurements \(\mathbf{b}\) with solution \(u(x,y)=\sin(\pi x)\sin(\pi y)\).
methods appear to agree on the rough locations of regions with higher uncertainty. For instance, when \(\sigma_{u}=0.01\), For example, when \(\sigma_{u}=0.01\), both methods exhibit large peaks around \((x,y)=(-1,0),(-0.3,0.7),(0.9,0.9),(1,-1)\). Similarly, major peaks can be seen around at \((x,y)=(0.8,0.8),(-0.8,-0.8),(-0.8,0.8),(0.8,-0.8)\) for both methods when \(\sigma_{u}=0.1\).

The relative error of the mean estimates of \(u\) and \(k\) and walltime for the B-PINN methods are presented in Table 4.7. The mean approximations of \(k\) for both B-PINNs show approximation errors around \(1\%\) and \(5\%\) of the ground truth for \(\sigma_{u}=0.01\) and \(\sigma_{u}=0.1\), respectively. Similarly, the sample mean approximation of the forward surrogate \(u\) achieved a relative error less than \(2\%\) and \(5\%\), respectively. The EKI B-PINN approximates the forward solution and physical parameter reasonably well, with mean estimates comparable to the HMC. Nonetheless, the EKI B-PINN provides inference approximately \(18\) times faster than the HMC B-PINN.

Figure 4.4: Example 4.3: (Top Row) the sample mean of forward surrogate for EKI and HMC B-PINNs for different noise levels \(\sigma_{u}\). (Middle Row) The standard deviation of the forward surrogate based on EKI and HMC B-PINNs for different noise levels. (Bottom Row) The absolute difference between the ground truth solution \(u(x,y)=\sin(\pi y)\sin(\pi y)\) and the sample mean of the forward approximation for different noise levels.
### Kraichnan-Orszag system

We next consider the Kraichnan-Orszag model [50] consisting of three coupled nonlinear ODEs describing the temporal evolution of a system composed of several interacting inviscid shear waves:

\[\frac{du_{1}}{dt}-au_{2}u_{3}=0, \tag{11}\] \[\frac{du_{2}}{dt}-bu_{1}u_{3}=0,\] (12) \[\frac{du_{3}}{dt}+(a+b)u_{1}u_{2}=0,\] (13) \[u_{1}(0)=1.0,u_{2}=0.8,u_{3}(0)=0.5, \tag{14}\]

where \(u_{1}\), \(u_{2}\), and \(u_{3}\) are solutions to the above system of ODEs and \(a\) and \(b\) are unknown physical parameters. We choose \(a=b=1\) in this example. We place \(12\) equally spaced data points over \(t\in[1,10]\), and observe \(u_{1}\) and \(u_{3}\) at each of these locations, and \(u_{2}\) at the first \(7\) locations, and thus, \(N_{u}=31\) for this example. We also utilize \(N_{f}=300\) residual points (with \(100\) points equally spaced over \(t\in[0,10]\) for each equation in the system of ODEs). Initial conditions for the ODE are assumed to be unknown, thus \(N_{b}=0\). As in [50], we place a Gaussian prior on \(a\) and \(b\) such that \(a\sim\mathcal{N}(0,2)\) and \(b\sim\mathcal{N}(0,2)\). Our goal is to estimate \(u_{1}\), \(u_{2}\), \(u_{3}\), and parameters \(a\) and \(b\) with uncertainty estimates.

Table 8 shows the mean and one standard deviation of the estimated \(a\) and \(b\). For the small noise level \(\sigma_{u}=0.01\), the parameter estimates capture the truth \(b\) within one standard deviation and \(a\) within two standard deviations of the mean. For \(\sigma_{u}=0.1\) noise level, both estimates find \(b\) within two standard deviations of the mean. However,\(a\) is less well approximated for both B-PINNs due to bias in the mean approximation. Furthermore, Figure 5 presents the approximation of the forward solutions \((u_{1},u_{2},u_{3})\) with the means and one standard deviation bands. For \(\sigma=0.01\), the forward predictions closely match the true solutions with narrow uncertainty bands for the B-PINNs. For the case \(\sigma=0.1\), the deviations between the mean prediction of \(u_{1}\) and the ground truth are more pronounced, particularly at locations where the standard deviation is also larger, suggesting both B-PINNs offer informative uncertainty estimates.

Table 9 presents the mean relative errors for parameter estimations and forward approximations. The mean approximations of the EKI B-PINN are accurate and comparable to those of the HMC B-PINN for both noise levels, with the exception of parameter \(a\) for \(\sigma_{u}=0.1\), where both B-PINNs show less accuracy. However, the EKI method is approximately \(25\) times faster than the HMC method in this example.

\begin{table}
\begin{tabular}{c|c|c c} \hline \hline  & & \(a\) (mean \(\pm\) std) & \(b\) (mean \(\pm\) std) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & \(0.953\pm 0.023\) & \(1.003\pm 0.018\) \\  & HMC & \(0.978\pm 0.025\) & \(1.001\pm 0.016\) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & \(0.847\pm 0.033\) & \(1.027\pm 0.026\) \\  & HMC & \(0.826\pm 0.053\) & \(1.029\pm 0.032\) \\ \hline \hline \end{tabular}
\end{table}
Table 8: Example 4: sample mean and standard deviation of parameters \(a\) and \(b\) for EKI B-PINN and HMC B-PINN for \(\sigma_{u}=0.01,0.1\) noise levels. The true values of \(a=1,b=1\).
\begin{table}
\begin{tabular}{c|c|c c c c c|c|c|c} \hline  & & \(e_{u_{1}}\) & \(e_{u_{2}}\) & \(e_{u_{3}}\) & \(e_{a}\) & \(e_{b}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & 1.04\% & 2.96\% & 1.63\% & 4.67\% & 0.35\% & 2.93 seconds & 71 \\  & HMC & 0.66\% & 1.45\% & 1.06\% & 2.25\% & 0.14\% & 93.34 seconds & - \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & 3.49\% & 5.57\% & 4.45\% & 15.26\% & 2.73\% & 3.51 seconds & 85 \\  & HMC & 3.96\% & 6.49\% & 5.19\% & 17.42\% & 2.88\% & 94.31 seconds & - \\ \hline \end{tabular}
\end{table}
Table 4.9: Example 4.4: relative errors \(e_{u}\) of the forward solution \(u\), \(e_{a}\) of parameter \(a\), and \(e_{b}\) of parameter \(b\) for the noise levels \(\sigma_{u}=0.01,0.1\), as well as average walltime. Average EKI iterations are also reported.

Figure 4.5: Example 4.4: the ground truth, sample mean, standard deviation, and measurements \(\mathbf{u}\) for \(u_{1}\) (top row), \(u_{2}\) (middle row), and \(u_{3}\) (bottom row) for the EKI B-PINN and HMC B-PINN with noise levels \(\sigma_{u}=0.01,0.1\).
### Burgers' Equation

We now consider the following 1D time-dependent Burgers' equation motivated by [42]:

\[u_{t}+uu_{x}-\nu u_{xx} =0, x\in[-1,1], \tag{4.13}\] \[u(x,0) =-\sin(\pi x), x\in[-1,1],\] (4.14) \[u(-1,t) =u(1,t)=0, \tag{4.15}\]

where \(\nu=\frac{0.1}{\pi}\) is an unknown physical parameter, and \(t\in[0,\frac{3}{\pi}]\). A reference solution \(u\) for (4.13) is found in [5], where we evaluate on a \(256\times 100\) equally spaced grid over the domain. We randomly sampled \(N_{u}=100\) measurements \(\mathbf{u}\) from the solution grid and placed \(N_{b}=75\) equally spaced boundary points over the boundary, which is visualized with the solution in Figure 6. Residual points \(\mathbf{f}\) were sampled using Latin hypercube sampling over \((x,t)\in[-1,1]\times[0,\frac{3}{\pi}]\) with \(N_{f}=100\). Instead of directly approximating \(\nu\), we approximate the transformed parameter \(\log\nu\) and place the prior \(\log\nu\sim\mathcal{N}(0,3)\). Additionally, for this example, we choose the standard deviation of the likelihood for the residual \(\sigma_{f}=0.1\). In this example, we aim to infer \(\nu\) and \(u\) with corresponding uncertainty estimates.

Table 10 shows the posterior mean and one standard deviation of the estimated parameter \(\nu\). Both methods provide accurate mean approximations of the parameter, and the true \(\nu\) is contained within the two standard deviation bounds for \(\sigma_{u}=0.01\), and lies just outside the one standard deviation confidence interval from the EKI B-PINN. For \(\sigma_{u}=0.1\), EKI B-PINN does provide a slightly better mean estimate, but it underestimates the uncertainty compared to the HMC.

\begin{table}
\begin{tabular}{c|c|c} \hline  & & \(\nu\) (\(\times 10^{3}\)) (\(\mathrm{mean}\pm\mathrm{std}\)) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & \(32.449\pm 0.602\) \\  & HMC & \(30.957\pm 1.027\) \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & \(33.740\pm 0.808\) \\  & HMC & \(34.174\pm 3.523\) \\ \hline \end{tabular}.

\end{table}
Table 10: Example 4.5: the sample mean and standard deviation of parameter \(\nu\) for EKI B-PINN and HMC B-PINN for \(\sigma_{u}=0.01,0.1\) noise levels. The true value \(\nu=\frac{0.1}{\pi}\approx 0.031831\).

Figure 6: Example 4.5: measurements of forward solution \(\mathbf{u}\) and boundary measurements \(\mathbf{b}\) with solution \(u(x,t)\) over \([-1,1]\times[0,\frac{3}{\pi}]\).
Figure 7 displays the posterior mean and standard deviation of \(u\) estimated by both methods, as well as their corresponding approximation error. As expected, the uncertainty in both B-PINNs increases with the measurement noise level. Notably, larger uncertainty is clustered around \(x=0\), and greater errors are observed in the same region in the error plot, which suggests that the uncertainty estimates are informative as an error indicator.

Table 11 provides mean relative error for \(\nu\) and \(u\). Both methods can achieve reasonably good accuracy with relative errors of less than 2% for both \(\nu\) and \(u\) when \(\sigma_{u}=0.01\) and 8% for \(\sigma=0.1\). The inference results from EKI are obtained approximately 10 times faster than the HMC method.

### Diffusion-Reaction Equation with Source Localization

Finally, we consider the following source inversion problem for a two-dimensional linear diffusion

Figure 7: Example 5: (Top Row) the sample mean of forward surrogate for EKI and HMC B-PINNs for different noise levels \(\sigma_{u}\). (Middle Row) The standard deviation of the forward surrogate based on EKI and HMC B-PINNs for different noise levels. (Bottom Row) The absolute difference between the ground truth solution \(u(x,y)\) and the sample mean of the forward approximation for different noise levels.
reaction equation inspired by the example [48]:

\[-\lambda\Delta u-f_{2} =f_{1},\quad(x,y)\in[0,1]^{2}, \tag{4.16}\] \[u(0,y)=u(1,y) =0,\] (4.17) \[u(x,0)=u(x,1) =0, \tag{4.18}\]

where \(\lambda=0.02\) is the known diffusivity and \(f_{1}(x)=0.1\sin(\pi x)\sin(\pi y)\) is a known forcing. The goal is to infer the location parameters of the field \(f_{2}\), corresponding to three contaminant sources of the unknown location in the following equation:

\[f_{2}(\mathbf{x})=\sum_{i=1}^{3}k_{i}\exp\left(-\frac{||\mathbf{x}-\mathbf{x}_ {i}||^{2}}{2\ell^{2}}\right). \tag{4.19}\]

Here, \(\mathbf{k}=(k_{1},k_{2},k_{3})=(2,-3,0.5)\) and \(\ell=0.15\) are known constants and parameters \(\mathbf{x}_{1}=(0.3,0.3)\), \(\mathbf{x}_{2}=(0.75,0.75)\), \(\mathbf{x}_{3}=(0.2,0.7)\) are parameters to be recovered. The prior distributions on the locations \(\mathbf{x}_{i}=(x_{i},y_{i})\) are chosen to be log-normal, that is \(\log(x_{i})\sim\mathcal{N}(0,1)\) and \(\log(y_{i})\sim\mathcal{N}(0,1)\) for \(i=1,2,3\). The measurements \(\mathbf{u}\) are generated by solving (4.16) using the _solvepde_ function in Matlab with 1893 triangle meshes and sampling randomly \(N_{u}=100\) points from among the nodes and \(N_{b}=100\) equally spaced points along the domain's boundary, shown in Figure 4.8. In addition, \(N_{f}=100\) residual points were obtained via Latin Hypercube Sampling within the domain.

The posterior mean and standard deviation of the locations \(\mathbf{x}_{1}\), \(\mathbf{x}_{2}\) and \(\mathbf{x}_{3}\) in Table 4.12 show an accurate estimation of most of the centers of contaminants, with all

\begin{table}
\begin{tabular}{c|c|c c|c|c} \hline  & & \(e_{u}\) & \(e_{\nu}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & 1.04\% & 1.94\% & 5.90 seconds & 71 \\  & HMC & 0.77\% & 2.75\% & 62.57 seconds & - \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & 3.48\% & 6.00\% & 6.72 seconds & 116 \\  & HMC & 3.07\% & 7.36\% & 62.92 seconds & - \\ \hline \end{tabular}
\end{table}
Table 4.11: Example 4.5: relative errors \(e_{u}\) of the forward solution \(u\) and \(e_{\nu}\) of parameter \(\nu\) for the noise levels \(\sigma_{u}=0.01,0.1\), as well as the average walltime. Average iterations for EKI B-PINN are also reported.

Figure 4.8: Example 4.6: measurements of forward solution \(\mathbf{u}\) and boundary measurements \(\mathbf{b}\) with solution \(u(x,t)\).
values falling within contained within two standard deviations and most within one standard deviation of the mean estimates. The corresponding solution obtained for both B-PINNs, with the first two moments and approximation error shown in Figure 4.9, shows high accuracy regarding the mean prediction. Although the uncertainty estimates for the surrogate solutions by both B-PINNs somewhat differ for both noise levels, the error estimates for the surrogate solutions are much smaller than the error estimates for the surrogate solutions.
levels, their magnitudes are comparable over the domain for both noise levels. Furthermore, when \(\sigma_{u}=0.01\), both methods identify some similar bulk regions with higher uncertainty, particularly near \((x,y)=(0.5,0.25)\) and \((0.75,0.75)\). Similarly for \(\sigma_{u}=0.1\), the regions of higher uncertainty show rough similarities between the two methods, with larger uncertainties clustering near the right diagonal of the domain.

Table 13 compares the mean relative errors of physical parameters \(\mathbf{x}_{i}\) and forward solution \(u\) for two noise levels, as well as the walltime and number of EKI iterations. Consistent with the previous examples, EKI-based inference is approximately 25-fold faster than that of the HMC.

## 5 Summary

This paper presents a new and efficient inference method for B-PINNs, utilizing Ensemble Kalman Inversion (EKI) to infer physical and neural network parameters. We demonstrate the applicability and performance of our proposed approach using several benchmark problems. Interestingly, while EKI methodology theoretically only provides unbiased inference with Gaussian priors and linear operators, our results show that it can still provide reasonable inference in non-linear and non-Gaussian settings, which is consistent with findings in other literature [8, 19]. In all examples, our proposed method delivers comparable inference results to HMC-based B-PINNs, but with around 8-30 times speedup. Furthermore, EKI B-PINNs can provide informative uncertainty estimates for physical model parameters and forward predictions at different noise levels comparable to the results of HMC B-PINNs. In cases where more detailed uncertainty quantification is necessary, our proposed approach can serve as a good initialization for other inference algorithms, such as HMC, to produce better results with reduced computational costs.

Besides, it is worth noting that our study did not investigate the case of a large dataset or residual points. In such cases, naive approaches to EKI would be computationally expensive due to the cubic scaling of EKI with the size of the dataset. In such a case, mini-batching techniques proposed in [32] for EKI to train NNs with large datasets, may help to overcome this challenge. Finally, we acknowledge that the EKI requires storing an ensemble of \(J\) neural network parameter sets, which can be memory-demanding for large ensemble sizes and large networks. In this situation, Dimension reduction techniques may help address this issue. We plan to investigate this strategy in future work.

\begin{table}
\begin{tabular}{c|c|c c c c c c|c|c} \hline  & & \(e_{u}\) & \(e_{x_{1}}\) & \(e_{y_{1}}\) & \(e_{x_{2}}\) & \(e_{y_{2}}\) & \(e_{x_{3}}\) & \(e_{y_{2}}\) & Walltime & Iterations \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.01\)} & EKI & 0.72\% & 0.28\% & 2.88\% & 0.60\% & 0.11\% & 0.73\% & 1.46\% & 2.66 seconds & 59 \\  & HMC & 0.72\% & 0.28\% & 2.88\% & 0.60\% & 0.11\% & 0.73\% & 1.46\% & 68.50 seconds & - \\ \hline \multirow{2}{*}{\(\sigma_{u}=0.1\)} & EKI & 1.30\% & 0.08\% & 1.02\% & 0.53\% & 0.30\% & 2.11\% & 1.24\% & 2.97 seconds & 67 \\  & HMC & 1.43\% & 0.45\% & 0.17\% & 0.23\% & 0.06\% & 5.08\% & 2.62\% & 68.58 seconds & - \\ \hline \end{tabular}
\end{table}
Table 13: Example 4.6: relative errors \(e_{u}\) of the forward solution \(u\) and \(e_{x_{i}}\), \(e_{y_{i}}\) for source locations \(x_{i}\) and \(y_{i}\) for \(i=1,2,3\) respectively for the noise levels \(\sigma_{u}=0.01,0.1\), as well as average walltime. Average EKI iterations for the EKI B-PINN are also reported.