Title: PDE-NET: LEARNING PDES FROM DATA

Abstract: Partial differential equations (PDEs) play a prominent role in many disciplines such as applied mathematics, physics, chemistry, material science, computer science, etc. PDEs are commonly derived based on physical laws or empirical observations. However, the governing equations for many complex systems in modern applications are still not fully known. With the rapid development of sensors, computational power, and data storage in the past decade, huge quantities of data can be easily collected and efficiently stored. Such vast quantity of data offers new opportunities for data-driven discovery of hidden physical laws. Inspired by the latest development of neural network designs in deep learning, we propose a new feed-forward deep network, called PDE-Net, to fulfill two objectives at the same time: to accurately predict dynamics of complex systems and to uncover the underlying hidden PDE models. The basic idea of the proposed PDE-Net is to learn differential operators by learning convolution kernels (filters), and apply neural networks or other machine learning methods to approximate the unknown nonlinear responses. Comparing with existing approaches, which either assume the form of the nonlinear response is known or fix certain finite difference approximations of differential operators, our approach has the most flexibility by learning both differential operators and the nonlinear responses. A special feature of the proposed PDE-Net is that all filters are properly constrained, which enables us to easily identify the governing PDE models while still maintaining the expressive and predictive power of the network. These constrains are carefully designed by fully exploiting the relation between the orders of differential operators and the orders of sum rules of filters (an important concept originated from wavelet theory). We also discuss relations of the PDE-Net with some existing networks in computer vision such as Network-In-Network (NIN) and Residual Neural Network (ResNet). Numerical experiments show that the PDE-Net has the potential to uncover the hidden PDE of the observed dynamics, and predict the dynamical behavior for a relatively long time, even in a noisy environment.

1. INTRODUCTION:

Differential equations, especially partial differential equations(PDEs), play a prominent role in many disciplines to describe the governing physical laws underlying a given system of interest. Traditionally, PDEs are derived based on simple physical principles such as conservation laws, minimum energy principles, or based on empirical observations. Important examples include the Navier-Stokes equations in fluid dynamics, the Maxwell's equations for electromagnetic propagation, and the Schrödinger's equations in quantum mechanics. However, many complex systems in modern applications (such as many problems in climate science, neuroscience, finance, etc.) still have eluded mechanisms, and the governing equations of these systems are only partially known. With the rapid development of sensors, computational power, and data storage in the last decade, huge quantities of data can be easily collected and efficiently stored . Such vast quantity of data offers new opportunities for data-driven discovery of potentially new physical laws. Then, one may ask the following interesting and intriguing question: can we learn a PDE model (if there exists one) from a given data set and perform accurate and efficient predictions using the learned model? One of earlier attempts on data-driven discovery of hidden physical laws is by Bongard & Lipson (2007) and Schmidt & Lipson (2009) . 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 dynamical system. Recently, Brunton et al. (2016) , Schaeffer (2017) , Rudy et al. (2017) and Wu & Zhang (2017) propose an alternative approach using sparse regression. They construct a dictionary of simple functions and partial derivatives that are likely to appear in the unknown governing equations. Then, they take advantage of sparsity promoting techniques to select candidates that most accurately represent the data. When the form of the nonlinear response of a PDE is known, except for some scalar parameters, Raissi & Karniadakis (2017) presented a framework to learn these unknown parameters by introducing regularity between two consecutive time step using Gaussian process. More recently, Raissi et al. (2017) introduced a new class of universal function approximators called the physics informed neural networks which is capable of discovering nonlinear PDEs parameterized by scalars.

These recent work greatly advanced the progress of the problem. However, symbolic regression is expensive and does not scale very well to large systems. The sparse regression method requires to fix certain numerical approximations of the spatial differentiations in the dictionary beforehand, which limits the expressive and predictive power of the dictionary. Although the framework presented by Raissi & Karniadakis (2017) ; Raissi et al. (2017) is able to learn hidden physical laws using less data than the approach of sparse regression, the explicit form of the PDEs are assumed to be known except for a few scalar learnable parameters. Therefore, extracting governing equations from data in a less restrictive setting remains a great challenge.

The main objective of this paper is to accurately predict the dynamics of complex systems and to uncover the underlying hidden PDE models (should they exist) at the same time, with minimal prior knowledge on the systems. 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 (He et al., 2016a; b) , have close relationship with PDEs (Chen et al., 2015; E, 2017; Haber & Ruthotto, 2017; Sonoda & Murata, 2017; Lu et al., 2017) . Furthermore, the deeper is the network, the more expressive power the network possesses, which may enable us to learn more complex dynamics arose from fields other than computer vision. 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 applied mathematics so that we can learn the governing PDEs of the dynamics and make accurate predictions at the same time. Note that our work is closely related to Chen et al. (2015) where the authors designed their network based on discretization of quasilinear parabolic equations. However, it is not clear if the dynamics of image denoising has to be governed by PDEs, nor did the authors attempt to recover the PDE (should there exists one).

In this paper, we design a deep feed-forward network, named PDE-Net, based on the following generic nonlinear evolution PDE

u t = F (x, u, ∇u, ∇ 2 u, . . .), x ∈ Ω ⊂ R 2 , t ∈ [0, T ].

The objective of the PDE-Net is to learn the form of the nonlinear response F and to perform accurate predictions. Unlike the existing work, the proposed network only requires minor knowledge on the form of the nonlinear response function F , and requires no knowledge on the involved differential operators (except for their maximum possible order) and their associated discrete approximations. The nonlinear response function F can be learned using neural networks or other machine learning methods, while discrete approximations of the differential operators are learned using convolution kernels (i.e. filters) jointly with the learning of the response function F . If we have a prior knowledge on the form of the response function F , we can easily adjust the network architecture by taking advantage of the additional information. This may simplify the training and improve the results. We will also discuss relations of the PDE-Net to some existing networks in computer vision such as Network-In-Network (NIN) and ResNet. Details are given in Section 2.

In Section 3 and Section 4, we conduct numerical experiments on a linear PDE (convection-diffusion equation) and a nonlinear PDE (convection-diffusion equation with a nonlinear source). We generate data set for each PDE using high precision numerical methods and add Gaussian noise to mimic real situations. Our numerical results show that the PDE-Net can uncover the hidden equations of the observed dynamics, and can predict the dynamical behavior for a relatively long time, even in a noisy environment.

A particular novelty of our approach is that we impose appropriate constraints on the learnable filters in order to easily identify the governing PDE models while still maintaining the expressive and pre-dictive power of the network. This makes our approach different from existing deep convolutional networks which mostly emphasis on the prediction accuracy of the networks, as well as all the existing approaches of learning PDEs from data which assume either the form of the response function is known or have fixed approximations of the differential operators. In other words, our proposed approach not only has vast flexibility in fitting observed dynamics and is able to accurately predict its future behavior, but is also able to reveal the hidden equations driving the observed dynamics. The constraints on the filters are motivated by the earlier work of Cai et al. (2012) ; Dong et al. (2017) where general relations between wavelet frame transforms and differential operators were established. In particular, it was observed in Dong et al. (2017) that we can relate filters and finite difference approximation of differential operators by examining the orders of sum rules of the filters (an important concept in wavelet theory and closely related to vanishing moments of wavelet functions). These constraints on the filters may also be useful in network designs for machine learning tasks in computer vision.

2. PDE-NET: A FLEXIBLE DEEP ARCHTECTURE TO LEARN PDES FROM DATA:

Given a series of measurements of some physical quantities {u(t, •) :

t = t 0 , t 1 , • • • } on the spatial domain Ω ⊂ R 2 , with u(t, •) : Ω → R,

we want to discover the governing PDEs of the data. We assume that the observed data are associated with a PDE that takes the following general form:

u t (t, x, y) = F (x, y, u, u x , u y , u xx , u xy , u yy , . . .), (x, y) ∈ Ω ⊂ R 2 , t ∈ [0, T ].

(1) Our objective is to design a feed-forward network, named the PDE-Net, that approximates the PDE (1) in the way that: 1) we can predict the dynamical behavior of the equation for as long time as possible; 2) we are able to reveal the form of the response function F and the differential operators involved. There are two main components of the PDE-Net that are combined together in the same network: one is automatic determination on the differential operators involved in the PDE and their discrete approximations; the other is to approximate the nonlinear response function F . In this section, we start with discussions on the relation between convolutions and differentiations in discrete setting.

2.1. CONVOLUTIONS AND DIFFERENTIATIONS:

A comprehensive analysis on the relations between convolutions and differentiations within variational and PDE framework were laid out by Cai et al. (2012) and Dong et al. (2017) , where the authors established general connections between PDE based approach and wavelet frame based approach for image restoration problems. We demonstrate one of the key observations of their work using a simple example. Consider the 2-dimensional Haar wavelet frame filter bank contains one low-pass filter h 00 and three high pass filters h 10 , h 01 and h 11 :

h 00 = 1 4 1 1 1 1 , h 10 = 1 4 1 -1 1 -1 , h 01 = 1 4 1 1 -1 -1 , h 11 = 1 4 1 -1 -1 1 .

The associated Haar wavelet frame transform on an image u is defined by

W u = {h ij [-•] u : 0 ≤ i, j ≤ 1}, where

is the circular convolution. It is easy to verify using Taylor's expansion that the high frequency coefficients of the Haar wavelet frame transform on u are discrete approximations of differential operators:

h 10 [-•] u ≈ 1 2 δ x u x , h 01 [-•] u ≈ 1 2 δ y u y , h 11 [-•] u ≈ 1 4 δ x δ y u xy .

Here, δ x and δ y represent the horizontal and vertical spatial grid size respectively. For simplicity of notation, we use regular character to denote both discrete and continuum functions, since there should be no confusion within the context.

A profound relationship between convolutions and differentiations was presented in Dong et al. (2017) , where the authors discussed the connection between the order of sum rules of filters and the orders of differential operators. Note that the order of sum rules is closely related to the order of vanishing moments in wavelet theory (Daubechies, 1992; Mallat, 1999) . We first recall the definition of the order of sum rules.

Definition 2.1 (Order of Sum Rules). For a filter q, we say q to have sum rules of order α = (α 1 , α 2 ), where α ∈ Z 2 + , provided that

k∈Z 2 k β q[k] = 0 (2)

for all β ∈ Z 2 + with |β| < |α| and for all β ∈ Z 2 + with |β| = |α| but β = α. If (2) holds for all β ∈ Z 2 + with |β| < K except for β = β 0 with certain β 0 ∈ Z 2 + and |β 0 | = J < K, then we say q to have total sum rules of order K\{J + 1}.

The following proposition from Dong et al. (2017) links the orders of sum rules with orders of differential operator. Propositin 2.1. Let q be a filter with sum rules of order α ∈ Z 2 + . Then for a smooth function

EQUATION

where C α is the constant defined by

C α = 1 α! k∈Z 2 k α q[k].

If, in addition, q has total sum rules of order K\{|α| + 1} for some K > |α|, then

EQUATION

According to Proposition 2.1, an αth order differential operator can be approximated by the convolution of a filter with α order of sum rules. Furthermore, according to (4), 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 > |α| + k, k 1. For example, the filter h 10 in the Haar wavelet frame filter bank has a sum rules of order (1, 0), and a total sum rules of order 2\{2}. Thus, up to a constant and a proper scaling, h 10 corresponds to a discretization of ∂ ∂x with first order approximation. The filer h 11 has a sum rules of order (1, 1), and a total sum rules of order 3\{3}. Thus, up to a constant and a proper scaling, h 11 corresponds to a discretization of ∂ 2 ∂x∂y with first order approximation. Finally, consider filter

q = 1 0 -1 2 0 -2 1 0 -1 .

It has a sum rules of order (1, 0), and a total sum rules of order 3\{2}. Thus, up to a constant and a proper scaling, q corresponds to a discretization of ∂ ∂x with second order approximation. Now, we introduce the concept of moment matrix for a given filter that will be used to constrain filters in the PDE-Net. For an N × N filter q, define the moment matrix of q as

EQUATION

for i, j = 1, 2, . . . , N . We shall call the (i, j)-element of M (q) the (i -1, j -1)-moment of q for simplicity. Combining (5) and Proposition 2.1, one can easily see that filter q can be designed to approximate any differential operator at any given approximation order by imposing constraints on M (q). For example, if we want to approximate ∂u ∂x (up to a constant) by convolution q u where q is a 3 × 3 filter, we can consider the following constrains on M (q): 0 0 1 or

EQUATION

)

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

M (q) = 0 0 0 1 0 0 0 0 0 , the corresponding filter can be uniquely determined, in which case we call it a "frozen" filter. In the PDE-Net which shall be introduced in the next subsection, all filters are learned subjected to partial constraints on their associated moment matrices.

It is worth noticing that the approximation property of a filter is limited by its size. Generally speaking, large filters can approximate higher order differential operators or lower order differential operators with higher approximation orders. Taking 1-dimensional case as an example, 3-element filters cannot approximate the fifth order differential operator, whereas 7-element filters can. In other words, the larger are the filters, the stronger is the representation capability of filters. However, larger filters lead to more memory overhead and higher computation cost. It is a wisdom to balance the trade-off in practice.

2.2. ARCHITECTURE OF PDE-NET:

Given the evolution PDE (1), we consider forward Euler as the temporal discretization. One may consider more sophisticated temporal discretization which leads to different network architectures.

For simplicity, we focus on forward Euler in this paper.

ONE δt-BLOCK::

Let ũ(t i+1 , •) be the predicted value of u at time t i+1 based on the value of u at t i . Then, we have

EQUATION

Here, the operators D 0 and D ij are convolution operators with the underlying filters denoted by q 0 and q ij , i.e. D 0 u = q 0 u and D ij u = q ij u. The operators D 10 , D 01 , D 11 , etc. approximate differential operators, i.e. D ij u ≈ ∂ i+j u ∂ i x∂ j y . The operators D 0 and D 00 are average operators. The purpose of introducing these average operators in stead of using the identity is to improve stability of the network and enables it to capture more complex dynamics. Other than the assumption that the observed dynamics is governed by a PDE of the form (1), we assume that the highest order of the PDE is less than some positive integer. Then, the task of approximating F is equivalent to a multivariate regression problem, which can be approximated by a point-wise neural network (with shared weights across the computation domain Ω) or other classical machine learning methods. Combining the approximation of differential operators and the nonlinear function F , we achieve an approximation framework of (7) which will be referred to as a δt-block (see Figure 1 ). Note that if we have a prior knowledge on the form of the response function F , we can easily adjust the network architecture by taking advantage of the additional information. This may simplify the training and improve the results.

PDE-NET (MULTIPLE δt-BLOCKS)::

One δt-block only guarantees the accuracy of one-step dynamics, which does not take error accumulation into consideration. This may cause severe instability in prediction. To improve the stability of the network and enable long-term prediction, we stack multiple δt-blocks into a deep network, and call this network the PDE-Net (see Figure 2 ). The importance of stacking multiple δt-blocks will be demonstrated in Section 3.

The PDE-Net can be easily described as: (1) stacking one δt-block multiple times; (2) sharing parameters in all δt-blocks. Given an input data u(t i , •), training a PDE-Net with n δt-blocks needs to minimize the accumulated error ||u(t i+n , •) -ũ(t i+n , •)|| 2 2 , where ũ(t i+n , •) is the output from the PDE-Net (i.e. n δt-blocks) with input u(t i , •). Thus, the PDE-Net with bigger n owns a longer time stability. Note that sharing parameters is a common practice in deep learning, which decreases the number of parameters and leads to significant memory reduction (Goodfellow et al., 2016) . 

LOSS FUNCTION AND CONSTRAINTS::

Consider the data set {u j (t i , •) : i, j = 0, 1, . . .}, where j indicates the j-th solution path with a certain initial condition of the unknown dynamics. We would like to train the PDE-Net with n δt-blocks. For a given n ≥ 1, every pair of the data {u j (t i , •), u j (t i+n , •)}, for each i and j, is a training sample, where u j (t i , •) is the input and u j (t i+n , •) is the label that we need to match with the output from the PDE-Net. We select the following simple 2 loss function for training:

L = i,j l ij , where l ij = ||u j (t i+n , •) -ũj (t i+n , •)|| 2 2 ,

where ũj (t i+n , •) is the output of the PDE-Net with u j (t i , •) as the input.

All the filters involved in the PDE-Net are properly constrained using their associated moment matrices. Let q 0 and q ij be the underlying filters of D 0 and D ij . We impose the following constrains (M (q 0 )) 1,1 = 1, (M (q 00 )) 1,1 = 1

and for i + j > 0 (M (q i,j )) k1,k2 = 0 k 1 + k 2 ≤ i + j + 2, (k 1 , k 2 ) = (i + 1, j + 1), (M (q i,j )) k1,k2 = 1 (k 1 , k 2 ) = (i + 1, j + 1).

For example, for 3 × 3 filters, we have

M (q 0 ) = M (q 00 ) = 1 and M (q 10 ) = 0 0 1 , M (q 01 ) = 0 1 0 , M (q 11 ) = 0 0 0 0 1 0 , . . . .

To demonstrate the necessity of learnable filters, we will compare the PDE-Net having the aforementioned constrains on the filters with the PDE-Net having frozen filters. To differentiate the two cases, we shall call the PDE-Net with frozen filters "the Frozen-PDE-Net".

To further increase the expressive power and flexibility of the PDE-Net, we may associate multiple filters to approximate a given differential operator. However, in order not to mess up the identifiability of the underlying PDE model, we may select only one of the filters to provide correct approximation to the given differential operator in the way as described above. The rest of the filters are constrained in the way that they only contribute to modify the local truncation errors. For example, consider two 3 × 3 filters {q 0 , q 1 } and constrain their moment matrices as follows

M (q 0 ) = 0 0 1 , M (q 1 ) = 0 0 0 .

Then, q 0 u + q 1 u is potentially a better approximation to u x (up to a constant) than q 0 u. However, for simplicity, we only use one filter to approximate a given differential operator in this paper.

NOVELTY OF THE PDE-NET::

Different from fixing numerical approximations of differentiations in advance in sparse regression methods (Schaeffer, 2017; Rudy et al., 2017) , using learnable filters makes the PDE-Net more flexible, and enables more robust approximation of unknown dynamics and longer time prediction (see numerical experiments in Section 3 and Section 4). Furthermore, the specific form of the response function F is also approximated from the data, rather than assumed to be known in advance (such as (Raissi & Karniadakis, 2017) ). On the other hand, by inflicting constrains on moment matrices, we can identify which differential operators are included in the underlying PDE which helps with identifying the nonlinear response function F . This grants transparency to the PDE-Net and the potential to reveal hidden physical laws. Therefore, the proposed PDE-Net is distinct from the existing learning based method to discover PDEs from data, as well as networks designed in deep learning for computer vision tasks.

2.3. INITIALIZATION AND TRAINING:

In the PDE-Net, parameters can be divided into three groups:

• filters to approximate differential operators;

• the parameters of the point-wise neural network to approximate F ;

• hyper-parameters, such as the number of filters, the size of filters, the number of layers, etc.

The parameters of the point-wise neural network are shared across the computation domain Ω, and are initialized by random sampling from a Gaussian distribution. For the filters, we initialize them by freezing them to their corresponding differential operators. For example, if a filter is to approximate ∂ ∂x , we freeze it by constraining its (1, 0)-moment to 1 and other moments to 0. During the training process, we release the filters by switching to the constrains described in Section 2.2.

Instead of training an n-layer PDE-Net directly, we adopt layer-wise training, which improves the training speed. To be more precise, we start with training the PDE-Net on the first δt-block, and then use the results of the first δt-block as the initialization and restart training the PDE-Net on the first two δt-blocks. Repeat until we complete all n blocks. Note that all the parameters in each of the δt-block are shared across layers. In addition, we add a warm-up step before the training of the first δt-block. The warm-up step is to obtain a good initial guess of the parameters of the point-wise neural network that approximates F by using frozen filters.

2.4. RELATIONS TO SOME EXISTING NETWORKS:

In recent years, a variety of deep neural networks have been introduced with great success in computer vision. The structure of the proposed PDE-Net is similar to some existing networks such as the Network-In-Network (NIN) (Lin et al., 2013) and the deep Residual Neural Network (ResNet) (He et al., 2016a; b) .

The NIN is an improvement over the traditional convolutional neural networks. One of the special designs of NIN is the use of multilayer perceptron convolution (mlpconv) layers instead of the ordinary convolution layers. An mlpconv layer contains the convolutions and small point-wise neural networks. Such design can improve the ability of the network to extract nonlinear features from shallow layers. The inner structure of one δt-block of the PDE-Net is similar to the mlpconv layer, and the multiple δt-blocks structure is similar to the NIN structure, except for the pooling and ReLU operations.

On the other hand, each δt-block of the PDE-Net has two paths (see Figure 1 and Figure 2 ): one is for the averaged quantity of u and the other is for the increment F . This structure coincides with the "residual block" introduced in ResNet. In fact, there has been a substantial study on the relation between ResNet and dynamical systems recently (E, 2017; Haber & Ruthotto, 2017; Sonoda & Murata, 2017) .

3. NUMERICAL STUDIES: CONVECTION-DIFFUSION EQUATIONS:

Convection-diffusion equations are classical PDEs that are used to describe physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: diffusion and convection (Chandrasekhar, 1943) . Convection-diffusion equations are widely applied in many scientific areas and industrial fields, such as pollutants dispersion in rivers or atmosphere, solute transferring in a porous medium, and oil reservoir simulation. In practical situations, usually the physical and chemical properties on different locations cannot be the same (called anisotropy in physics), thus it is more reasonable that convection coefficients and diffusion coefficients are variables instead of constants.

3.1. SIMULATED DATA, TRAINING AND TESTING:

We consider a 2-dimensional linear variable-coefficient convection-diffusion equation on The computation domain Ω is discretized using a 50 × 50 regular mesh. Data is generated by solving problem (8) using a high precision numerical scheme with pseudo-spectral method for spatial discretization and 4th order Runge-Kutta for temporal discretization (with time step size δt = 0.01).

EQUATION

We assume periodic boundary condition and the initial value u 0 (x, y) is generated from

EQUATION

where N = 9, λ k,l , γ k,l ∼ N (0, 1 50 ), and k and l are chosen randomly. In order to mimic real world scenarios, we add noise to the generated data. For each sample sequence u(x, y, t), t ∈ [0, 0.2], the noise is added as u(x, y, t) = u(x, y, t) + 0.01 × M W (10)

where M = max x,y,t {u(x, y, t)}, W ∼ N (0, 1) and N (0, 1) represents the standard normal distribution.

• The PDE-Net with 7 × 7 filters significantly outperforms the PDE-Net with 5 × 5 filters in terms of the length of reliable predictions (see Figure 3 and 4 ). To reach an O(1) error, the length of prediction for the PDE-Net with 7 × 7 filters is about 10 times of that for the PDE-Net with 5 × 5 filters. 

DISCOVERING THE HIDDEN EQUATION:

For the linear problem, identifying the PDE amounts to finding the coefficients {c ij : i + j ≤ 4} that approximate {f ij : i + j ≤ 4}. The coefficients {c ij : i + j ≤ 2} of the trained PDE-Net are shown in Figure 6 . Note that {f 11 } ∪ {f ij : 2 < i + j ≤ 4} are absent from the PDE (8), and the corresponding coefficients learned by the PDE-Net are indeed close to zero. In order to have a more concise demonstration of the results, we only show the image of {c ij : i + j ≤ 2} in Figure 6 .

Comparing the first three rows of Figure 6 , the coefficients {c ij } learned by the PDE-Net are close to the true coefficients {f ij } except for some oscillations due to the presence of noise in the training data. Furthermore, the last row of Figure 6 indicates that having multiple δt-blocks helps with estimation of the coefficients. However, having larger filters does not seem to improve the learning of the coefficients, though it helps tremendously in prolonging predictions of the PDE-Net. 

FURTHER EXPERIMENTS:

To further demonstrate how well the learned PDE-Net generalizes, we generate initial values following (9) with highest frequency equal to 12, followed by adding noise (10). Note that the maximum allowable frequency in the training set is 9. The results of long-time prediction and the estimated dynamics are shown in Figure 7 . Although oscillations are observed in the prediction, the estimated dynamic still captures the main pattern of the true dynamic.

The PDE ( 8) is of second order. In our previous experiments, we assumed that the PDE does not exceed the 4th order. If we know that the PDE is of second order, we will be able to have a more accurate estimation of the variable coefficients of the convection and diffusion terms. However, the prediction errors are slightly higher since we have fewer trainable parameters. Nonetheless, since we are using a more accurate prior knowledge on the unknown PDE, the variance of the prediction errors are smaller than before. These results are summarized in Figure 8 (green curves) and Figure 9 .

To further demonstrate the importance of the moment constraints on the filters in the PDE-Net, we trained the network without any moment constraints and skipped any steps that utilize the knowledge of the relation between the filters and differential operators (i.e. we skipped warm-up and the initialization using finite difference filters). For simplicity, we call the PDE-Net train in this way as the Freed-PDE-Net. The prediction errors of the Freed-PDE-Net are shown as the red curves in Figure 8 . Since without moment constraints, we do not know the correspondence of the filters with differential operators. Therefore, we cannot identify the correspondence of the learned variable coefficients either. We plot all the 15 variable coefficients (assuming the underlying PDE is of order ≤ 4) in Figure 10 . As one can see that the Freed-PDE-Net is better in prediction than the PDE-Net since it has more trainable parameters than the PDE-Net. However, we are unable to identify the PDE from the Free-PDE-Net. 

4. NUMERICAL STUDIES: DIFFUSION EQUATIONS WITH NONLINEAR SOURCE:

When modeling physical processes like particle transportation or energy transfer, in addition to convection and diffusion, we have to consider source/sink terms. In some problems, the source/sink plays an important role. For example, when convection-diffusion equations are used to describe the distribution and flow of pollutants in water or atmosphere, identifying the intensity of pollution source is equivalent to finding the source term, which is important for environmental pollution control problems. 

EQUATION

where c = 0.3 and f s (u) = 15 sin(u). The computation domain Ω is discretized using a 50 × 50 regular mesh. Data is generated by solving problem (11) using forward Euler for temporal discretization (with time step size δt = 0.0009) and central differencing for spatial discretization on 100 × 100 mesh, and then restricted to the 50 × 50 mesh. We assume zero boundary condition and Under review as a conference paper at 2018 the initial value u 0 (x, y) is generated by u 0 (x, y) = u 0 (x, y) x(2π-x)y(2π-y)

(2π) 4

, where u 0 is obtained from (9) with maximum allowable frequency N = 6. Same as the numerical setting in Section 3, Gaussian noise is added to each sample sequence u(x, y, t), t ∈ [0, 0.2] as described by (10).

Suppose we know a priori that the underlying PDE is a convection-diffusion equation of order no more than 2 with a nonlinear source depending on the variable u. Then, the response function F takes the following form

F = 1≤i+j≤2 f ij (x, y) ∂ i+j u ∂x i ∂y j + f s (u).

Each δt-block of the PDE-Net can be written as

ũ(t n+1 , •) = D 0 u(t n , •) + δt • (c 01 D 01 u + c 10 D 10 u + c 11 D 11 u + c 20 D 20 u + c 02 D 02 u) + fs (u),

where {D 0 , D ij : 1 ≤ i + j ≤ 2} are convolution operators and {c ij : 1 ≤ i + j ≤ 2} are 2D arrays which approximate functions f ij (x, y) on Ω. The approximation is achieved using piecewise quadratic polynomial interpolation with smooth transitions at the boundaries of each piece. The approximation of fs is obtained by piecewise 4th order polynomial approximation over a regular grid of the interval [-30, 30] with 40 grid points. The training and testing strategy is exactly the same as in Section 3. In our experiments, the size of the filters is 7 × 7. The total number of trainable parameters for each δt-block is approximately 1.2k.

4.2. RESULTS AND DISCUSSIONS:

This section presents numerical results of the trained PDE-Net using the data set described in Section 4.1. We will observe how the trained PDE-Net performs in terms of prediction of dynamical behavior and identification of the underlying PDE model.

PREDICTING LONG-TIME DYNAMICS:

We demonstrate the ability of the trained PDE-Net in prediction, which in the language of machine learning is the ability to generalize. The testing method is exactly the same as the method described in Section 3. Comparisons between PDE-Net and Frozen-PDE-Net are shown in Figure 11 , where we can clearly see the advantage of learning the filters. Long-time predictions of the PDE-Net is shown in Figure 12 and we visualize the predicted dynamics in Figure 13 . To further demonstrate how well the learned PDE-Net generalizes, we generate initial values following (9) with highest frequency equal to 10, followed by adding noise (10). Note that the maximum allowable frequency in the training set is only 6. The results of long-time prediction and the estimated dynamics are shown in Figure 14 . All these results show that the learned PDE-Net performs well in prediction. 

5. CONCLUSION AND DISCUSSION:

In this paper, we designed a deep feed-forward network, called the PDE-Net, to discover the hidden PDE model from the observed dynamics and to predict the dynamical behavior. The PDE-Net consists of two major components which are jointly trained: to approximate differential operations by convolutions with properly constrained filters, to approximate the nonlinear response by deep neural networks or other machine learning methods. The PDE-Net is suitable for learning PDEs as general as in (1). However, if we have a prior knowledge on the form of the response function F , we can easily adjust the network architecture by taking advantage of the additional information. This may simplify the training and improve the results. As an example, we considered a linear variable-coefficient convection-diffusion equation. The results show that the PDE-Net can uncover the hidden equation of the observed dynamics, and predict the dynamical behavior for a relatively long time, even in a noisy environment. Furthermore, having deep structure (i.e. multiple δt-blocks) and larger learnable filters can improve the PDE-Net in terms of stability and can prolong reliable predictions. As part of the future work, we will try the proposed framework on real data sets. One of the important directions is to uncover hidden variables which cannot be measured by sensors directly, such as in data assimilation. Another interesting direction which is worth exploring is to learn stable and consistent numerical schemes for a given PDE model based on the architecture of the PDE-Net. 

:



Figures and Tables:

Figure 1: The schematic diagram of a δt-block.

Figure 2: The schematic diagram of the PDE-Net: multiple δt-blocks.

where a(x, y) = 0.5(cos(y) + x(2π -x) sin(x)) + 0.6, b(x, y) = 2(cos(y) + sin(x)) + 0.8, c = 0.2 and d = 0.3.

Figure 3: Prediction errors of the PDE-Net (orange) and Frozen-PDE-Net (blue) with 5 × 5 (first row) and 7 × 7 (second row) filters. In each plot, the horizontal axis indicates the time of prediction in the interval (0, 60 × δt] = (0, 0.6], and the vertical axis shows the normalized errors. The banded curves indicate the 25% & 75% percentile of the normalized errors among 560 test samples.

Figure 4: Long-time prediction for the PDE-Net with 7 × 7 filters. The horizontal axis ranges in (0, 5]. Time step δt = 0.01.

Figure 5: Images of the true dynamics and the predicted dynamics. The first row shows the images of the true dynamics. The second row shows the images of the predicted dynamics using the PDE-Net having 3 δt-blocks with 5 × 5 and 7 × 7 filters. Time step δt = 0.01.

Figure 6: First row: the true coefficients of the equation. From the left to right are coefficients of u, u x , u y , u xx , u xy and u yy . Second row: the learned coefficients by the PDE-Net with 6 δt-blocks and 5 × 5 filters. Third row: the learned coefficients by the PDE-Net with 6 δt-blocks and 7 × 7 filters. Last row: the errors between true and learned coefficients v.s. number of δt-blocks (1, 2, . . . , 13) with different sizes of filters (blue for 5 × 5 and orange for 7 × 7).

Figure 7: Testing with higher frequency initializations (linear convection-diffusion equation). First row: long-time prediction. Second row: estimated dynamics. Here, δt = 0.01.

Figure 8: Prediction errors of the PDE-Net assuming the underlying PDE has order ≤ 4 (orange), order ≤ 2 (green) and Freed-PDE-Net (red) with 7 × 7 filters. In each plot, the horizontal axis indicates the time of prediction in the interval (0, 80 × δt] = (0, 0.8], and the vertical axis shows the normalized errors. The banded curves indicate the 25% & 75% percentile of the normalized errors among 560 test samples.

Figure 9: First row: the true coefficients of the equation. From the left to right are coefficients of u, u x , u y , u xx , u xy and u yy . Second row: the learned coefficients by the PDE-Net assuming the order of the PDE is ≤ 4 (same as the third row of Figure6). Third row: the learned coefficients by the PDE-Net assuming the order of the PDE is ≤ 2. Last row: the errors between true and learned coefficients v.s. number of δt-blocks (1, 2, . . . , 13) for PDE-Net assuming the PDE is of order ≤ 4 (orange) and ≤ 2 (green).

Figure 10: The images of all the variable coefficients learned from the Freed-PDE-Net.

Figure 11: Prediction errors of the PDE-Net (orange) and Frozen-PDE-Net (blue) with 7 × 7 filters. In each plot, the horizontal axis indicates the time of prediction in the interval (0, 0.6], and the vertical axis shows the normalized errors. The banded curves indicate the 25% & 75% percentile of the normalized errors among 560 test samples.

Figure 12: Long-time prediction for PDE-Net with 7 × 7 filters in (0, 2].

Figure 13: Images of the true dynamics and the predicted dynamics. The first row shows the images of the true dynamics. The second row shows the images of the predicted dynamics using the PDE-Net having 3 δt-blocks with 7 × 7 filters. Here, δt = 0.01.

Figure 14: Testing with higher (diffusion equation with a nonlinear source). First row: long-time prediction. Second row: estimated dynamics.Here, δt = 0.01.

Figure 15: First row: the true coefficients {f ij : 1 ≤ i + j ≤ 2} of the equation. Second row: the learned coefficients {c ij : 1 ≤ i + j ≤ 2} by the PDE-Net with 3 δt-blocks and 7 × 7 filters.

Peer review (iclr2018):
 The paper explores the use of deep learning machinery for the purpose of identifying dynamical systems specified by PDEs. The paper advocates the following approach: One assumes a dynamic PDE system involving differential operators up to a given order. Each differential operator term is approximated by a filter whose values are properly constrained so as to correspond to finite difference approximations of the corresponding term. The paper discusses the underlying wavelet-related theory in detail. A certain form of the dynamic function and/or source terms is also assumed. An explicit Euler scheme is adopted for time discretization. The parameters of the system are learned by minimizing the approximation error at each timestep. In the experiments reported in the paper the reference signal is provided by numerical simulation of a ground truth system and the authors compare the prediction quality of different versions of their system (eg, for different kernel size). Overall I find the paper good, well written and motivated. The advocated approach should be appealing for scientific applications of deep learning where not only the quality of approximation but also the interpretability of the identified model is important. Some suggestions for improvement: * The paper doesn*t discuss the spatial boundary conditions. Please clarify this. * The paper adopts a hybrid approach that lies in between the classic fully analytical PDE approach and the data driven machine learning. I would like to see a couple more experiments comparing the proposed approach with those extremes. (1) in the first experiment, the underlying model order is 2 but the experiment allows filters up to order 4. Can you please report if generalization quality improves if the correct order 2 is specified? (2) On the other side, what happens if no sum (vanishing order) constraints are enforced during model training? This abandons the interpretability of the model as approximating a PDE of given order but I am curious to see what is the generalization error of this less constrained system. Nit: bibtex error in Weinan (2017) makes the paper appear as E (2017).