Adaptive operator learning for infinite-dimensional Bayesian inverse problems

Consider a steady physical system described by the following PDEs:

where $\mathcal{H}$ denotes the general partial differential operator defined in the domain $\Omega\in\mathbb{R}^{d}$, $\mathcal{B}$ is the boundary operator on the boundary $\partial\Omega$, $m\in\mathcal{M}$ represents the unknown parameter field and $u\in\mathcal{U}$ represents the state field of the system.

Let $y\in\mathbb{R}^{N_{y}}$ denotes a set of discrete and noisy observations at specific locations in $\Omega$. Suppose the state $u$ and $y$ are connected through an observation system $\mathcal{O}:\mathcal{U}\rightarrow\mathbb{R}^{N_{y}}$,

where $\mathbf{\eta}\sim\mathcal{N}(0,\Sigma_{\eta})$ is a Gaussian with mean zero and covariance matrix $\Sigma_{\eta}$, which models the noise in the observations. Combining the PDE model (1) and the observation system (2) defines the parameter-to-observation map $\mathcal{G}=\mathcal{O}\circ\mathcal{F}:\mathcal{M}\rightarrow\mathbb{R}^{N_{y}}$, i.e,

Here $\mathcal{F}:\mathcal{M}\rightarrow\mathcal{U}$ is the solution operator, or the parameter-to-state map, of the PDE model (1).

The following least squared functional plays an important role in such inverse problems:

where $\|\cdot\|_{\Sigma_{\eta}}=\|\Sigma_{\eta}^{-\frac{1}{2}}\cdot\|$ denotes the weighted Euclidean norm in $\mathbb{R}^{N_{y}}$. In cases where the inverse problem is ill-posed, optimizing $\Phi$ in $\mathcal{M}$ is not a well-behaved problem, and some type of regularization is necessary. Bayesian inference is another method to consider. In the Bayesian framework, $(m,y)$ is viewed as a jointly varying random variable in $\mathcal{M}\times\mathbb{R}^{N_{y}}$. Given the prior $\nu_{0}$ on $m$, the solution to the inverse problem is to determine the distribution of $m$ conditioned on the data $y$, i.e., the posterior $\nu$ given by an infinite dimensional version of Bayes’ formula as

where $Z(y)$ is the model evidence defined as

In general, the main challenge of infinite-dimensional Bayesian inverse problems lies in well-posedness of the problem and numerical methodologies. To guarantee the well-posedness, the prior is frequently considered to be a Gaussian random field, which guarantees the existence of the posterior distribution[4, 39]. To obtain the finite posterior distributions, one can use Karhunen-Loeve (KL) expansions or direct spatial discretization methods. The posterior distribution can then be approximated using numerical techniques like Markov Chain Monte Carlo (MCMC)[40] and variational inference (VI)[41]. It should be emphasized that each likelihood function evaluation requires a forward model $\mathcal{G}$ (or $\mathcal{F}$) evaluation. The computation of the forward model can be very complicated and expensive in some real-world scenarios, making the computation challenging. As a result, it is critical to replace the forward model with a low-cost surrogate model. In this paper, we apply deep operator learning to construct the surrogate in order to substantially reduce computational time.

In this section, we employ the neural operator DeepOnet as the surrogate, which is fast to evaluate and can speed up in the posterior evaluations. The basic idea is to approximate the true forward operator $\mathcal{F}$ with a neural network $\mathcal{F}_{\theta}:\mathcal{M}\rightarrow\mathcal{U}$, where $\mathcal{M},\mathcal{U}$ are spaces defined before and $\theta$ are the parameters in the neural network. This neural operator can be interpreted as a combination of encoder $\mathcal{E}$, approximator $\mathcal{A}$ and reconstructor $\mathcal{R}$ [42] as depicted in Figs.1 and 2, i.e.,

Here, the encoder $\mathcal{E}$ maps $m$ into discrete values $\{m(\mathbf{x}_{i})\}_{i=1}^{N_{m}}$ in $\mathbb{R}^{N_{m}}$ at a fixed set of sensors $\{\mathbf{x}_{i}\}_{i=1}^{N_{m}}\in\Omega$, i.e.,

The encoded data is then approximated by the approximator $\mathcal{A}:\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}^{p}$ through a deep neural network. Given the encoder and approximator, we can define the branch net $\mathbf{\beta}:\mathcal{M}\rightarrow\mathbb{R}^{p}$ as the composition $\mathbf{\beta}(m)=\mathcal{A}\circ\mathcal{E}(m)$. The decoder $\mathcal{R}:\mathbb{R}^{p}\rightarrow\mathcal{U}$ maps the results $\{\beta_{i}\}_{i=1}^{p}$ to $\mathcal{U}$ with the form

where $t_{i}$ are the outputs of the trunk net as depicted in Fig.2. Note that by this formula, the trunk net approximates the basis of the solution space and the branch net approximates the coefficients of the expansion, which together approximate the spectral expansion of the solution space.

Combined with the branch net and trunk net, the operator network approximation $\mathcal{F}_{\theta}(m)(\mathbf{x})$ is obtained by finding the optimal $\theta$, which minimizes the following loss function:

where $\Theta$ is the parameter space. It should be noted that the loss function cannot be computed exactly and is usually approximated by Monte Carlo simulation by sampling the space $\mathcal{M}$ and the input sample space $\Omega$. That is, we take $N_{prior}$ i.i.d samples $m_{1},m_{2},\cdots,m_{N_{prior}}\sim\nu_{0}$ at $N_{\mathbf{x}}$ points $\mathbf{x}_{j}^{1},\cdots,\mathbf{x}_{j}^{N_{\mathbf{x}}}$, leading to the following empirical loss

After the operator network has been trained, an approximation of the forward model $\mathcal{G}$ can be constructed by adding the observation operator $\mathcal{O}$, i.e., $\widehat{\mathcal{G}}=\mathcal{O}\circ\mathcal{F}_{\theta}$. We then can obtain the surrogate posterior

where $\nu_{0}$ is again the prior of $m$ and $\widehat{\Phi}(m;y)$ is the approximate least-squares data misfit defined as

The main advantage of the surrogate method is that once an accurate approximation is obtained, it can be evaluated many times without resorting to additional simulations of the full-order forward model. However, using approximate models directly may introduce a discrepancy or modeling error, exacerbating an already ill-posed problem and leading to a worse solution[13]. Specifically, we can define an $\epsilon$-feasible set $\mathcal{M}(\epsilon):=\{m\in\mathcal{M}|\|\mathcal{G}(m)-\widehat{\mathcal{G}}(m)\|\leq\epsilon\},$ and the associated posterior measure $\nu(\mathcal{M}(\epsilon))$ as

Then, the complement of the $\epsilon$-feasible set is given by $\mathcal{M}^{\bot}(\epsilon)=\mathcal{M}\setminus\mathcal{M}(\epsilon)$, which has posterior measure $\nu(\mathcal{M}^{\bot}(\epsilon))=1-\nu(\mathcal{M}(\epsilon))$. We can obtain an error bound between $\widehat{\nu}$ and $\nu$ in the Kullback-Leibler distance:

It is important to note that in order for the approximate posterior $\widehat{\nu}$ to converge to the exact posterior $\nu$, the posterior measure $\nu(\mathcal{M}^{\bot}(\epsilon))$ must tend to zero. One way to achieve this goal is to enable the surrogate model $\widehat{\mathcal{G}}$ trained sufficiently in the entire input space such that the model error is small enough. However, a significant amount of data and training time are frequently required to effectively train the surrogate model. Indeed, the surrogate model only needs to be accurate within the posterior distribution space, not the entire prior space[32, 13, 43]. To maintain accurate results while lowering the computational costs, an adaptive algorithm should be developed. In the following section, we will look at how to design a framework for adaptively reducing the surrogate’s modeling error.

Developing stable and reliable adaptive surrogate modeling methods presents several challenges, especially when dealing with infinite-dimensional Bayesian inverse problems. One major challenge lies in the need to maintain the accuracy of the surrogate model in the high-density regions of the posterior distribution, where the true posterior is most concentrated. According to Theorem 1, the approximate posterior will be close to the true posterior when the surrogate is accurate in the high density region of the posterior distribution. On the other hand, the accuracy of the surrogate model, especially in the context of an operator network like DeepOnet, depends heavily on the quality and distribution of the training set. Ideally, this training set should be well-distributed across the posterior distribution in order to capture the key characteristics of the problem. However, the high density region of the posterior distribution is unknown until the observations are given. One feasible approach, as discussed earlier, is the CES framework [34], which first involves a calibration step where posterior samples are obtained using the full-order model, and then these samples can be used to train the surrogate. Nonetheless, this strategy can be very computationally expensive, especially when new observational data emerges, as the entire process must be repeated. A natural question is how to create an effective local approximation that maintains a good balance between accuracy and computational cost. To achieve this, the new adaptive framework should be designed to automatically select new training points during the posterior exploration process. These selected points are then used to fine-tune the emulator generated by DeepOnet. Instead of relying on the full-order model to explore the posterior space initially, the training points for the emulator are sequentially updated by integrating samplers that utilize the emulator itself. Based on this purpose, we separate the adaptive algorithms into offline and online stages. Offline, a small amount of samples from the prior distribution are utilized to train the DeepOnet in an acceptable length of time. The purpose of offline computation is to obtain a rough but comprehensive approximation of the DeepOnet model. The online stage consists of two major steps. First, decide whether to adaptively update the surrogate model with the posterior computation algorithm. Second, if the surrogate model needs to be updated, choose a new training set locally. Our new approach is an online strategy that depends on a specific realization of the observed data and an associated posterior approximation, in contrast to the direct DeepOnet strategy which maintains the surrogate model $\widehat{\mathcal{G}}$ unchanged. We expect that the inversion results computed by adaptive DeepOnet will have better accuracy than those produced by the offline direct DeepOnet strategies because using the data centers attention on regions of high posterior probability and causes a localization effect in the construction of surrogate. We will explore this conjecture in numerical results below.

Our approach, as depicted in Fig.3 is indeed the summary of the previous efforts. The procedure is broken down into the following steps, which are modular in nature and can be approached in various ways:

• •

  Initialization (offline): Build a surrogate $\mathcal{F}_{\theta}$ using the initial training dataset $\mathcal{D}$ with a relatively small sample size. This model is used as the initial pre-trained model.

• •

  Posterior computation: Using some numerical techniques to approximate, or draw samples from the approximate posterior $\widehat{\nu}$ induced by $\widehat{\mathcal{G}}=\mathcal{O}\circ\mathcal{F}_{\theta}$.

• •

  Refinement (online): Choose a criteria to determine whether the refinement is needed. If refinement is needed, then selecting new training points from $\widehat{\nu}$ to refine the training dataset $\mathcal{D}$ and the surrogate $\mathcal{F}_{\theta}$.

• •

  Repeat the above procedure several times until the stop criteria is met.

The rest of this part turns this outline into a workable algorithm by explaining when and where to refine, as well as how to choose new training data to improve approximations. To create an initialization surrogate, we can start by creating a small training dataset $\mathcal{D}$ from the prior distribution. In detail, suppose we have a collection of model evaluations $\mathcal{D}=\{(m_{i},\mathcal{F}(m_{i}))\}$. We then can use these points to train an operator network $\mathcal{F}_{\theta}$ and obtain a surrogate $\widehat{\mathcal{G}}=\mathcal{O}\circ\mathcal{F}_{\theta}$, which can be evaluated repeatedly at negligible cost. Subsequently, we decide where and when to refine the surrogate by conducting exploration using the current surrogate. Specifically, assume that $\widehat{\nu}_{t}$ is the approximate posterior at $t$ step. To assess the accuracy of the surrogate model in the approximate posterior space, we define the following “local” model error assist with $\widehat{\nu}_{t}$:

A natural approach is to use this error as an indicator: if the error exceeds a predefined tolerance, the surrogate model should be refined. Unfortunately, the need for high-dimensional integration makes directly evaluating this error in high-dimensional spaces prohibitively expensive. Notice that our primary goal is to maintain the surrogate accurate along the posterior computation trajectories. To achieve this, we can now obtain $T$ samples $\mathbf{M}^{(t)}:=\{m_{k}^{(t)}\}_{k=1}^{T}$ from $\widehat{\nu}_{t}$ using a posterior computation method, such as particle methods or MCMC-based sampling methods. Importantly, this process only relies on information from the surrogate model $\mathcal{\widehat{G}}_{t}$ and the data $y$, without requiring any information from the full-order model $\mathcal{G}$. To prevent the current samples from deviating too far from the true posterior trajectory, we define an anchor point $r^{t}$ from the obtained samples $\mathbf{M}^{(t)}$ using the full-order model $\mathcal{G}$ as follows:

This allows us to define an error indicator $e_{\mathcal{D}}(t)$ based on the data-fitting term at the anchor point $r^{t}$ within the posterior sample set:

This error indicator not only measures the model error but also ensures that the samples do not deviate significantly from the true posterior trajectory. When the relative error $|\frac{e_{\mathcal{D}}(t)-e_{\mathcal{D}}(t-1)}{e_{\mathcal{D}}(t)}|$ exceeds a tolerance of $\epsilon$, the surrogate must be refined in $\widehat{\nu}_{t}$. Otherwise, the refinement process stops.

The following task is to develop sampling strategies for surrogate refinement. We only need to ensure that the surrogate is accurate along the posterior evaluation trajectories, so we can generate a small set of important (adaptive) points near the anchor point and then refine the surrogate. This ensures that the refined surrogate $\widehat{\mathcal{G}}_{t}$ is locally accurate in $\widehat{\nu}_{t}$ while also reducing computational costs. This topic is relevant to optimal experimental design. A number of existing strategies can be used to accomplish this goal. In this work, we propose a greedy algorithm for efficiently selecting the most important samples from $\widehat{\nu}_{t}$. Specifically, we first draw a large set of $K$ samples $\Gamma=\{m_{1},m_{2},\cdots,m_{K}\}$ from $\widehat{\nu}_{t}$. A greedy algorithm is then used to select a subset $\gamma^{Q}=\{\widetilde{m}_{1},\widetilde{m}_{2},\cdots,\widetilde{m}_{Q}\}\subset\Gamma$ of “important” points one by one from $\Gamma$. Assuming the current selected point sets are $\gamma^{j}(j<Q)$, then the newly selected point $\widetilde{m}_{j+1}$ must be near the anchor point $r^{t}$. Furthermore, the surrogate solution for this point should have the greatest distance with the set $\widehat{\mathcal{G}}_{t}^{j}:=\{\widehat{\mathcal{G}}_{t}(\widetilde{m}_{i})\}_{i=1}^{j}$, i.e.,

where $d\big{(}\widehat{\mathcal{G}}_{t}(\cdot),\widehat{\mathcal{G}}_{t}^{j}\big{)}=\max_{\tilde{m}\in\gamma^{j}}\left\|\widehat{\mathcal{G}}_{t}(\cdot)-\hat{\mathcal{G}}_{t}(\tilde{m})\right\|_{2}$ is the distance between the value $\widehat{\mathcal{G}}_{t}(\cdot)$ and the set $\widehat{\mathcal{G}}_{t}^{j}$. Here $\lambda$ is the factor used to control the balance between these two distances and is chosen to be 1 for convenience. The process only calculates the distance between points and predicted values of the surrogate model on the preselected data set, resulting in a negligible calculation time. Then, we have $\gamma^{j+1}=\gamma^{j}\cup\{\widetilde{m}_{j+1}\}$. We can improve the DeepOnet model $\mathcal{F}_{\theta}$ by collecting the new training set $\gamma^{Q}$. Specifically, during the online training phase, we initialize the neural network parameter from the previously trained model $\mathcal{F}_{\theta}$. We expect a significant speedup in solving Eq. (5) with this initialization, which can be viewed as an example of transfer learning.

The advantages of this method are obvious. Specifically, the selected adaptive points will have varying features in the surrogate solution spaces and will be close to the anchor point $r^{t}$. This is more of an adversarial process; in order to guarantee local accuracy, we want the new points to be close to the $r^{t}$. In order to enhance generalization capacity, we also want them to incorporate as much data as the current surrogate did. The only remaining question is how to approximate the approximate posterior distribution $\widehat{\nu}_{t}$. To this end, traditional MCMC-based sampling methods and particle-based approaches[44, 45, 36, 46, 35] can be applied. In this paper, we have chosen to exclusively focus on the Unscented Kalman Inversion (UKI) method [38] in order to conceptually validate that our proposed adaptive operator learning framework. The details of the UKI algorithm will be presented in the following section. It is important to emphasize that UKI’s strategy relies on Gaussian approximations. While UKI may not provide highly accurate posterior approximations in many cases, such as non-Gaussian posteriors, we have still chosen to use UKI for sampling due to its computational efficiency. As a gradient-free, particle-based method, UKI requires only $2N_{m}+1$ full-order model evaluations per iteration and typically converges within $\mathcal{O}(10)$ iterations, making it computationally less expensive compared to MCMC-type methods. The numerical examples will demonstrate that even with the lower computational complexity of UKI, our newly designed framework still achieves significant improvements in computational efficiency. Additionally, our algorithm is easily integrated with MCMC or other particle-based methods[47, 48, 49, 50, 36, 37, 51, 52, 53], and when used with other samplers that require a larger number of samples, our framework yields even greater computational efficiency gains.

In this section, we give a brief review of the UKI algorithm discussed in [38]. The UKI is derived within the Bayesian framework and is considered to approximate the posterior distribution using Gaussian approximations on the random variable $m|y$ via its ensemble properties. We consider the following stochastic dynamical system:

where $m_{n+1}$ is the unknown discrete parameter vector, and $y_{n+1}$ is the observation vector, the artificial evolution error $\omega_{n+1}$ and observation error $\eta_{n+1}$ are mutually independent, zero-mean Gaussian sequences with covariances $\Sigma_{\omega}$ and $\Sigma_{\eta}$, respectively. Here $\alpha\in(0,1]$ is the regularization parameter, $r_{0}$ is the initial arbitrary vector.

Let $Y_{n}:=\{y_{1},y_{2},\cdots,y_{n}\}$ denote the observation set at time $n$. In order to approximate the conditional distribution $\nu_{n}$ of $m_{n}|Y_{n}$, the iterative algorithm starts from the prior $\nu_{0}$ and updates $\nu_{n}$ through the prediction and analysis steps: $\nu_{n}\rightarrow\widetilde{\nu}_{n+1}$, and then $\widetilde{\nu}_{n+1}\rightarrow\nu_{n+1}$, where $\widetilde{\nu}_{n+1}$ is the distribution of $m_{n+1}|Y_{n}$. In the prediction step, we assume that $\nu_{n}=\mathcal{N}(r_{n},C_{n})$, then under Eq. (12), $\widetilde{\nu}_{n+1}$ is also Gaussian with mean and covariance:

In the analysis step, we assume that joint distribution of $\{m_{n+1},y_{n+1}\}|Y_{n}$ can be approximated by a Gaussian distribution

where

Conditioning the Gaussian in Eq.(14) to find $m_{n+1}|\{Y_{n},y_{n+1}\}=m_{n+1}|Y_{n+1}$ gives the following expressions for the mean $r_{n+1}$ and covariance $C_{n+1}$ of the approximation to $\nu_{n+1}$:

By assuming all observations are identical to $y$ (i.e., $Y_{n}=y$), Eqs.(13)-(16) define a conceptual algorithm for using Gaussian approximation to solve BIPs. To evaluate integrals appearing in Eq. (13), UKI employs the unscented transform described below.

We obtain the following UKI algorithm in Algorithm 1 by applying the aforementioned quadrature rules. UKI is a derivative-free algorithm that applies a Gaussian approximation theorem iteratively to transport a set of particles to approximate given distributions. As a result, it only needs $2N_{m}+1$ model evaluations per iteration, making it simple to implement and inexpensive to compute. However, for highly-nonlinear problems, UKI may encounter the intrinsic difficulty of using Gaussian distributions to approximate the posterior and different samplers can be applied.

The overview of our UKI-based adaptive operator learning strategy is provided by Algorithm 2. To summarize, we begin with an offline pre-trained DeepOnet and refine the surrogate until the stop criteria is satisfied using local training data from the current approximate posterior distribution obtained by UKI. In other words, the UKI provides the DeepOnet with valuable candidate training points to refine, while the DeepOnet efficiently returns approximate data misfit information for the UKI to explore the parameter space further. Together, they establish a mutual learning system that enables them to grow and learn from one another over time. In particular, in order to approximate the induced posterior distribution $\widehat{\nu}_{t}$ and produce a series of intermediate Gaussian approximations $\mathcal{N}(r_{t+1}^{n},C_{t+1}^{n}),\,\,n=1,\cdots,T$, we first perform UKI with $T$ steps during the refinement process. The least squared error $e_{t+1}^{n}:=\Phi(r_{t+1}^{n};y)$ with $r_{t+1}^{n}$ will be calculated using Eq.(3) to form $U_{t+1}=\{(r_{t+1}^{n},C_{t+1}^{n},e_{t+1}^{n})\}_{n=1}^{T}$. It should be noted that not every step in this process is valid because the model error will eventually blow up during the iteration process. As a result, if the refinement process is terminated, we will select the last valid inversion result from this set as:

If the surrogate requires more refinement, the pair $(r_{t+1},C_{t+1})$ is used to select new training points and serve as the initial vector for the subsequent UKI iteration.

We review the computational efficiency of our method. Since the pre-trained operator network can be applied as surrogates for a class of various BIPs with models that are governed by the same PDEs but have various types of observations and noise modes. Thus, for a given inversion task, the main computational cost centers on the online forward evaluations and the online fine-tuning. On the other hand, the online retraining only takes a few seconds each time, so it can be ignored in comparison to the forward evaluation time. In these situations, the forward evaluations account for the majority of the computational cost. For our algorithm, the maximum number of online high-fidelity model evaluations is $N_{DeepOnet}=(Q+T)I_{max}$, where $Q$ is the number of adaptive samples for refinement, $T$ is the maximum number of iterations for UKI using our approach, and $I_{max}$ is the number of adaptive refinement. While $N_{\tiny{FEM}}=(2N_{m}+1)T_{\tiny{FEM}}$ represents the total evaluations for UKI using the FEM solver, $N_{m}$ represents the discrete dimension of the parameter field, and $T_{FEM}$ represents the maximum iteration number for UKI. Consequently, the asymptotic speeds increase can be computed as

Note that the efficiency of our method basically depends on $Q,I_{max}$ since $T$ is typically small. First of all, the number of adaptive samples $Q$ will be sufficiently small (e.g. $\mathcal{O}(10)$) compared to the discrete dimension $N_{m}$ (e.g. $\mathcal{O}(10^{2})-\mathcal{O}(10^{3})$), resulting in a significant reduction in computational cost as it determines the number of forward evaluations. Second, the total number of adaptive retraining $I_{max}$ is determined by the inverse tasks, which are further subdivided into in-distribution data (IDD) and out-of-distribution (OOD) cases. The IDD typically refers to ground truth that is located in the high density region of the prior distribution. Alternatively, the OOD refers to the ground truth that is located far from the high density area of the prior. The original pre-trained surrogate can be accurate in nearly all of the high probability areas for the IDD case. As a result, our framework will converge quicker. In the case of OOD, our framework requires a significantly higher number of retraining cycles in order to reach the high density area of the posterior distribution. However, $I_{max}$ for both inversion tasks can be small. Consequently, our method can simultaneously balance accuracy and efficiency and has the potential to be applied to dynamical inversion tasks. In other words, once the initial surrogate is trained, we can use our adaptive framework to modify the estimate at a much lower computational cost.

It is important to note that the UKI’s ensemble properties are used to approximate the posterior distribution using Gaussian approximations. Specifically, the sequence in Eq.(16) obtained by the full-order model $\mathcal{G}$ will converge to the equilibrium points of the following equations under certain mild conditions in the linear case [38]:

We can actually demonstrate that in the linear case, the mean vector and covariance matrix obtained by our approach will be close to those obtained by the true forward if the surrogate $\widehat{\mathcal{G}}$ is close to the true forward $\mathcal{G}$. Consider the following: $\mathop{\operator@font Range}\nolimits(\widehat{\mathcal{G}})=\mathbb{R}^{N_{m}}$, and $\widehat{\mathcal{G}}$ is linear. Using $\widehat{\mathcal{G}}$ as a surrogate, the corresponding sequence of $r_{n},C_{n}$ in Eq.(16) then converges to the following equations

In the following, we will demonstrate that if the surrogate $\widehat{\mathcal{G}}$ is near the true forward model $\mathcal{G}$, then the $\widehat{r}_{\infty},\widehat{C}_{\infty}$ ought to be near the true ones as well. We shall need the following assumptions.

We can obtain the following lemma.

Note that these assumptions are reasonable and can be found in many references [12, 32]. We will then supply the main theorem based on these assumptions.

In the first example, we consider the following Darcy flow problem:

Here, the source function $f(\mathbf{x})$ is defined as