PDE-Net 2.0: Learning PDEs from Data with A Numeric-Symbolic Hybrid Deep Network

Earlier attempts on data-driven discovery of hidden physical laws include [4, 5]. Their main idea is to compare numerical differentiations of the experimental data with analytic derivatives of candidate functions, and apply the symbolic regression and the evolutionary algorithm to determining the nonlinear dynamic system. When the form of the nonlinear response function of a PDE is known, except for some scalar parameters, [6] presented a framework to learn these unknown parameters by introducing regularity between two consecutive time step using Gaussian process. Later in [7], a PDE constraint interpolation method was introduced to uncover the unknown parameters of the PDE model. An alternative approach is known as the sparse identification of nonlinear dynamics (SINDy) [8, 9, 10, 11, 12, 13]. The key idea of SINDy is to first construct a dictionary of simple functions and partial derivatives that are likely to appear in the equations. Then, it takes the advantage of sparsity promoting techniques (e.g. $\ell_{1}$ regularization) to select candidates that most accurately represent the data. In [14], the authors studied the problem of sea surface temperature prediction (SSTP). They assumed that the underlying physical model was an advection-diffusion equation. They designed a special neural network according to the general solution of the equation. Comparing with traditional numerical methods, their approach showed improvements in accuracy and computation efficiency.

Recent work greatly advanced the progress of PDE identification from observed data. However, SINDy requires to build a sufficiently large dictionary which may lead to high memory load and computation cost, especially when the number of model variables is large. Furthermore, the existing methods based on SINDy treat spatial and temporal information of the data separately and does not take full advantage of the temporal dependence of the PDE model. Although the framework presented by [6, 7] is able to learn hidden physical laws using less data than the approach based on SINDy, the explicit form of the PDEs is assumed to be known except for a few scalar learnable parameters. The approach of [14] is specifically designed for advection-diffusion equations, and cannot be readily extended to other types of equations. Therefore, extracting governing equations from data in a less restrictive setting remains a great challenge.

The main objective of this paper is to design a transparent deep neural network to uncover hidden PDE models from observed complex dynamic data with minor prior knowledge on the mechanisms of the dynamics, and to perform accurate predictions at the same time. The reason we emphasize on both model recovery and prediction is because: 1) the ability to conduct accurate long-term prediction is an important indicator of accuracy of the learned PDE model (the more accurate is the prediction, the more confident we have on the underlying recovered PDE model); 2) the trained neural network can be readily used in applications and does not need to be re-trained when initial conditions are altered. Our inspiration comes from the latest development of deep learning techniques in computer vision. An interesting fact is that some popular networks in computer vision, such as ResNet[15, 16], have close relationship with ODEs/PDEs and can be naturally merged with traditional computational mathematics in various tasks [17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. However, existing deep networks designed in deep learning mostly emphasis on expressive power and prediction accuracy. These networks are not transparent enough to be able to reveal the underlying PDE models, although they may perfectly fit the observed data and perform accurate predictions. Therefore, we need to carefully design the network by combining knowledge from deep learning and numerical PDEs.

The proposed deep neural network is an upgraded version of our original PDE-Net [1]. The main difference is the use of a symbolic network to approximate the nonlinear response function, which significantly relaxes the requirement on the prior knowledge on the PDEs to be recovered. During training, we no longer need to assume the general type of the PDE (e.g. convection, diffusion, etc.) is known. Furthermore, due to the lack of prior knowledge on the general type of the unknown PDE models, more carefully designed constraints on the convolution filters as well as the parameters of the symbolic network are introduced. We refer to this upgraded network as PDE-Net 2.0.

Assume that the PDE to be recovered takes the following generic form

PDE-Net 2.0 is designed as a feed-forward network by discretizing the above PDE using forward Euler in time and finite difference in space. The forward Euler approximation of temporal derivative makes PDE-Net 2.0 ResNet-like [15, 18, 21], and the finite difference is realized by convolutions with trainable kernels (or filters). The nonlinear response function $F$ is approximated by a symbolic neural network, which shall be referred to as $SymNet$. All the parameters of the $SymNet$ and the convolution kernels are jointly learned from data. To grant full transparency to the PDE-Net 2.0, proper constraints are enforced on the $SymNet$ and the filters. Full details on the architecture and constraints will be presented in Section 2.

Data-driven discovery of hidden physical laws and model reduction have a lot in common. Both of them concern on representing observed data using relatively simple models. The main difference is that, model reduction emphasis more on numerical precision rather than acquiring the analytic form of the model.

It is common practice in model reduction to use a function approximator to express the unknown terms in the reduced models, such as approximating subgrid stress for large-eddy simulation[27, 28, 29] or approximating interatomic forces for coarse-grained molecular dynamic systems[30, 31]. Our work may serve as an alternative approach to model reduction and help with analyzing the reduced models.

The particular novelties of our approach are that we impose appropriate constraints on the learnable filters and use a properly designed symbolic neural network to approximate the response function $F$. Using learnable filters makes the PDE-Net 2.0 more flexible, and enables more powerful approximation of unknown dynamics and longer time prediction (see numerical experiments in Section 3 and Section 4). Furthermore, the constraints on the learnable filters and the use of a deep symbolic neural network enable us to uncover the analytic form of $F$ with minor prior knowledge on the dynamic, which is the main advantage of PDE-Net 2.0 over the original PDE-Net. In addition, the composite representation by the symbolic network is more efficient and flexible than SINDy. Therefore, the proposed PDE-Net 2.0 is distinct from the existing learning based methods to discover PDEs from data.

Inspired by the dynamic system perspective of deep neural networks [17, 18, 19, 20, 21, 22], we consider forward Euler as the temporal discretization of the evolution PDE (1), and unroll the discrete dynamics to a feed-forward network. One may consider more sophisticated temporal discretization which naturally leads to different network architectures [21]. For simplicity, we focus on forward Euler in this paper.

In the original PDE-Net [1], the learnable filters are properly constrained so that we can easily identify their correspondence to differential operators. PDE-Net 2.0 adopts the same constrains as the original version of the PDE-Net. For completeness, we shall review the related notions and concepts, and provide more details.

A profound relationship between convolutions and differentiations was presented by [32, 33], where the authors discussed the connection between the order of sum rules of filters and the orders of differential operators. Note that the definition of convolution we use follows the convention of deep learning, which is defined as

This is essentially correlation instead of convolution in the mathematics convention. Note that if $f$ is of finite size, we use periodic boundary condition.

In practical implementation, the filters are normally finite and can be understood as matrices. For an $N\times N$ filter $q$ ($N$ being an odd number), assuming the indices of $q$ start from $-\frac{N-1}{2}$, (3) can be written in the following simpler form

The following proposition from [33] links the orders of sum rules with orders of differential operators.

According to Proposition 2.1, an $\alpha$th order differential operator can be approximated by the convolution of a filter with $\alpha$ order of sum rules. Furthermore, according to (5), one can obtain a high order approximation of a given differential operator if the corresponding filter has an order of total sum rules with $K>|\alpha|+k,k\geqslant 1$. For example, consider filter

It has a sum rules of order $(1,0)$, and a total sum rules of order $3\backslash\{2\}$. Thus, up to a constant and a proper scaling, $q$ corresponds to a discretization of $\frac{\partial}{\partial x}$ with second order accuracy.

For an $N\times N$ filter $q$, define the moment matrix of $q$ as

where

We shall call the $(i,j)$-element of $M(q)$ the $(i,j)$-moment of $q$ for simplicity. For any smooth function $f:\mathbb{R}^{2}\to\mathbb{R}$, we apply convolution on the sampled version of $f$ with respect to the filter $q$. By Taylor’s expansion, one can easily obtain the following formula

From (8) we can see that filter $q$ can be designed to approximate any differential operator with prescribed order of accuracy by imposing constraints on $M(q)$.

For example, if we want to approximate $\frac{\partial u}{\partial x}$ (up to a constant) by convolution $q\circledast u$ where $q$ is a $3\times 3$ filter, we can consider the following constrains on $M(q)$:

Here, $\star$ means no constraint on the corresponding entry. The constraints described by the moment matrix on the left of (9) guarantee the approximation accuracy is at least first order, and the one on the right guarantees an approximation of at least second order. In particular, when all entries of $M(q)$ are constrained, e.g.

the corresponding filter can be uniquely determined. In the PDE-Net 2.0, all filters are learned subjected to partial constraints on their associated moment matrices, with at least second order accuracy.

The symbolic neural network $SymNet_{m}^{k}$ of the PDE-Net 2.0 is introduced to approximate the multivariate nonlinear response function $F$ of (1). Neural networks have recently been proven effective in approximating multivariate functions in various scenarios [36, 37, 38, 39, 40]. For the problem we have in hand, we not only require the network to have good expressive power, but also good transparency so that the analytic form of $F$ can be readily inferred after training. Our design of $SymNet_{m}^{k}$ is motivated by EQL/EQL^÷ proposed by [41, 42].

Figure 3 shows the symbolic neural network with two hidden layers, i.e. $SymNet_{m}^{2}$, where $f$ is a dyadic operation unit, e.g. multiplication or division. In this paper, we focus on multiplication, i.e. we take $f(a,b)=a\times b$. Different from EQL/EQL^÷, each hidden layer of the $SymNet_{m}^{k}$ directly takes the outputs from the preceding layer as inputs, rather than a linear combination of them. Furthermore, it adds one additional variable (i.e. $f(\cdot,\cdot)$) at each hidden layer. To better understand $SymNet_{m}^{k}$, we present an example in Algorithm 1 showing how $SymNet_{6}^{2}$ is constructed. In particular, when $b^{i}=0,i=1,2,3$,

and $W^{3}=(0,0,0,0,0,0,-1,-1)$, then $SymNet_{6}^{2}(u,u_{x},u_{y},v,v_{x},v_{y})=-uu_{x}-vu_{y}$ which is the right-hand-side of $\frac{\partial u}{\partial t}$ of the Burgers’ equation without viscosity.

The $SymNet_{m}^{k}$ can represent all polynomials of variables $(x_{1},x_{2},\ldots,x_{m})$ with the total number of multiplications not exceeding $k$. If needed, one can add more operations to the $SymNet_{m}^{k}$ to increase the capacity of the network.

For example, $\sum_{i=1}^{m}x_{i}+\sum_{i=1}^{k}x_{i}x_{i+1}$ and $\prod_{i=1}^{k+1}x_{i}$ with $k<m$ are all members of $\mathcal{P}^{k}[x_{1},\cdots,x_{m}]$. The elements in $\mathcal{P}^{k}[x_{1},\cdots,x_{m}]$ are of simple forms when $k$ is relatively small. The following proposition shows that $SymNet_{m}^{k}$ can represent all polynomials of variables $(x_{1},x_{2},\ldots,x_{m})$ with the total number of multiplications not exceeding $k$. Note that the actual capacity of $SymNet_{m}^{k}$ is larger than $\mathcal{P}^{k}[x_{1},\cdots,x_{m}]$, i.e. $\mathcal{P}^{k}[x_{1},\cdots,x_{m}]$ is a subset of the set of functions that $SymNet_{m}^{k}$ can represent.

Proof: We prove this proposition by induction. When $k=1$, the conclusion obviously holds. Suppose the conclusion holds for $k$. For any polynomial $P\in\mathcal{P}^{k+1}[x_{1},\cdots,x_{m}]$, we only need to consider the cases when $P$ has a total number of multiplications greater than 1.

We take any monomial of $P$ that has degree greater than 1, which we suppose take the form $x_{1}x_{2}\cdot A$ where $A$ is a monomial of variable $(x_{1},\ldots,x_{m})$. Then, $P$ can be written as $P=x_{1}x_{2}A+Q$. Define new variable $x_{m+1}=x_{1}x_{2}$. Then, we have

By the induction hypothesis, there exists a set of parameters such that $P=SymNet_{m+1}^{k}(x_{1},\cdots,x_{m+1})$.

We take the linear transform between the input layer and the first hidden layer of $SymNet_{m}^{k+1}$ as

Then, the output of the first hidden layer is $x_{1},x_{2},\cdots,x_{m},x_{1}x_{2}$. If we use it as the input of $SymNet_{m+1}^{k}$, we have

which concludes the proof. $\square$

We use the following example to show the advantage of $SymNet$ over SINDy. Consider two variables $u,v$ and all of their derivatives of order $\leq 2$:

Suppose the polynomial to be approximated is $P=-uu_{x}-vu_{y}+u_{xx}$. For $k=l=3$, the size of the dictionary of SINDy is $\binom{15}{3}=455$ and the computation of linear combination of the elements requires $909$ flops. The memory load of $SymNet_{12}^{3}$, however, is 15 and an evaluation of the network requires 180 flops. Therefore, $SymNet$ can significantly reduce memory load and computation cost when input data is large. Note that when $k$ is large and $l$ small, $SymNet_{m}^{k}$ is worse than SINDy. However, for system identification problems, we normally wish to obtain a compact representation (i.e. smaller $k$). Thus, $SymNet$ takes full advantage of this prior knowledge and can significantly save on memory and computation cost which is crucial in the training of the PDE-Net 2.0.

We adopt the following loss function to train the proposed PDE-Net 2.0:

where the hyper-parameters $\lambda_{1}$ and $\lambda_{2}$ are chosen as $\lambda_{1}=0.001$ and $\lambda_{2}=0.005$. Now, we present details on each of the term of the loss function and introduce pseudo-upwind as an additional constraint on PDE-Net 2.0.

In the PDE-Net 2.0, parameters can be divided into three groups: 1) moment matrices to generate convolution kernels; 2) the parameters of $SymNet$; and 3) hyper-parameters, such as the number of filters, the size of filters, the number of $\delta t$-Blocks and number of hidden layers of $SymNet$, regularization weights $\lambda_{1},\lambda_{2},\lambda_{3}$, etc. The parameters of the $SymNet$ are shared across the computation domain $\Omega$, and are initialized by random sampling from a Gaussian distribution. For the filters, we initialize $D_{01},D_{10}$ with second order pseudo-upwind scheme and central difference for all other filters. For example, if the size of the filters were set to be $5\times 5$, then the initial values of the convolution kernels $q_{01},q_{02}\in\mathbb{R}^{5\times 5}$ are

We use layer-wise training to train the PDE-Net 2.0. We start with training the PDE-Net 2.0 on the first $\delta$t-block with a batch of data, and then use the results of the first $\delta$t-block as the initialization and restart training on the first two $\delta$t-blocks with another batch. Repeat this procedure until we complete all $n$ $\delta$t-blocks. Note that all the parameters in each of the $\delta$t-block are shared across layers. In addition, we add a warm-up step before the training of the first $\delta$t-block by fixing filters and setting regularization term to be 0 (i.e. $\lambda_{1}=\lambda_{2}=0$). The warm-up step is to obtain a good initial guess of the parameters of $SymNet$.

To demonstrate the necessity of having learnable filters, we will compare the PDE-Net 2.0 containing learnable filters with the PDE-Net 2.0 having fixed filters. To differentiate the two cases, we shall call the PDE-Net 2.0 with fixed filters the “Frozen-PDE-Net 2.0”. Note that for Frozen-PDE-Net 2.0, the filters are fixed to be the initial values we choose to train the regular PDE-Net 2.0. This is a natural choice since when we know apriori that the PDE is Burgers’ equation, it would be a stable finite difference scheme. However, intuitively speaking, freezing any finite difference approximations of the differential operators during training of PDE-Net 2.0 is not ideal, because you cannot possibly know which numerical scheme to use without knowing the form of the PDE. Therefore, for inverse problem, it is better to learn both the PDE model and the discretization of the PDE model simultaneously. This assertion is supported by our emperical comparisons between frozen and regular PDE-Net 2.0 in Table 1 and 3.

In this section we consider a 2-dimensional Burger’s equation with periodic boundary condition on $\Omega=[0,2\pi]\times[0,2\pi]$,

with $(t,x,y)\in[0,4]\times\Omega,$ where $\nu=0.05$

The training data is generated by a finite difference scheme on a $128\times 128$ mesh and then restricted to a $32\times 32$ mesh. The temporal discretization is 2nd order Runge-Kutta with time step $\delta t=\frac{1}{1600}$, the spatial discretization uses a 2nd order upwind scheme for $\nabla$ and the central difference scheme for $\Delta$. The initial value $u_{0}(x,y),v_{0}(x,y)$ takes the following form,

where $w_{0}(x,y)=\sum_{|k|,|l|\leq 4}\lambda_{k,l}\cos(kx+ly)+\gamma_{k,l}\sin(kx+ly)$, $\lambda_{k,l},\gamma_{k,l}\sim\mathcal{N}(0,1),c\sim\mathcal{U}(-2,2)$. Here, $\mathcal{N}(0,1),\mathcal{U}(-2,2)$ represents the standard normal distribution and uniform distribution on $[-2,2]$ respectively. We also add noise to the generated data:

where $M=\max_{x,y,t}\{U(t,x,y)\}$, $W\sim\mathcal{N}(0,1)$.

Suppose we know a priori that the order of the underlying PDE is no more than 2, we can use two $SymNet_{12}^{5}$ to approximate the right-hand-side nonlinear response function of (10) component-wise. Let $U=(u,v)^{\top}$. We denote the two $SymNet_{12}^{5}$ as $Net_{u}$ and $Net_{v}$ respectively. Then, each $\delta t$-block of the PDE-Net 2.0 can be written as

where $\{D_{ij}:0\leq i+j\leq 2\}$ are convolution operators.

During training and testing, the data is generated on-the-fly. The size of the filters that will be used is $5\times 5$. The total number of parameters in $Net_{u}$ and $Net_{v}$ is 336, and the number of trainable parameters in moment matrices is 105 for $5\times 5$ filters ($6\times 5\times 5-45$ (constraint on moment)). During training, we use BFGS, instead of SGD, to optimize the parameters. We use 28 data samples per batch to train each $\delta t$-block and we only construct the PDE-Net 2.0 up to 9 layers, which requires totally 420 data samples during the whole training procedure.

We first demonstrate the ability of the trained PDE-Net 2.0 to recover the analytic form of the unknown PDE model. We use the symbolic math tool in python to obtain the analytic form of $SymNet$. Results are summarized in Table 1. As one can see from Table 1 that we can recover the terms of the Burgers’ equation with good accuracy, and using learnable filters helps with the identification of the PDE model. Furthermore, the terms that are not included in the Burgers’ equation all have relatively small weights in the $SymNet$ (see Figure 7).

We also demonstrate the ability of the trained PDE-Net 2.0 in prediction, i.e. the ability to generalize. After the PDE-Net 2.0 with $n$ $\delta t$-blocks ($1\leq n\leq 9$) is trained, we randomly generate 1000 initial guesses based on (11) and (12), feed them to the PDE-Net 2.0, and measure the relative error between the predicted dynamics (i.e. the output of the PDE-Net 2.0) and the actual dynamics (obtained by solving (10) using high precision numerical scheme). The relative error between the true data $U$ and the predicted data $\tilde{U}$ is defined as $\epsilon=\frac{\|\tilde{U}-U\|_{2}^{2}}{\|U-\bar{U}\|_{2}^{2}},$ where $\bar{U}$ is the spatial average of $U$. The error plots are shown in Figure 4. Some of the images of the predicted dynamics are presented in Figure 5 and the errors maps are presented in Figure 6. As we can see that, even trained with noisy data, the PDE-Net 2.0 is able to perform long-term prediction. Having multiple $\delta t$-blocks indeed improves predict accuracy. Furthermore, PDE-Net 2.0 performs significantly better than Frozen-PDE-Net 2.0. This suggests that when we do not know the PDE, we cannot possibly know how to properly discretize it. Therefore, for inverse problems, it is better to learn both the PDE model and its discretization simultaneously.

This subsection demonstrates the importance of enforcing sparsity on the $SymNet$ and using pseudo-upwind. As we can see from Figure 7 that having sparsity constraints on the $SymNet$ helps with suppressing the weights on the terms that do not exist in the Burgers’ equation. Furthermore, Figure 8 and Figure 9 show that having sparsity constraint on the $SymNet$ or using pseudo-upwind can significantly reduce prediction errors.