Quantifying total uncertainty in physics-informed neural networks for solving forward and inverse stochastic problems

Introduction

How to make the best use of existing data while exploiting the information from classical mathematical models or even empirical correlations developed within a discipline is an important issue as data-driven modeling is emerging as a powerful paradigm for physical and biological systems. For example, in geophysics, researchers have been using the remote sensing data collected from multi-spectral satellites and the top-of-atmospheric reflectance model as a calibration of the data to study the soil salinization , or estimating the Earth heat loss based on the heat flow measurements and a model of the hydrothermal circulation in the oceanic crust . Data can be used to provide closures in nonlinear models or to estimate parameters or functions in mathematical models. Moreover, mathematical models can be used as additional knowledge to formulate “informative priors” in statistical estimation methods or be encoded in specially designed machine learning tools so that a smaller amount of data is required for inference of system identification. There has been recent progress for both forward (inference) and inverse (identification) problems using different methods. For example, for the forward problem, some of the popular choices of machine learning tools are Gaussian process and deep neural networks (DNNs) . For inverse problems, similar methods have been advanced, e.g., Bayesian estimation and variational Bayes inference , and have been proposed for a wide variety of objectives, from parameter estimation to discovering partial differential equations to learning constitutive relationships .

In this work we focus on the DNNs, and in particular the physics-informed neural networks (PINNs) for forward and inverse problems, first introduced in . However, in those works the mathematical models were deterministic differential equations, so here we consider stochastic differential equations that model either random micro-structure in a medium or the lack of complete knowledge (“uncertainty”), e.g., of the material property. There have only been very few works published on solving stochastic differential equations using DNNs, e.g., , for forward problems. Here we study the special and perhaps the most complex case where some of the physics is known, namely via the stochastic differential equations, and the parameter in the equation is represented as a stochastic process, introducing parametric uncertainty. First, we solve the forward stochastic Poisson equation, where there is uncertainty associated with the driving force. Subsequently, we consider the inverse stochastic elliptic equation, where the diffusivity is modeled as a random process. In the latter case, we have only partial information of the diffusivity from scattered sensors but we have much more data available for the solution, and we aim to infer the stochastic processes of not only the solution but also the diffusivity, and quantify their uncertainties given the randomness in the data. An additional uncertainty is due to the DNN approximation, which we will refer to as the approximation uncertainty. Taken together, we refer to the parametric uncertainty and the approximation uncertainty as the total uncertainty. To the best of our knowledge, the current work is the first to address total uncertainty in solving stochastic forward and inverse problems using DNNs.

In particular, in this paper we combine the arbitrary polynomial chaos (aPC) with PINNs for both the forward and the inverse stochastic problems. One of the most popular methods for uncertainty quantification studies is the polynomial chaos because it has been very effective in representing correlated stochastic fields. However, aPC is more suitable for building the orthogonal basis from arbitrary random space, without the need of any assumption on the distribution of the data. Therefore, in the current work, we employ the aPC to develop a combined method that we call NN-aPC, where we use the DNNs to learn each individual mode of the aPC expansion. More importantly, after training, the proposed method can be used to predict new realizations of the solution based only on very few measurements.

Treatment of the DNN approximation uncertainty has been addressed using different methods in the past. The traditional way to estimate uncertainty in DNNs is using the Bayes’ theorem, e.g., the Bayesian neural networks (BNNs) . BNNs are standard DNNs with prior probability distributions placed over their weights, and given observed data inference is then performed on weights. Because the inference is not tractable in general, variational inference is often used to approximate the inference . However, these models have very high additional computational cost because they require more parameters for the same network size and more time for the DNN parameters to converge. Recently, Gal et al. developed a new way to quantify uncertainty in DNNs by using dropout , which is largely used as a regularization technique to address the problem of over-fitting. Gal et al. showed that a DNN with dropout is mathematically equivalent to approximating a probabilistic deep Gaussian process , no matter what network architecture and non-linearities are used. Moreover, dropout does not induce much computation overhead and thus has been used as a practical tool to obtain uncertainty estimation effectively in real applications including language modelling , computer vision and medical applications . In this paper, dropout is used to to estimate the uncertainty in approximating each aPC mode. Based on the magnitude of this uncertainty we set up an active learning strategy and deploy additional sensors to obtain more measurements of the quantity of interest (QoI), in order to improve the predictability of PINNs.

The organization of this paper is as follows. In Section 2, we set up the data-driven forward and inverse problems. In Section 3, we introduce the PINNs for solving deterministic differential equations, followed by our main algorithm, the NN-aPC, and the method of dropout for uncertainty. In Section 4, we provide a detailed study of the accuracy and performance of the NN-aPC method for solving both the forward and inverse stochastic diffusion equation and demonstrate the effectiveness of active learning via dropout-induced uncertainty. Finally, we conclude with a brief discussion in Section 5.

Problem Setup

Suppose we have a stochastic differential equation:

We consider two types of problems here: first, a forward problem, where we know exactly the distribution of $k(x;\omega)$ everywhere in the domain $\mathcal{D}$ and $u(x;\omega)$ is our QoI; and second, an inverse problem, where we assume that we have incomplete information on $k(x;\omega)$ but some extra knowledge on $u(x;\omega)$, and we are interested in inferring the full stochastic profile of $k(x;\omega)$. In practice, both problems are data-driven, since the information usually comes from data collected via sensor measurements. Here, we summarize the different scenarios of the sensors placement for each type of the problems:

Forward problem: The $u$-sensors are placed only at the boundary $\Gamma$ to provide boundary condition, while the $k$-sensors are virtual (since we know the distribution of $k$), thus we can have as many $k$-sensors as we want and they can be placed anywhere in $\mathcal{D}$.

Inverse problem: In addition to having $u$-sensors at the boundary $\Gamma$, we have a limited number of extra $u$-sensors that can be placed in the domain $\mathcal{D}$, whereas we only have a limited number of $k$-sensors.

In this paper, we address both types of problems but we will focus more on solving the inverse problem.

Methodology

In this part, we briefly review using DNNs to solve deterministic differential equations , and its generalization for solving deterministic inverse problems in . To this end, re-consider Equation 1 but replace the random input $\omega$ and approximate it with a finite set of parameters, leading to a parameterized differential equation:

where $u(x)$ is the solution and $\eta$ denotes the parameters.

A DNN, denoted by $\hat{u}(x;\theta)$, is constructed as a surrogate of the solution $u(x)$, and it takes the coordinate $x$ as the input and outputs a vector that has the same dimension as $u$. Here we use $\theta$ to denote the DNN parameters that will be tuned at the training stage, namely, $\theta$ contains all the weights $\bm{w}$ and biases $\bm{b}$ in $\hat{u}(x;\theta)$. For this surrogate network $\hat{u}$, we can take its derivatives with respect to its input by applying the chain rule for differentiating compositions of functions using the automatic differentiation, which is conveniently integrated in many machine learning packages such as Tensorflow . The restrictions on $\hat{u}$ is two-fold: first, given the set of scattered data of the $u(x)$ observations, the network should be able to reproduce the observed value, when taking the associated $x$ as input; second, $\hat{u}$ should comply with the physics imposed by Equation 2. The second part is achieved by defining a residual network:

which is computed from $\hat{u}$ straightforwardly with automatic differentiation. This residual network $\hat{f}$, also named the physics-informed neural network (PINN), shares the same parameters $\theta$ with network $\hat{u}$ and should output the constant 0 for any input $x\in\mathcal{D}$. Figure 1 shows a sketch of the PINN. At the training stage, the shared parameters $\theta$ (and also $\eta$, if it is also to be inferred) are fine-tuned to minimize a loss function that reflects the above two constraints.

Suppose we have a total number of $N_{u}$ observations on $u$, collected at location $\{x_{u}^{(i)}\}_{i=1}^{N_{u}}$, and $N_{c}$ is the number of collocation points $\{x_{f}^{(i)}\}_{i=1}^{N_{c}}$ where we evaluate the residual $\hat{f}(x_{f}^{(i)};\theta,\eta)$. We shall use $(x^{*},y^{*})$ to represent a single instance of training data, where the first entry $x^{*}$ denotes the input and the second entry $y^{*}$ denotes the anticipated output (also called “label”). The workflow of solving a differential equation with PINN can be summarized as follows:

2 NN-aPC: Combining arbitrary polynomial chaos with neural networks

We generalize the PINN method to solve stochastic differential equations for both forward and inverse problems, i.e., we aim to infer continuous random processes. Assume that a sensor will generate a sequence of measurements after being installed, and when the data is recorded, all sensors are read simultaneously. We denote the measurements from all the sensors at the same instant by a snapshot of the sensor data. Although the measurement results change from one measurement to the next due to randomness, it is reasonable to believe that every snapshot of sensor data corresponds to the same random event in the random space. We also assume that when the number of snapshots is big enough, the empirical distribution approximates the true distribution.

Let us consider Equation 1. Suppose we have $N_{k}$ sensors for $k(x;\omega)$ placed at $\{x_{k}^{(i)}\}_{i=1}^{N_{k}}$, $N_{u}$ sensors for $u(x;\omega)$ placed at $\{x_{u}^{(i)}\}_{i=1}^{N_{u}}$, and $N_{f}$ collocation points at $\{x_{f}^{(i)}\}_{i=1}^{N_{f}}$ that are used to calculate the residual of Equation 1. A total number of $N$ snapshots of measurements are made from all these sensors. Let $k_{s}^{(i)}$ and $u_{s}^{(i)}$ ($s=1,2,...,N$) be the $s$-th measurement of $k$ and $u$ at location $x_{k}^{(i)}$ and $x_{u}^{(i)}$ respectively, and $\omega_{s}$ is the random instance at the $s$-th measurement, i.e., $k_{s}^{(i)}=k(x_{k}^{(i)};\omega_{s})$ and $u_{s}^{(i)}=u(x_{u}^{(i)};\omega_{s})$. The training data set can be represented by

The proposed NN-aPC method consists of the following steps:

building the NN-aPC as a surrogate model of aPC modes and train the network for each mode.

The trained NN-aPC can then be used to calculate the statistics of our QoI and to predict new instances of the continuous trajectories of the QoI, with newly collected sensor data. We will explain each of three steps and the prediction procedure below.

As the first step, we find a lower dimensional random space spanned by a set of hidden random variables for the dimension reduction of our QoI. The most convenient way to do this is via the principal component analysis (PCA). Naturally, we would analyze the data of $k$, which is the source of randomness in Equation 1. Let $K$ be the $N_{k}\times N_{k}$ covariance matrix for the sensor measurements on $k$, i.e.,

Let $\lambda_{l}$ and $\phi_{l}$ be the $l$-th largest eigenvalue and its associated normalized eigenvector of $K$. Therefore, PCA yields

where $\bm{\Phi}=[\phi_{1},\phi_{2},...,\phi_{N_{k}}]$ is an orthonormal matrix and $\bm{\Lambda}=\text{diag}(\lambda_{1},\lambda_{2},...\lambda_{N_{k}})$ is a diagonal matrix. Let $\bm{k}_{s}=[k_{s}^{(1)},k_{s}^{(2)},...,k_{s}^{(N_{k})}]^{T}$ be the results of the $k$ measurements of the $s$-th snapshot, then

is an uncorrelated random vector, and hence $\bm{k}_{s}$ can be rewritten as a reduced dimensional expansion

where $k_{0}(x_{k}^{(i)})$ is the mean of $k$ measurements at $x_{k}^{(i)}$, $k_{l}(x)$ is the $l$-th mode function of $k(x;\omega)$ whose value at $x_{k}^{(i)}$ coincides with the $i$-th entry of the eigenvector $\phi_{l}$, and $\xi_{s,l}$ is the $l$-th entry of the random vector $\bm{\xi}_{s}$. We want to extend the range of $k_{l}$ to the entire domain to approximate the continuous samples of $k(x;\omega)$.

2.2 Arbitrary polynomial chaos

Assume that we have a set of $M$-dimensional samples of random vectors

with hidden probability measure $\rho(\bm{\xi})$. Given a sufficiently large number of snapshots, we can approximate the underlying probability measure $\rho(\bm{\xi})$ by the discrete measure $\nu_{S}(\bm{\xi})$,

Specifically, the basis $\{\psi_{\alpha}(\bm{\xi})\}_{\alpha=0}^{P}$ are constructed using the recursive algorithm

where $\psi^{*}_{\alpha}(\bm{\xi}):=\prod_{i=1}^{M}\bm{\xi}_{i}^{\alpha_{i}}$ represents the multivariate monomial basis function; the coefficients $w_{\beta}^{\alpha}$ are determined by imposing the orthonormal condition with respect to the discrete measure $\nu_{S}$, i.e.,

With the polynomial basis $\{\psi_{\alpha}(\bm{\xi})\}$ that are automatically adapted to the distribution of $\bm{\xi}$, we can write any function $g(x;\bm{\xi})$ in the form of the aPC expansion,

where the functions $g_{\alpha}(x)$ are called the aPC modes of $g$ and can be calculated by

2.3 Learning stochastic modes

The key to our method is to train DNNs that predict the stochastic modes of our QoI. In this section, we focus on the inverse problem where we have to learn both the modes of $u$ and $k$. (Solving a forward problem is similar and more straightforward, and will be briefly discussed in Section 4.1.1.) Two disjoint DNNs are constructed, i.e., the network $\widehat{u_{\alpha}}$, which takes the coordinate $x$ as the input and outputs a $(P+1)\times 1$ vector of the aPC modes of $u$ evaluated at $x$, and the network $\widehat{k_{i}}$ that also takes the coordinate $x$ as the input and outputs a $(M+1)\times 1$ vector of the $k$ modes (we take $k_{0}$ in Equation 11 as the 0-th mode of $k$). Then, we can approximate $k$ and $u$ at the $s$-th snapshot by

In practice, for $\widehat{k_{i}}$, we separate the mean from the rest of the modes and learn it with a small scale DNN, and for $\widehat{u_{\alpha}}$, we group the modes corresponding to the same order of aPC expansion together and learn each group of modes with a separate DNN, as depicted in Figure 3. This is due to fact that the mean and the modes of different orders often correspond to vastly different scales.

The loss function is defined as a sum of the mean squared errors (MSE):

So far, we have specified the training data, constructed the DNNs and formalized the loss function, and we now ready to train the DNNs.

2.4 Predicting stochastic realizations

In real applications, the training set could be generated from historical data, and the training process should be performed at the offline stage. The fine-tuned model shall be used to predict new random instances provided with new snapshots of sensor data, at the online stage. Suppose we have a snapshot of sensor data $\{\{(x_{k}^{(i)},k_{\text{new}}^{(i)})\}_{i=1}^{N_{k}},\{(x_{u}^{(i)},u_{\text{new}}^{(i)})\}_{i=1}^{N_{u}}\}$; the first step is to extract the hidden random variables $\bm{\xi}_{\text{new}}$ from $\{k_{\text{new}}^{(i)}\}_{i=1}^{N_{k}}$ using Equation 9. Then, for any assigned location $x$, we can predict the modal functions for both $k(x;\omega)$ and $u(x;\omega)$ from the trained DNNs $\widehat{k_{i}}$ and $\widehat{u_{\alpha}}$. Finally, the prediction of $k$ and $u$ for the new random instance can be made via Equation 18 and 19, respectively.

3 Dropout for uncertainty

Although DNNs can be used to approximate any measurable function accurately, standard DNNs do not capture model uncertainty. Dropout is one convenient way to quantify the approximation uncertainty in DNNs. The key idea of dropout is to drop units from the DNN independently and randomly with a pre-selected probability $p\in(0,1)$. In the original work, dropout was only used during training, while no units were dropped at test time, i.e., the prediction of the DNN for the unknown data was deterministic. In the dropout for uncertainty, the units are also dropped at test time, resulting in stochastic predictions each time. The loss in the dropout inference is the summation of the original loss and a $l_{2}$ regularization term over the DNN parameters:

where $N_{t}$ is the number of training points, $\widehat{y_{i}}$ is the prediction, $y_{i}$ is the true value, $l(\cdot,\cdot)$ is the loss for a single prediction, $\theta_{i}$ is any weight and bias, and $\lambda$ is the $l_{2}$ regularization rate. During prediction, the mean of the output is directly estimated by the Monte Carlo (MC) method,

where $\mathcal{NN}_{t}$ is the dropped neural network at the $t$-th prediction. The output variance is also estimated from these MC outputs.

Figure 4 shows an example of using the dropout DNN for regression. We plan to use the dropout strategy in the NN-aPC method to estimate the uncertainty of our DNN model and as a guidance for active learning.

Numerical Examples

We first demonstrate the effectiveness of solving stochastic differential equations with the NN-aPC method for the forward and inverse problems.

Consider the following one-dimensional stochastic Poisson equation with homogeneous boundary conditions:

Here $\Omega$ is the random space, the forcing term $f(x;\omega)\sim\mathcal{GP}(f_{0}(x),\text{Cov}(x,x^{\prime}))$ is a Gaussian random process with mean $f_{0}(x)=10\sin(\pi x)$ and a squared exponential covariance function

where the standard deviation $\sigma=1.0$ and the correlation length $l_{c}=0.5$.

In this and the following examples, all data are generated by the MC sampling method. Specifically, we sample $N=1000$ snapshots of continuous $f(x;\omega)$ trajectories $\{f_{s}=f(x;\omega_{s})\}_{s=1}^{N}$ and extract from $\{f_{s}\}_{s=1}^{N}$ the values where the $N_{f}$ (virtual) sensors are located. For every $f(x;\omega)$ trajectory, we solve for its corresponding solution trajectories $u(x;\omega)$ using the finite difference method, and will use the statistics of these $u$ trajectories as our reference. To evaluate the performance of the trained model, we collect another $N_{s}=500$ snapshots of continuous $f(x;\omega)$ and $u(x;\omega)$ trajectories independently from the training data. The $N_{s}$ pairs of $(f,u)$ trajectories form our test sample set. Similarly, we extract the $f$-sensor data from every snapshot in the test set as the input at the predicting stage. We shall use the same $N$ and $N_{s}$ in the following tests, if not explicitly mentioned.

We place $N_{f}=13$ sensors of $f(x;\omega)$ in the $$domain (the sensors are equidistant) and keep$6$principal random variables corresponding to$99\%$stochastic energy after performing PCA. The solution$u(x;\omega)$is approximated with a first-order aPC expansion. Figure 5 shows the scattered plots of the measurements from the first three$f$-sensors and the first three arbitrary polynomial basis. It is evident that the raw data from the measurements are correlated while their induced polynomials are not, thus the induced polynomials would serve as a valid set of basis in the random space.

The DNNs used to approximate the modes of $u$ are constructed as in Figure 3, where we use an isolated small scale DNN of 2 hidden layers with 4 neurons per hidden layer to approximate the mean profile, and a DNN of 4 hidden layers with 32 neurons per hidden layer to model the modes. The tanh function is selected as the default activation function due to it is second order differentiable. Then, a DNN for the residual can be constructed via auto-differentiation and arithmetic operations. The training set $\mathcal{S}_{t}$ is

where the $u$ data is collected only at the boundaries to provide boundary conditions. The loss function is slightly modified based on Equation 20 to add a $l_{2}$ regularization term. At the training stage, we choose the $l_{2}$ regularization rate $\lambda=0.001$, and use the Adam optimizer with learning rate $0.001$ to train our model for $20000$ epochs. Figure 6 shows the predicted mean and standard deviation of the solution $u$ versus the reference. Figure 7 shows our DNN prediction of three $u$ modes where the reference modes are calculated by Equation 17. We can see that the NN-aPC method makes accurate predictions of the mean and standard deviation of the solution $u(x;\omega)$ and learns the arbitrary polynomial chaos modes.

The trained model is then used to predict the QoI at any location $x_{o}$ given new snapshots of sensor data, with only one forward evaluation of the DNN with $x_{o}$ as the input, as described in Section 3.2.4. Figure 8(a) illustrates the prediction of solution for three different snapshots of the sensor data in the test samples; our prediction recovers the true solution very well. In Figure 8(b), we use an increasing value of $N_{f}$ to study the effect of the number of $f$-sensors on the accuracy of the predicted solutions; we observe that when more $f$-sensors are deployed, we can achieve better accuracy.

1.2 Inverse problem: stochastic elliptic equation

We solve the one-dimensional stochastic elliptic equation as an inverse problem, where we have some extra information on the solution $u(x;\omega)$ but incomplete information of the diffusion coefficient $k(x;\omega)$. The equation reads

In this example, we use a constant forcing term $f(x)=10$. The randomness comes from the diffusion coefficient $k(x;\omega)$, for which we only have limited information at the locations where we place the $k$-sensors. Here in this example, $k$ is modeled and sampled from a non-Gaussian random process such that

where the mean $k_{0}(x)=\sin(3\pi x/2)/5$ and the covariance function has the same form as in Equation 27, where we set the standard deviation $\sigma=0.1$ and correlation length $l_{c}=1.0$. We use the same strategy as in Section 4.1.1 to generate the training and testing samples, and the training set is constructed as in Equation 6.

For this case, two groups of DNNs, i.e., $\widehat{k_{i}}$ and $\widehat{u_{\alpha}}$ are built to calculate the modes for $k$ and $u$. Again, we use the Adam optimizer with learning rate 0.001 to train the DNNs for 50000 epochs. In Figure 9(a), we use the 1st-order aPC expansion and we study the impact of using different $l_{2}$ regularization rate $\lambda$ and different shapes of DNNs. We only change the DNNs that learn the stochastic modes, while the DNNs that learn the mean profiles are fixed to have 2 hidden layers and 4 neurons per hidden layer; this is the default setting for the future examples as well. In these plots, we compare the averaged relative $L_{2}$ error of predicting $k$ and $u$ in the test set. The results indicate that a moderate choice of $\lambda$ (0.0005) gives us the most accurate predictions, since a too small/large choice of $\lambda$ causes over-/under- fitting. A suitable choice of DNN shape (4 hidden layers with 32 neurons per hidden layer) produces the best trained model. For the rest of this numerical example, we shall adopt these optimal DNN setting.

In Figure 9(b), we use the 1st-order aPC expansion and compare the averaged relative $L_{2}$ error in predictions when different numbers of $k$- and $u$-sensors are deployed to collect the training data. In general, the proposed method makes more accurate predictions after training with data collected from more $k$- and $u$-sensors. One reason is that a larger number of sensors supports a greater variety of input $x$ in the training data, thus feeding more information to the model to reduce the probability of over-fitting. Also, more $k$-sensors allows for a higher effective random dimension of the aPC expansion for better approximation. Figure 9(b) also shows that a 2nd-order aPC expansion helps to improve predictions. We note that this is not the case when we use only three $k$-sensors. The bottleneck here is insufficient random dimension, so without enough training information, adopting the 2nd-order aPC expansion doubles the number of $\widehat{u_{\alpha}}$ net outputs and would only increase the risk of over-fitting.

We use a combination of 7 $u$-sensors and 4 $k$-sensors. Figure 10(a) and 11(a) compare the mean and standard deviation of $u$ and $k$ calculated by the trained DNNs when we use the 1st- and 2nd-order aPC expansions. Figure 10(b) shows the first four common aPC modes of $u$ for both expansions, where we can see that the 2nd-order aPC expansion helps to improve the accuracy of the predicted lower-order modes. The 2nd-order aPC expansion would also help with learning the $k$ modes, as depicted in Figure 11(b), even if the $\widehat{k_{i}}$ network is disjoint with the $\widehat{u_{\alpha}}$ network. Table 1 lists the relative $L_{2}$ error of the trained model in calculating the mean, standard deviation and all the common modes of $k$ and $u$. We can see that using a 2nd-order aPC significantly improves the accuracy. We also note that in Figure 11(a), due to the placement of $k$-sensors, the measurements from each sensor will yield almost the same mean (1.0) and standard deviation (0.1), but the trained $k$ net would reveal the non-trivial wavy structure of its mean and standard deviation in the entire domain, containing information more than just the training data could provide. This indicates that there is information fusion of all three types of training data and the stochastic differential equation, rather than a simple interpolation. Figure 12 shows that the higher-order aPC modes can also be captured even if they exist in a much smaller scale compared to the lower-order more energetic modes.

Figure 13 shows the prediction of $k$ and $u$ for three arbitrary snapshots in the test sample set. Every pair of predictions are made based on a single time measurement from 4 $k$- and 7 $u$-sensors, and they agree with the true reference value. Again, the predicting stage takes little computational time since only one forward evaluation of the well-trained DNNs is needed for every input $x$, no matter how many snapshots of predictions are to be made.

2 Solving deterministic differential equations with PINN and dropout

In this part, we first show that dropout reduces over-fitting in solving differential equations. Then, more importantly, we show that the dropout-induced uncertainty serves as useful guidance for active learning.

To show how dropout reduces over-fitting, we implement the dropout neural networks to solve both a forward Poisson equation and an inverse elliptic equation. They are the deterministic version of Equation 26 and 28. For the forward Poisson equation, we choose $f(x)=9\pi^{2}\sin(3\pi x/2)/4$ as the forcing term. A dropout neural network with 6 hidden layers and 100 neurons per hidden layer is constructed to model the solution $u(x)$. The training data consists of 2 $u$ measurements at both boundaries and 6 $f$-sensors in the domain. For the inverse elliptic equation, we choose, as before, $f(x)=10$ and $k(x)=\exp(\sin(3\pi x/2)/5)$ as the hidden diffusion coefficient, and we use 5 $k$-sensors and 7 $u$-sensors. Two separate DNNs are constructed: a small scale regular DNN that has 2 hidden layers with 4 neurons per hidden layer to model the function $u(x)$, and a large scale dropout neural network with 6 hidden layers and 100 neurons per hidden layer to model the function $k(x)$. The dropout rate in both examples is fixed at 0.01.

For both the forward and inverse problem, we train the DNNs as described in Section 3.3, using an Adam optimizer with learning rate $0.001$ for 30000 epochs. Due to the lack of sufficient number of sensors, there is a big chance that over-fitting occurs at the training stage if we use just regular DNNs. Figure 14 shows a comparison of the results when we train the networks with and without using dropout. As we can see in the plots, the results from a regular DNN are very different from each other, showing irregular jumps of large amplitudes, while the results from dropout DNNs are similar to each other, and they are closer to the truth. This shows that dropout works as an effective means of reducing over-fitting.

2.2 Estimating DNN approximation uncertainty

More importantly, the uncertainty introduced by dropout serves as a useful guidance for active learning. Again, we solve the same inverse elliptic equation but this time we are provided with additional $k$-sensors. At the training stage, we add an $l_{2}$ regularization term with $\lambda=10^{-6}$ to the loss function. At the predicting stage, we evaluate the dropout network for $10000$ times to estimate the mean and standard deviation. Next, a new $k$-sensor is placed where the standard deviation reaches its maximum. If it happens that the location to add the new sensor is close to an existing $k$-sensor by a threshold distance of $\rho$, we do not add the new sensor, but instead, we count the nearest existing sensor twice as if we have added a virtual sensor at the same location of the existing one. In practice, here we choose $\rho=0.03$. We have designed an automatic iterative procedure for active learning that reduces the cost of adding new sensors by efficiently re-using the old ones.

Figure 15 shows the initial prediction of $k$ and the prediction after iterating the above algorithm for 15 steps. In Figure 15(b), the sensors are clustered where the curvature of $k$ is big, which is consistent with our intuition that we should put more sensors where the function changes rapidly. In Figure 16, we see that during these 15 iterations, only 8 new sensors are deployed while the relative error of $k$ predictions reduces from more than $5\%$ to less than $1\%$.

3 Active learning for inverse stochastic elliptic problems

We consider the inverse stochastic elliptic problem as described in Equation 28 for active learning, but in this example, $\log(k(x;\omega))$ is modeled by a Gaussian random process with correlation length $l_{c}=0.5$. To start with, we have 1000 snapshots of data from 3 $k$-sensors, 7 $u$-sensors and 21 $f$-sensors that are equidistantly distributed in the physical domain, and our goal is to infer $k(x;\omega)$ in the entire domain. Suppose we are then provided with additional sensors of $k$; we shall allocate them according to the uncertainty induced by the dropout neural networks. In practice, we model the modes of $k(x;\omega)$ with a dropout neural network that has 4 hidden layers with 128 neurons per hidden layer, and a dropout rate 0.01. The solution $u(x;\omega)$ is expanded with the 1st-order aPC expansion, while the modes are modeled by a regular DNN with 4 layers and 32 neurons per hidden layer. The reasons for implementing dropout only on $\widehat{k_{i}}$ are as follows:

Dropout can be used to efficiently reduce over-fitting, and the over-fitting issue of $k(x,\omega)$ is worse than that of $u(x,\omega)$.

As an inverse problem, our main goal is to identify $k(x,\omega)$ and then identify where we should add more $k$-sensors to enhance the accuracy of prediction.

We use an Adam optimizer with learning rate 0.0005 to train the networks for 50000 epochs. The mean and standard deviation of the $k$ modes model are evaluated from 10000 independent evaluations of the dropout neural network. We place the additional $k$-sensor where the standard deviation of the first mode attains its maximum. This is because the first mode is associated with the largest eigenvalue, therefore bringing the largest impact to the stochastic structure, thus it is most important to learn the first mode accurately.

We carry out the active learning steps until the $k$ prediction error does not decrease after a new sensor is added. To better illustrate the process of adding new sensors, in Figure 17(a), 17(b) and 17(c) we show the learned first mode of $k$ and its associated standard deviation from dropout uncertainty. Three steps of active learning are displayed. Note that the shapes of the first modes do not stay the same due to the fact that every time a new $k$-sensor is added, the principal components of $K$ in Equation 8 are therefore changing. Figure 18(a) and Figure 18(b) show the comparison of the predicted standard deviation of $k$ and $u$ in the first three steps and the last step, respectively. It is evident that adding extra $k$-sensors automatically according to the dropout-induced uncertainty will improve the accuracy of standard deviation prediction for both $k$ and $u$. Finally, the trained model is used to predict continuous trajectories of $k$ and $u$ in the test samples. We can conclude from Table 2 that adding extra $k$-sensors based on active learning with the dropout helps us to make better predictions.

Summary

We have presented a new approach to quantify the parametric uncertainty in the physics-informed neural networks (PINNs) by employing the arbitrary polynomial chaos (aPC) expansion to represent the stochastic solution. We use the data collected from sensor measurements to build a set of arbitrary polynomial basis and learn the modal functions of the aPC expansion through the PINNs, i.e., DNNs that encode the underlying stochastic differential equation. The proposed data-driven method can be used to solve forward problems, but more importantly, it deals with stochastic inverse problems. In the classical inverse problem, typically all the information available is for the solution, and we aim to identify the parameters. Here we selected to solve the inverse problems, where we have partial information available both for the solution and the parameter, which is a stochastic process. Once the model is trained with existing sensor data, i.e., historical data, it can be used to predict new instances of trajectories of the quantity of interest (solution or parameter) at very small additional computational cost.

We aim at quantifying two different types of uncertainties, i.e., the parametric uncertainty due to the stochastic equation, as well as the approximation uncertainty of the PINN. The latter represents how well the PINN is trained and how robust it is at the predicting stage. To this end, we adopt the dropout strategy for estimating the approximation uncertainty. Dropout is typically used to deal with over-fitting problems but it can also be exploited to quantify the approximation uncertainty at no additional cost. In our examples, we use DNNs to learn the modal functions of the stochastic parametric modes and the dropout strategy to quantify their associated uncertainty. Based on this, we propose an iterative method of actively learning where to place new sensors to enhance the approximation accuracy of PINNs. The numerical results exhibit the effectiveness of such an active learning strategy, not only in placing new sensors but also in making better use of the existing ones.

There are other possible methods of quantifying parametric and approximation uncertainties. For example, we can abandon aPC and use directly the stochastic data as the input. One possible approach is to consider the random space together with the physical space. Hence, standard DNNs, and in particular the deterministic PINNs developed in can be directly used to stochastic differential equations. However, we did not choose this approach for two reasons: first, it does not lead to explicit expressions for stochasiticity of the quantity of interest (QoI); and second, different from sampling in the physical space where we can always mark the location by their coordinates, marking random instances with random variables is much harder. Moreover, the dimension of the random space is usually much higher than that of the physical space, requiring a proper dimension reduction procedure to be carried out before actually solving the problem. Although weaker, the high dimensionality is also a limitation of the aPC approach because it requires the high-dimensional DNN outputs, which could make the training process harder. Other promising methods that may be able to deal with high dimensional stochastic problems are the generative adversarial networks (GANs) , and we plan to systematically investigate this line of research in our future work.

Acknowledgement

This work is supported by DARPA N66001-15-2-4055, Air Force FA9550-17-1-0013, ARL W911NF-12-2-0023, NSF of China (No. 11671265) and the Science Challenge Project (No. TZ2018001). In addition, we would like to thank Liu Yang for his generous advice.

References

References