Physics-Informed Generative Adversarial Networks for Stochastic Differential Equations

Introduction

Generative adversarial networks (GANs) have achieved remarkable success within short time for diverse tasks of generating synthetic data, such as images , texts , and even music . GANs can learn probability distributions from data, an attribute suggestive of its potential application to modeling the inherent stochasticity and extrinsic uncertainty in physical and biological systems. However, to the best of our knowledge, there is no work explicitly encoding the known physical laws into the framework of GANs so far in the spirit of physics-informed neural networks first introduced in .

The specific data-driven approach to modeling physical and biological systems depends crucially on the amount of data available as well as on the complexity of the system itself, as illustrated in Figure 1. The classical paradigm for which many different numerical methods have been developed over the last fifty years is shown on the the top of Figure 1, where we assume that the only data available are the boundary and initial conditions while the specific governing partial differential equation (PDE) and associate parameters are precisely known. On the other extreme (lower plot), we may have a lot of data, e.g. in the form of time series, but we may not know the governing physical law, e.g. the underlying PDE, at the continuum level; many problems in social dynamics fall under this category although work so far has focused on recovering known PDEs from data only, e.g. see . Perhaps the most interesting category is sketched in the middle plot, where we assume that we know the physics partially, e.g. in an advection-diffusion-reaction system the reaction terms may be unknown, but we have several scattered measurements in addition to the boundary and initial conditions that we can use to infer the missing functional terms and other parameters in the PDE and simultaneously recover the solution. It is clear that this middle category is the most general case, and in fact it is representative of the other two categories, if the measurements are too few or too many. This is the mixed case that we address in this paper but with the significantly more complex scenario, where the solution is a stochastic process due to stochastic excitation or an uncertain material property, e.g. permeability or diffusivity in a porous medium. Hence, we employ stochastic differential equations (SDEs) to represent these stochastic solutions and other stochastic fields.

Taking inspiration from the work on physics-informed neural networks for deterministic PDEs , here we employ GANs for stochastic problems in the second category of Figure 1. We wish to encode the known physics, more specifically, the form of the stochastic differential equation (SDE), into the architecture of GANs while at the same time exploit the feature of GANs to model and learn the unknown stochastic terms in the equations from data. This approach represents a seamless integration of models and data both for inference and system identification but also for the aforementioned mixed case, where we have insufficient data for both, and hence we wish to infer both the system and the state. This is an emerging new paradigm in machine learning research addressing engineering applications, where typically the physics is very complex and we only have partial measurements of forcing or material properties. To this end, recent advances include using Gaussian process regression and deep neural networks (DNNs) to solve forward problems as well as Bayesian estimation and DNNs for inverse problems. However, published works are mostly dealing with deterministic systems. There have only been very few works published on data-driven methods for SDEs, e.g., , for forward problems; for the inverse problem, Zhang et.al. have recently proposed a DNN based method that learns the modal functions of the quantity of interest, which could also be the unknown system parameters. While robust, this method suffers from the “ curse of dimensionality” in that the number of polynomial chaos expansion terms grows exponentially as the effective dimension increases, leading to prohibitive computational costs for modeling stochastic systems in high dimensions. In this paper, we propose physics-informed GANs (PI-GANs) to solve SDEs. Our method is flexible in that without changing the general framework, it is capable of solving a wide range of problems, from forward problems to inverse problems and in between, i.e., mixed problems. Moreover, PI-GANs do not suffer from the “curse of dimensionality”. As a result, we can possibly tackle SDEs involving stochastic processes with high effective dimensions. In addition, PI-GANs can make use of data collected from multiple groups consisting of multiple sensors with no alignment between data in the groups.

The organization of this paper is as follows. In Section 2, we set up the data-driven problems. In Section 3, we give a brief review of GANs and the specific version we use, namely Wasserstein GANs with gradient penalty (WGAN-GP). Our main algorithms are introduced in Section 4, followed by a detailed study of the performance of our methods in Section 5. We first consider the learning of stochastic processes from a limited number of realizations and discuss the currently under-studied issue of overfitting in WGAN-GP. Subsequently, we present solutions of stochastic PDEs for different types of available data and different dimensions. We conclude in Section 6 with a brief summary and discussion of the current limitations of the method.

Problem setup

To illustrate the main idea of PI-GANs we consider the following stochastic differential equation:

where $\mathcal{N}_{x}$ is a general differential operator, $\mathcal{D}$ is the $d$-dimensional physical domain in $\rm I\!R^{d}$, $\Omega$ is the probability space, and $\omega$ is a random event. The coefficient $k(x;\omega)$ and the forcing term $f(x;\omega)$ are modeled as random processes, and thus the solution $u(x;\omega)$ will depend on both $k(x;\omega)$ and $f(x;\omega)$. $B_{x}$ is the boundary condition operator acting on the domain boundary $\Gamma$. As a pedagogical example, in this paper we consider the one-dimensional stochastic elliptic equation, which retains most of the main features of more complex SDEs:

where $k(x;\omega)$ and $f(x;\omega)$ are independent stochastic processes, and $k(x;\omega)$ is strictly positive. For simplicity, we impose homogeneous Dirichlet boundary conditions on $u(x;\omega)$.

We consider a general scenario for Equation (1), where we have a limited number of measurements from scattered sensors for the stochastic processes. Specifically, we place sensors at $\{x^{k}_{i}\}_{i=1}^{n_{k}}$, $\{x^{u}_{i}\}_{i=1}^{n_{u}}$, $\{x^{f}_{i}\}_{i=1}^{n_{f}}$ and $\{x^{b}_{i}\}_{i=1}^{n_{b}}$ to collect “snapshots” of $k(x;\omega)$, $u(x;\omega)$, $f(x;\omega)$ and $b(x;\omega)$, where $n_{k}$, $n_{u}$, $n_{f}$ and $n_{b}$ are the numbers of sensors for $k(x;\omega)$, $u(x;\omega)$, $f(x;\omega)$ and $b(x;\omega)$, respectively. Here, one “snapshot” represents a simultaneous read of all the sensors, and we assume that the data in one snapshot correspond to the same random event $\omega$ in the random space $\Omega$, while $\omega$ varies for different snapshots. Note that each snapshot is the concatenation of snapshots from $k(x;\omega)$, $u(x;\omega)$, $f(x;\omega)$ and $b(x;\omega)$, thus it is actually a vector of size $(n_{k}+n_{u}+n_{f}+n_{b})$. Suppose we have a group of $N$ snapshots, then we denote the accessible data set as $\{T(\omega^{(j)})\}_{j=1}^{N}$ defined as

The corresponding terms in Equation (3) are omitted if we put no sensors for that specific process.

In this paper, we assume that we always have a sufficient number of sensors for $f(x;\omega)$ in Eqn 2. The type of the problems varies depending on the number of sensors on $k(x;\omega)$ and $u(x;\omega)$: as we decrease the number of sensors on $k(x;\omega)$ while increase the number of sensors on $u(x;\omega)$, the problems gradually transform from forward to mixed and finally to inverse problems.

Moreover, we could have more than one group of snapshots of measurements, so in this case, our accessible data are

where $M$ is the number of groups, $t$ is the index for the groups, and $N_{t}$ is the number of snapshots in group $t$. The random event $\omega^{(t,j)}$ denotes the random instance of the $j$-th snapshot in group $t$. Note that the sensor setups in different groups could be different, and we use $\{x^{k,t}_{i}\}_{i=1}^{n_{k,t}}$ to denote the position setup of $n_{k,t}$ sensors for $k$ in group $t$ (similarly for other terms). More importantly, snapshots from different groups could be collected independently. For example, in addition to the existing old sensors, we may place some new sensors and collect and utilize data from both the new and the old sensors. In this setup, we can make use of both the newly collected data and the previously collected data solely from the old sensors.

A brief review of GANs and WGANs

Before moving to our main algorithms, we briefly review GANs and WGANs. Consider the problem of learning a distribution $P_{r}$ on $\rm I\!R^{d}$, given data sampled from $P_{r}$. Suppose we have a DNN $G_{\theta}(\cdot)$ parameterized by $\theta$ that takes the random variable $z\in\rm I\!R^{m}$ as input and outputs a sample $G_{\theta}(z)\in\rm I\!R^{d}$. The random input $z$ is sampled from a prescribed distribution $Q_{z}$ (e.g., uniform, Gaussian), and we denote the distribution of $G_{\theta}(z)$ as $P_{g}$. We want to approximate $P_{r}$ with $P_{g}$.

GANs deal with this problem by defining a two-player zero-sum game between the generator $G_{\theta}(\cdot)$ and the discriminator $D_{\rho}(\cdot)$, which is another DNN parameterized by $\rho$. The discriminator takes a sample $x\in\rm I\!R^{d}$ as input and aims to tell if it is sampled from $P_{r}$ or $P_{g}$. Meanwhile, the generator tries to “deceive” the discriminator by mimicking the true distribution $P_{r}$. In vanilla GANs , the two-player zero-sum game is defined as follows:

Correspondingly, the loss functions for the generator and discriminator are

If the discriminator is optimal, $L_{g}$ measures the Jensen-Shannon (JS) divergence between $P_{r}$ and $P_{g}$ up to multiplication and addition by constants:

where $M=(P_{r}+P_{g})/2$, and $KL(\cdot||\cdot)$ is the Kullback-Leibler divergence . By training the generator and discriminator iteratively, ideally we can make $P_{g}$ approach $P_{r}$ in the sense of the JS divergence.

However, the JS divergence does not always provide a usable gradient for the generator, especially when the two distributions concentrate on low dimensional manifolds . As a consequence, training vanilla GANs is quite a delicate process and could be unstable . To fix this problem, Wasserstein GANs with clips on weights (weight-clipped WGANs) and WGANs with gradient penalty (WGAN-GP) were proposed. Similar to vanilla GANs, the two-play zero-sum game in WGANs is formulated as

The difference between weight-clipped WGANs and WGAN-GP is mainly on the technique of imposing the Lipschitz constraint for $D_{\rho}$: weight-clipped WGANs force $\rho$ to be bounded during the training, while WGAN-GP adds a gradient penalty to the loss function for the discriminatorIn and , the discriminators are names as “critics”, but in this paper we still use the name of discriminators for consistency with other versions of GANs, including vanilla GANs.. To be specific, in WGAN-GP the loss functions for the generator and discriminator are defined as

where $P_{\hat{x}}$ is the distribution generated by uniform sampling on straight lines between pairs of points sampled from $P_{r}$ and $P_{g}$, and $\lambda$ is the gradient penalty coefficient.

To demonstrate the effectiveness of WGAN-GP, we apply it on four toy problems of approximating distributions in $\rm I\!R^{10}$, as illustrated in Figure 2. We consider two types of $P_{r}$: one is a uniform distribution on a hypercube in $\rm I\!R^{10}$, while the other one is a uniform distribution concentrated on a curve embedded in $\rm I\!R^{10}$. For both cases of $P_{r}$, we test with an input of either 1-D or 10-D standard Gaussian noise, i.e., $Q_{z}=N(0,1)$ or $Q_{z}=N(0,I_{10})$. In all four cases, the generator converges and generates samples with distribution $P_{g}$ close to the real distribution $P_{r}$, even for the cases where the dimension of $Q_{z}$ mismatches the support of $P_{r}$. Due to its robustness, we use WGAN-GP as our default version of GANs in this paper.

Methodology

As a pedagogical problem, let us first consider the problem of modeling a stochastic process $f(x;\omega)$ on the domain $D\in\rm I\!R^{d}$ with GANs. We use the total $N$ snapshots of $f(x;\omega)$ sensor data as our training set:

be the generated “fake” snapshots, where $\{\bm{\xi}^{(j)}\}_{j}$ are instances of $\bm{\xi}$ with $j$ as the index.

The discriminator $D_{\rho}(\cdot)$ is implemented as another feed forward DNN parameterized by $\rho$. The GAN strategy shall be applied by feeding both the real and generated snapshots into the discriminator, and train the generator and discriminator iteratively with the Adam optimizer according to Equation (9). The detailed algorithm is presented in Algorithm 1.

2 Solving stochastic differential equations with PI-GANs

Consider Equation (1), as illustrated in Figure 3, solving SDEs with PI-GANs consists of the following three steps.

In the third step, we incorporate our training data to conduct adversarial training. The training data are collected as a group of snapshots in Equation (3). The corresponding generated “fake” snapshots are

where $\{\bm{\xi}^{(j)}\}_{j}$ are instances of $\bm{\xi}$ with $j$ as the index. We could then feed the generated snapshots and real snapshots into discriminator, and train the generators and discriminators iteratively. With well trained generators, we can then calculate all the statistics with sample paths created by the generators.

If the snapshots are collected in $M$ groups ($M>1$), we will also need to generate $M$ groups of “fake” snapshots:

where $t$s are the indices for the groups, $\bm{\xi}^{(t,j)}$ is an instance of $\bm{\xi}$ for each $(t,j)$, $\{x^{k,t}_{i}\}_{i=1}^{n_{k,t}}$ is the position setup of $n_{k,t}$ sensors for $k$ in group $t$ (similarly for other terms). We will also use multiple discriminators $\{D_{\rho_{t}}(\cdot)\}_{t}^{M}$, with each discriminator focusing on one group of snapshots, while the generators need to “deceive” all the discriminators simultaneously.

We give a formal and detailed description of our method in Algorithm 2. For the case where we only have one group of data, we set $M=1$. For simplicity, here we set the weight in the generator loss function $a_{t}=1$ for each $t$, which works well, but the method of setting $\{a_{t}\}_{t=1}^{M}$ requires further study.

Note that our method does not explicitly distinguish the three types of problems described Section 2. Solving forward problems, inverse problems or mixed problems actually uses the same framework.

Numerical Results

The following settings are commonly shared by all the test cases. We use tanh as the activation function instead of the commonly used ReLU activation function, because piece-wise linear functions are not suitable for solving SDEs, where we may need to take high order derivatives. All the feed forward DNNs for generators in the following numerical experiments have 4 hidden layers of width 128. The sizes of the input layer into the generators vary and will be specified case by case. Most of the discriminators also have hidden layers of the same size as the generators, except the ones in Section 5.1 that have 4 hidden layers of width 64, and the one for the additional group of snapshots on $u(x;\omega)$ in Section 5.3 that has 4 hidden layers of width 16. To initialize the DNNs, we use the uniform Xavier initializer for weights and zero initializer for biases. The distribution $Q_{z}$ of the noise input into the generators is an independent multi-variate standard Gaussian distribution. For the hyper-parameters in loss functions and optimizers, we use the default values of $\lambda=0.1,n_{d}=5,\alpha=0.0001,\beta_{1}=0,\beta_{2}=0.9$, as in the toy problem in . The sensors are placed equidistantly in the domain. Our algorithms are implemented with Tensorflow.

In this section, we test the problem of approximating stochastic processes. Consider the following Gaussian processes with zero mean and squared exponential kernel:

We consider three different choices of correlation length $l$: $l=1,0.5,0.2$, and two sensor numbers for $f(x;\omega)$: 11 and 6, and fix the number of snapshots to be 1000. For $l=0.2$ we also consider a supplementary case, where we have $1\times 10^{4}$ snapshots. The training sample paths and positions of sensors are illustrated in Figure 4. For each case, we run the code three times with different random seeds. The batch size in all the cases is 1000. The input layer of the generators has width of 5, i.e., the input noise is a four-dimensional random vector.

The results are illustrated in Figure 6, from which we conclude that:

When we fix the number of snapshots to be 1000, as the correlation length decreases, the gap between our generated processes and the reference solutions becomes wider. This is because smaller $l$ results in higher effective dimension, and more subtle local behavior of the stochastic processes; however, this gap could be narrowed if we increase the number of snapshots.

The approximations in the cases with 11 sensors are better than the approximations in the cases with 6 sensors if we have sufficient training data. This is reasonable since we need more sensors to describe the stochastic processes with small correlation length.

Despite the fact that the input noise into the generator is a four-dimensional vector, the spectra of the generated processes approximate the spectra of the target processes with much higher effective dimensionality. This makes sense since the low dimensional manifold could fold and twist itself to fill the high dimensional region, as discussed earlier.

1.2 WGAN-GP versus vanilla GANs

We compare the performance of WGAN-GP and the vanilla GANs in approximating stochastic processes for the following two cases:

A Gaussian process in Equation (17) with $l=0.2$.

A stochastic process with fixed 0 boundary condition. Specifically, we consider

In both cases, we use $1\times 10^{4}$ snapshots collected from 11 sensors as training data. In Figure 7 we show the means and standard deviations of generated processes trained by WGAN-GP and the vanilla GANs. We can see that both versions of GANs produce good approximations for case 1. However, the vanilla GANs fail in case 2, while WGAN-GP still generates a good approximation. This agrees with the theory in that vanilla GANs are not suitable for approximating distributions concentrated on low dimensional manifolds (the two fixed boundaries in this case), since the JS divergence cannot provide a usable gradient for the generators.

1.3 Overfitting issues

How about the overfitting of generators? In our problem, the overfitting of generators comes in two types:

: Overfitting in the random space: The distribution of generated snapshots is biased towards the empirical distribution of the training snapshots. In the worst case, the generated snapshots are concentrated on or near the support of training data.

: Overfitting in the physical space: The generated stochastic processes become worse after extensive training, and tend to match the real processes only at the sensor locations and display large variations where there is no sensor.

We first report that we did not detect type-2 overfitting in our experiments. Actually, this type of overfitting would be reflected on the mean and standard deviation of generated processes, which fit the reference values pretty well in our experiments. We attribute this to the property of GANs that the target for the generator is to approximate a distribution rather than a single point on the sensors. As a result, the generator does not need to overfit a specific value on the sensors in order to decrease the loss.

Finally, we point out that type-1 overfitting is less harmful than underfitting: in our problems, even in the worst case where the generated distribution is concentrated on or near the support of training data, we can still recover the sample paths whose snapshots are concentrated on or near the training snapshots, and based on these sample paths we can still obtain decent estimations.

2 Forward problem

We consider the stochastic elliptic equation (Equation (2)), and $k(x;\omega)$ and $f(x;\omega)$ as the following independent stochastic processes:

In this case, we need 13 dimensions to retain $99\%$ energy of $f(x;\omega)$. We put 13 $k$-sensors , 21 $f$-sensors and 2 $u$-sensors on the boundary of physical domain $\mathcal{D}$. The positions of the sensors and some sample paths of $k(x;\omega)$, $u(x;\omega)$ and $f(x;\omega)$ are illustrated in Figure 9.

To study the influence of input noise dimension, we fix the number of training snapshots to be 1000, and vary the input noise dimension to be 2, 4, 20 and 50. Subsequently, we fix the input noise dimension as 20 and vary the number of training snapshots to be 300, 1000 and 3000 to study the influence of the number of training snapshots. During the training process, we keep the batch size to be the total number of training snapshots. One notable condition here is the independence of $k(x;\omega)$ and $f(x;\omega)$. To reflect this, we shuffle the alignment of snapshots from $k(x;\omega)$ and $f(x;\omega)$ in each training step. For each case, we run the code three times with different random seeds. We stop the training after $1\times 10^{5}$ steps, and then select 33 generators and generate $1\times 10^{4}$ sample paths from each generator in the same way as in Section 5.1 to calculate the following statistics.

Our main quantity of interest in this problem is the mean and standard deviation of $u(x;\omega)$. In Figure 10 we show the relative error of our inferences, and compare it with the relative error of the stochastic collocation method and the Monte Carlo full trajectory sample paths of $u(x;\omega)$. We can see that:

Our errors are comparable with the errors calculated from 1000 full trajectory sample paths of $u(x;\omega)$, showing the effectiveness of our method considering that we only have 1000 snapshots on sparsely placed sensors for $k(x;\omega)$, $f(x;\omega)$ and boundary of $u(x;\omega)$.

The stochastic collocation method gives a better solution, but it requires a full knowledge of $k(x;\omega)$ and $f(x;\omega)$, including the covariance kernel function, which is far beyond our accessible data.

When we increase the dimension of input noise, we can see that the error of mean does not change too much, but the error of the standard deviation decreases. We attribute this to the fact that although a low dimensional manifolds could twist itself to fill in the high dimensional region, higher dimensional manifolds produce better approximations by filling in the high dimensional region more efficiently.

With fixed input noise dimension, the error decreases as the number of training snapshots increases.

We also illustrate the spectra of the generated processes for $f(x;\omega)$ in Figure 11 to verify that the generated processes captured the covariance structure of $f(x;\omega)$. We can see that the spectra of our generated processes fit the reference solution well. With fixed number of training snapshots, as we increase the input noise dimension, the gap between our generated processes and the reference solutions narrows. With fixed input noise dimension, the gap narrows as we increase the number of training snapshots. These observations on the spectra are similar with those about the error of $u$.

We further check the correlation between the generated processes. The correlation coefficient $C(f_{1}(x;\omega),f_{2}(x;\omega))$ between two stochastic processes $f_{1}(x;\omega)$ and $f_{2}(x;\omega)$ ($x\in D$) is defined as

where $x_{i}$ is the uniform grid points in $D$. Here, we set $n=201$ for $C(k,f)$, while $n=199$ for $C(k,u)$ to exclude the boundary points for $u$.

From Table 1 we can see that our generated processes for $k(x;\omega)$ and $f(x;\omega)$ are weakly correlated, while the generated processes for $k(x;\omega)$ and $u(x;\omega)$ are strongly correlated. All the correlation coefficients are close to the reference values.

2.2 Case 2: a (relatively) high dimensional problem

In this case, we solve the stochastic differential equation where the correlation length of $f(x;\omega)$ is relatively small. We consider $k(x;\omega)$ and $f(x;\omega)$ as the following independent stochastic processes:

Note that we need 30 dimensions to retain $99\%$ energy of $f(x;\omega)$. We put 13 $k$-sensors, 41 $f$-sensors, and 2 $u$-sensors on the boundary of our domain of interest $\mathcal{D}$. In this case, we use $1\times 10^{4}$ snapshots as training data while we set the batch size to be 1000 and fix the input noise dimension to be 20. Again, we shuffle the alignment of snapshots of $k(x;\omega)$ and $f(x;\omega)$ during the training to reflect their independence.

We run the code three times with different random seeds and stop the training after $2\times 10^{5}$ steps and then select 33 generators and generate $1\times 10^{4}$ sample paths from each generator in the same way as in Section 5.1 to calculate the error of $u(x;\omega)$ as well as the spectra of $f(x;\omega)$, as illustrated in Figure 12. We can see that the spectra fit the reference values well. Although the error of our inferred mean is larger than that calculated from $1\times 10^{4}$ sample paths of $u(x;\omega)$, the error of standard deviation is comparable.

3 Inverse and mixed problems

In this section, we show that our method can manage a wide range of problems, from forward problems to inverse problems, and mixed problems in between. In particular, we solve the three types of problems governed by Equation (2), where $k(x;\omega)$ and $f(x;\omega)$ are independent processes as follows:

We consider the following four cases of sensor placement:

: 1 $k$-sensor, 13 $u$-sensors (including 2 on the boundary), 13 $f$-sensors.

: 5 $k$-sensors, 9 $u$-sensors (including 2 on the boundary), 13 $f$-sensors.

: 9 $k$-sensors, 5 $u$-sensors (including 2 on the boundary), 13 $f$-sensors.

: 13 $k$-sensors, 2 $u$-sensors on the boundary of $u$, 13 $f$-sensors.

Note that case 1 is an inverse problem, case 4 is a forward problem, while case 2 and 3 represent mixed problems. For each case, we use 1000 snapshots for training. Note that in cases 1-3, we cannot shuffle the alignment of snapshots from $k(x;\omega)$ and $f(x;\omega)$ as in Section 5.2 since both $k(x;\omega)$ and $f(x;\omega)$ are correlated with $u(x;\omega)$. For consistency, in this section, we do not shuffle the alignment of snapshots for any of the four cases. We set the batch size to be 1000 and the input noise dimension to be 20. For each case, we run the code three times with different random seeds. The training stops after $2\times 10^{5}$ steps, and we select 33 generators and generate $1\times 10^{4}$ sample paths from each generator in the same way as in Section 5.1 to calculate the mean and standard deviation of $k(x;\omega)$ and $u(x;\omega)$.

In Figure 13 we compare the relative errors with reference solutions calculated from Monte Carlo sample paths as well as the method proposed in . We can see that our errors are in the same order of magnitude with errors from 1000 sample paths, showing the effectiveness of our method in solving all three types of problems. Also, for the case of 5 $k$-sensors and 9 $u$-sensors, our method achieves comparable accuracy with the method in .

4 Multiple groups of training data

Finally, we test our method for the case where we have multiple groups of snapshots as training data. In particular, we perform our test based on the case 4 in Section 5.3. Apart from the sensors in that case, we put one additional $u$-sensor at position $x=0$, and collect another group of 1000 snapshots at this additional sensor. Hence, we have two groups of training data:

: 13 $k$-sensors, 2 $u$-sensors on the boundary of $u$, 13 $f$-sensors. 1000 snapshots.

Again, we emphasize that the two groups of snapshots cannot be aligned. Note that a single snapshot in group 2 is almost useless since we cannot align it with another snapshot in group 1. However, the ensemble of the snapshots in group 2 actually can tell the distribution for $u(x;\omega)$ at $x=0$.

To utilize the data in two groups, we apply Algorithm 2 with the number of discriminators $M=2$. In this case, we set the training batch size as 1000 while the input noise dimension as 20. We run the code three times with different random seeds. The training stops after $2\times 10^{5}$ steps, and we select 33 generators and generate $1\times 10^{4}$ sample paths from each generator in the same way as in Section 5.1 to calculate the mean and standard deviation of $k(x;\omega)$ and $u(x;\omega)$. In Figure 14, we plot the errors of the inferred $k(x;\omega)$ and $u(x;\omega)$ in this case as well as errors from case 4 in Section 5.3. Compared with the results obtained by only using one group of data, in this case, the error of the standard deviation of $u(x)$ decreases significantly, showing the capability of our method in learning from the ensemble of snapshots in group 2. To the best of our knowledge, our method is so far the only one that can manage this case.

Summary and future work

We proposed physics-informed generative adversarial networks (PI-GANs) as a data-driven method for solving stochastic differential equations (SDEs) based on a limited number of scattered measurements. PI-GANs are composed of a discriminator, which is represented by a simple feed forward deep neural network (DNN) and of generators, which are a combination of feed forward DNNs and a neural network induced by the SDE. We assumed that partial data are available in terms of different realizations of the stochastic process obtained simultaneously at different locations in the domain. We trained the generators and discriminator iteratively with the loss functions employed in WGAN-GP so that the joint distribution of generated processes approximates the target stochastic processes. We also proposed a more general architecture with multiple discriminators to deal with cases, where data are collected in multiple groups, i.e., data collected independently from different sets of sensors.

We first tested WGAN-GP in approximating Gaussian processes for different correlation lengths. As shown in Figure 6, we obtained good approximation of the generated stochastic processes to the target ones even for a mismatch between the input noise dimensionality and the effective dimensionality of the target stochastic processes. The approximations were improved by increasing the number of sensors and snapshots. We also compared WGAN-GP and vanilla GANs, and concluded that vanilla GANs are not suitable for approximating stochastic processes with deterministic boundary condition, as shown in Figure 7. We further studied the overfitting issue by monitoring the negative discriminator loss (Figure 8(a)) and Wasserstein distance between empirical distributions (Figure 8(b)). We found that overfitting occurs also in the generator in addition to the discriminator as previously reported.

Subsequently, we considered the solution of elliptic SDEs requiring approximations of three stochastic processes, namely the solution $u(x;\omega)$, the forcing $f(x;\omega)$, and the diffusion coefficient $k(x;\omega)$. Without changing the framework, we were able to solve a wide range of problems, from forward to inverse problems, and in between, i.e., mixed problems where we have incomplete information for both the solution and the diffusion coefficient. As shown in Figure 10 and Figure 13, we obtained both the means and standard deviations of the stochastic solution and the diffusion coefficient in good agreement with benchmarks. In the case of the forward problem, we studied the influence of dimensionality of the input noise into the generators in Figure 10, and found that although small dimensionality can also work, high dimensionality for the input noise leads to better results. Moreover, we tested PI-GANs for a relatively high dimensional problem with $f(x;\omega)$ of stochastic dimension $30$. In Figure 12 we showed that the inferred mean and standard deviation of $u(x;\omega)$ as well as the spectra of $f(x;\omega)$ match the reference values well. Finally, we applied PI-GANs consisting of two discriminators to a stochastic problem with two groups of snapshots available for training, where the second group includes snapshots from only a single sensor for $u(x;\omega)$. We could see in Figure 14 that the error decreases compared with the error of the case where we only have the first group of data. This demonstrates the capability of PI-GANs to utilize information from multiple groups of data and learn from an ensemble of snapshots, even when a single snapshot is useless.

We also point out some limitations of the current version of PI-GANs. Since the computational cost of training GANs is much higher than training a single feed forward neural network, the PI-GANs method has higher computational cost than the physics-informed neural networks for deterministic PDEs and SDEs . Also, in our numerical experiments, overfitting was detected in both discriminators and generators when training data are limited. Although overfitting of generators is less harmful than underfitting, we wish to address this issue systematically in future work. Moreover, in the current work we only take into consideration the uncertainty described in SDEs, but not from measurements, nor do we account for the uncertainty of the approximability of GANs as was done in for DNNs using the dropout method. In future work, we wish to endow PI-GANs with uncertainty quantification coming from diverse sources, and hence provide estimates of the total uncertainty of the unknown distribution. Finally, we comment on the computational cost due to high dimensionality of stochastic problems. In experiments not reported here, we estimated the computational cost of solving the forward SDE for 60 dimensions, i.e., twice the dimensionality in the case reported in Section 5.2.2. Specifically, we decreased the correlation length for $f(x;\omega)$ by half and doubled the sensors for $f(x;\omega)$, while keeping all the other settings the same as in Section 5.2.2. We found that the computational cost was approximately doubled, obtaining $1-3\%$ of relative error in the mean and standard deviation of the stochastic solution. This is consistent with other lower dimensional problems we presented in this paper, which suggests that the overall computational cost increases with a low-polynomial growth, hence PI-GANs can, in principle, tackle very high dimensional stochastic problems. We will systematically investigate the scalability of PI-GANs for more complex stochastic PDEs in very high dimensions in future work. Moreover, we shall aim to optimize PI-GANs in terms of the architecture, depth and width of the generators and discriminators as well as the other hyperparameters.

Acknowledgement

We would like to acknowledge support from the Army Research Office (ARO) W911NF-18-1-0301, and from the Department of Energy (DOE) DE-SC0019434 and DE-SC0019453.

References

References