Year: 2021

Abstract:
There have been growing interests in leveraging experimental measurements to discover the underlying partial differential equations (PDEs) that govern complex physical phenomena. Although past research attempts have achieved great success in data-driven PDE discovery, the robustness of the existing methods cannot be guaranteed when dealing with low-quality measurement data. To overcome this challenge, we propose a novel physics-encoded discrete learning framework for discovering spatiotemporal PDEs from scarce and noisy data. The general idea is to (1) firstly introduce a novel deep convolutional-recurrent network, which can encode prior physics knowledge (e.g., known PDE terms, assumed PDE structure, initial/boundary conditions, etc.) while remaining flexible on representation capability, to accurately reconstruct high-fidelity data, and (2) perform sparse regression with the reconstructed data to identify the explicit form of the governing PDEs. We validate our method on three nonlinear PDE systems. The effectiveness and superiority of the proposed method over baseline models are demonstrated.

heading: 1 INTRODUCTION
text: Deriving physical laws remains critical for understanding the dynamical patterns and complex spatiotemporal behaviors in nature. Well formulated governing equations facilitate scientific discovery and promote establishment of new disciplines. Traditionally, scientists seek to find governing equations by rigorously following the first principles, including the conservation laws, geometric brownian motion assumptions and knowledge-based deductions. However, these classical approaches show limited capability and slow progress in identifying clear analytical governing equations, and leave numerous complex systems under-explored. Thanks to the recent achievements in machine learning (ML), advances in governing equation/law discovery has been drastically accelerated.
The pioneer works (Bongard & Lipson, 2007; Schmidt & Lipson, 2009) apply symbolic regression to reveal the underlying differential equations that govern nonlinear dynamical systems without any prior background knowledge. Although this inspiring investigation implies the dawn of discovering the fundamental principles automatically in scientific ML, the scalability and overfitting issues herein have been the critical concerns. Recently, a ground-breaking work by (Brunton et al., 2016) introduces the sparse regression into the discovery of parsimonious governing equations (aka., Sparse Identification of Nonlinear Dynamics (SINDy)), based on the assumption that the dynamical systems are essentially controlled by only few dominant terms, which evolves to promote the sparsity in the candidate functions. Furthermore, the introduction of spatial derivative terms and the improvement of the sparsity-promoting algorithm as Sequential Threshold Ridge regression (STRidge) make SINDy applicable to general PDE discovery (i.e., the PDE-FIND algorithm) (Rudy et al., 2017). However, due to the demanding requirement of derivative estimation via numerical differentiation, these sparse regression methods largely rely on high-quality measurement data, which is barely accessible in the real-world sensing environment. Therefore, there is an urgent need to address the issue of discovering the underlying governing equations from low-resolution (LR) and noisy measurement data. Very recently, many works attempt to relieve the pain by reconstructing the high-resolution (HR) scientific data from LR measurements and then discovering the governing laws, including robust low-rank discovery (Li et al., 2020a), physics-informed spline learning (Sun et al., 2021) and physics-informed deep learning (Chen et al., 2021). Nevertheless, the dominant reconstruction part lacks explicit physical interpretability (i.e., grey-box model), which may induce the difficulty in optimization especially for high-dimensional PDE systems with delicate patterns.
Contributions. To overcome the challenges resulted from the noisy and LR measurements, we propose a novel physics-encoded deep learning framework (i.e., white-box model) for discovering the governing equations of physical systems. The main contributions of this paper are three-fold:
• We design a novel recurrent network architecture that is able to fit spatiotemporal measurement data for HR metadata generation with high accuracy and physical consistency. Notably, this network is characterized with capability to encode given physics knowledge (e.g., known terms, assumed PDE structure based on multiplicative feature operation, initial/boundary conditions) so that the given physics as model prior is respected rigorously.
• Based on the proposed reconstruction network, we develop a hybrid optimization approach to discover spatiotemporal PDEs from noisy and LR measurements. Specifically, to determine the explicit form of underlying governing equations, sparse regression is performed by using HR prediction inferred from the model. Finally, through inheriting the PDE terms and coefficients from the sparse discovery, we conduct the fine-tuning step with the completely physics-based network to further improve the accuracy of the scientific discovery.
• The proposed coupled framework establishes the physics-encoded data-driven predictive model from the noisy and LR measurement. To evaluate the performance of the model, we perform three numerical experiments on a variety of nonlinear systems. It is found that the proposed approach demonstrates excellent robustness and accuracy in spite of the high level of noise and the scarcity (spatially and temporally) in the measurements.

heading: 2 RELATED WORK
text: Dynamical system modeling. Simulating multi-dimensional PDE systems with neural networks (NNs) has been a renewed research topic, which can be dated back to last century (Lee & Kang, 1990; Lagaris et al., 1998). Like the traditional simulation methods, the recent breakthrough in NNbased modeling for physical systems also fall into two streams: continuous and discrete learning. The continuous learning scheme can avoid the notorious demand of time-step and grid sizes for numerical stability. Among the meshfree strategies, based on the incorporation of prior knowledge of PDE systems, there exist two directions seeking for the approximation of simulation with DNNs, i.e., the pure data-driven approaches (Han et al., 2018; Long et al., 2018; Wang et al., 2020) and physics-informed learning (Yu et al., 2017; Raissi et al., 2019; Rao et al., 2021b; Karniadakis et al., 2021). The former option relies on large amounts of high-quality data, while the physics-informed ML only requires scarce or even no labeled data due to the enhancement from physical constraints. Besides, it is worthwhile to mention that the recent studies on neural operators (Li et al., 2020b; Lu et al., 2021; Patel et al., 2021) have exhibited the great potential of learning the nonlinear, meshfree and infinite-dimensional mapping with DNNs for physical systems.
Regarding the discrete learning scheme, the mesh-based methods play a dominant role in simulating complex high-dimensional PDE systems. Herein, the convolutional architectures are mainly applied for handling regular grids in steady-state problems (Zhu & Zabaras, 2018; Zhu et al., 2019) and spatiotemporal systems (Bar-Sinai et al., 2019; Geneva & Zabaras, 2020; Kochkov et al., 2021; Ren et al., 2021), as well as irregular domains by considering elliptic coordinate mapping (Gao et al., 2021a). Moreover, with the theoretical development of graph neural networks (GNNs) (Kipf & Welling, 2016; Bronstein et al., 2017; Battaglia et al., 2018), researchers seek for a new direction for geometry-adaptative learning of nonlinear PDEs with arbitrary domains (Belbute-Peres et al., 2020; Sanchez-Gonzalez et al., 2020a; Gao et al., 2021c). In addition to the common mesh-based approaches, the particle-based frameworks have been attracting a lot of interest in scientific modeling for physical systems, especially by employing GNNs (Li et al., 2018; Ummenhofer et al., 2019; Sanchez-Gonzalez et al., 2020b).
Scientific data augmentation. Scientific data are typically sparse, noisy and incomplete. Therefore, data augmentation has been an active research area and received massive attention in the scientific ML community, e.g., data reconstruction from sparse measurements (Raissi, 2018; Callaham et al., 2019; Erichson et al., 2020; Rao et al., 2021a), super-resolution (SR) (Stengel et al., 2020; Gao et al., 2021b) and denoising (Fathi et al., 2018). Generally, based on the type of incorporated data, we categorize the existing data augmentation methods into two groups. The first one is to extract the coherent structure and correlation features as the basis by utilizing the sufficient high-
quality datasets. For instance, in the context of SR, the majority of NN-based frameworks (Liu et al., 2020; Esmaeilzadeh et al., 2020; Fukami et al., 2021) aim to reconstruct the high-resolution (HR) full-field data from low-resolution (LR) measurements by employing HR ground truth reference as the constraint. Nevertheless, the second type of data augmentation takes advantage of the physicsbased model and can accurately capture the high-fidelity dynamical patterns from the LR and noisy measurement data. The recent advances of physics-based learning models have shown the the great success of data augmentation for PDE systems with only sparse and noisy data (Raissi, 2018; Rao et al., 2021a) and even without any labeled data (i.e., only with PDEs) (Gao et al., 2021b).
Data-driven discovery. The data-driven discovery is also a renascent research attempt. The earliest work (Dzeroski & Todorovski, 1995) introduced the inductive logic programming to find the natural laws from experimental data. With the rapid advance of machine intelligence, the past three decades have seen a considerable amount of literatures growing up around the interplay between ML and scientific discovery. Essentially, two genres have been observed in distilling governing equations, with respect to the searching of candidate functions. Firstly, the natural and optimal solution to automatically identify the governing equations for dynamical systems is to learn a symbolic model from experimental data. Hence, symbolic regression (Schmidt & Lipson, 2009) and symbolic neural networks (SNN) (Sahoo et al., 2018; Kim et al., 2020) have been explored for inferring the physical models without any prior knowledge. In addition, GNN (Cranmer et al., 2020) has also been studied to identify the nontrivial relation due to its excellent inductive capability.
The second group manually defines a large library functions based on the prior knowledge, and utilizes the sparse-prompting algorithms to discover the parsimonious physics model efficiently. The most representative works are SINDy (Brunton et al., 2016) and PDE-FIND (Rudy et al., 2017) for ordinary differential equations (ODEs) and PDEs, respectively. Nevertheless, the standard sparse representation-based methods are limited to the high-fidelity noiseless measurement, which is usually difficult and expensive to obtain. Hence, many recent efforts have been devoted to discover PDEs from sparse/noisy data (Li et al., 2020a; Xu et al., 2021; Sun et al., 2021; Chen et al., 2021).

heading: 3 METHODOLOGY

heading: 3.1 PROBLEM DESCRIPTION
text: Discovering the explicit form of the governing equation(s) for a given system from the measurement data is of critical importance for scientists to understand some underexplored processes. To formalize this problem, which is also known as the data-driven equation discovery, let us assume the response of a multi-dimensional spatiotemporal system is described by
ut = F (x, t,u,u,ux,uy,∇u, · · · ) (1) where x = [x, y]T and t denote the spatial and temporal coordinates respectively, u denotes the state variable possibly consisting of multiple components (e.g., u has two components u and v) and the subscript denotes the partial derivatives, e.g., ux = ∂u/∂x. F is a function parameterizing the system dynamics with possible linear/nonlinear terms (e.g., u2, ux and ∇u, where the Nabla operator ∇ is defined as [∂/∂x, ∂/∂y]T). The objective of the data-driven equation discovery is to uncover the closed form of F provided time series measurements on some fixed spatial locations. Although this problem has been studied extensively over the past few years, the combinatorically large search space of the possible equation appears as one of the major obstacles that hinders the discovery of the governing equation from very limited measurements. The existence of noise, as well as the coarse-sampling of the measurement, also adds to the difficulty of this problem. To address these challenges, we proposed a novel physics-encoded recurrent network that is able to reconstruct the high-fidelity solution from the noisy and LR measurements with the help of prior physics knowledge. We successfully combine this network with the sparse regression for discovering the explicit form of the governing PDE of dynamical systems.

heading: 3.2 PHYSICS-ENCODED NETWORK FOR DATA RECONSTRUCTION
text: Obtaining the numerical derivatives for constructing the library matrix is crucial to the equation discovery with sparse regression. Due to the LR and the noise in measurements, directly applying
the traditional finite difference (FD) onto the raw measurements may lead to inaccurate derivatives, which would significantly complicate the equation discovery. To handle this issue, we introduce a recurrent network architecture, namely the Physics-encoded Recurrent Convolutional Neural Network (PeRCNN) shown in Fig. 1(1a) and (1b), to establish the data-driven predictive model from LR noisy measurements. PeRCNN consists of a fully convolutional (Conv) network as initial state generator (ISG) and a well designed Conv block, namely Π-block (product) for updating the state variables on high-resolution grid recurrently. Note that our network does not require any HR data as the “label” for training, which is fundamentally different from the classic super resolution task.
The network architecture mimics the forward Euler time stepping Û (k+1) = Û (k) + F(Û (k);θ)δt, where Û
(k) ∈ RH×W denotes the HR prediction on a H × W grid at kth time step, δt is the time spacing and θ is a set of trainable variables for approximating F . Since we seek the predictive model for HR solution, the ISG component is introduced to generate a HR initial state from the noisy and LR measurement ũ(0) for the recurrent computation. The recurrent Π-block is a major innovation of this network architecture to capture the dynamics of F . As shown in Fig. 1(1b), the input state variable to Π-block would firstly be mapped to the hidden space through multiple parallel Conv layers, whose output is then fused via the element-wise product operation. A Conv layer with 1 × 1 kernel (Lin et al., 2013) is subsequently used to linearly combine multiple channels into the output (i.e., approximated F). Mathematically, the Π-block seeks to approximate the function F via polynomial combination of solution u and its spatial derivatives, given by
F̂(Û (k)) = Nc∑ c=1
fc · (
Nl∏ l=1
K(c,l) ⊛ Û (k+1) )
(2)
where Nc and Nl denote the numbers of channels and parallel Conv layers respectively; ⊛ denotes the Conv operation; (c, l) indicates the Conv filter K of l-th layer and c-th channel; fc is the weight in 1 × 1 Conv layer while the bias is omitted for simplicity. This multiplicative representation of Π-block promotes the network expressiveness for nonlinear function compared with the additive representation (i.e., weighted sum of multiple output) commonly seen, such as in PDE-Net (Long et al., 2018). We demonstrate in Appendix Section A that Π-block is an universal polynomial approximator to the nonlinear function F . In addition, this network architecture features the capability to encode physics knowledge, such as the existing term in the PDE or the known initial and boundary conditions (I/BCs). As shown in Fig. 1(1b), a highway physics-based Conv layer1 can be created to account for existing terms known a priori in the PDEs. In Appendix Section D, we demonstrate that such a highway connection can accelerate the training speed and improve the accuracy of the fitted model. In addition, the knowledge of I/BCs of the system can be encoded into the network through the customized padding in Conv operation. For example, the periodic boundary condition (PBC) could be encoded into the model via circular padding on the prediction (e.g., [1,2,3,4] padded to [4,1,2,3,4,1]). Note that “physics-encoded” has the meaning of threefold: (1) PDE solution time marching that follows the recurrence of forward Euler integration, (2) encoding a priori PDE terms into the Π-block, and (3) encoding multiplicative polynomial-type PDE terms into the Π-block.
1This Conv layer with the corresponding FD filter would be freezed throughout the training.

heading: 3.3 SPARSE DISCOVERY
text: Based on the reconstructed HR solution (generated metadata), we further employ the sparse discovery algorithm (Rudy et al., 2017) to distill the underlying governing equations. The sparse regression method roots on a critical observation that the right hand side (RHS) of Eq. (1) for the majority of natural systems consists of only a few terms, i.e., the F demonstrates sparsity in the space of possible functions. Given a library of candidate functions of F , we can formulate the discovery of the governing equation as a regression problem. Let us consider the state variable with one single component u ∈ Rns×nt , which is defined on ns spatial locations and within nt time steps. After flattening the state variable into a column vector U ∈ Rns·nt×1, a library Θ(U) ∈ Rns·nt×s that encompasses a pool of s candidate functions (e.g., linear, nonlinear terms, derivatives, etc.) can be constructed to represent the governing equations of physical systems. Each column in Θ(U) represents a candidate function in F , such as Θ(U) = [1,U,U2, . . . ,Ux,Uy, . . . ]. Assuming that Θ(U) has a sufficiently rich column space, the governing equation described by Eq. (1) can be rewritten as a linear system by using a sparse coefficient vector Ξ ∈ Rs×1, namely,
Ut = Θ(U)Ξ. (3)
To compute a reasonable Ξ such that its sparsity is satisfied while the regression error is small, the Sequential Threshold Ridge regression (STRidge) algorithm (Rudy et al., 2017), among other effective sparsity-promoting methods such as the Iterative Hard Thresholding (IHT) method (Haupt & Nowak, 2006; Blumensath & Davies, 2009), is used in this paper due to its superior performance compared with other sparsity-promoting algorithms, such as LASSO (Tibshirani, 1996) and Sequentially Thresholded Least Squares (STLS) (Brunton et al., 2016). For a given tolerance that filters the entries of Ξ with small value, we can obtain a sparse representation of F . Iterative search with STRidge can be performed to find the optimal tolerance according to the selection criteria:
Ξ∗ = argmin Ξ {||Ut −Θ(U)Ξ||2 + γ||Ξ||0} (4)
where ||Ξ||0 is used to measure the sparsity of coefficient vector; ||Ut − Θ(U)Ξ||2 denotes the regression error; γ is the coefficient that balances the sparsity and the regression errors. Since the optimization objective has two components, we apply Pareto analysis to select an appropriate γ (see Appendix Section H). This procedure to find a sparse Ξ is also known as sparse regression.

heading: 3.4 HYBRID PEDL FRAMEWORK FOR EQUATION DISCOVERY
text: In this part, we present the hybrid/coupled scheme of the proposed PeRCNN and the sparse regression for discovering the governing PDEs. As shown in Fig. 1, the whole equation discovery process can be divided into three stages, including the metadata reconstruction, the sparse regression for PDE discovery and the fine-tuning of PDE coefficients. The proposed framework enables us to iterate over these three stages to keep refining the result until convergence. The purpose of each stage, as well as the technical details, is detailed as follows.
Data reconstruction. Since the experimental measurements are typically sparse and noisy, we firstly utilize the PeRCNN introduced in Section 3.2 to reconstruct the high-fidelity data in order to obtain accurate derivative terms in the stage of sparse regression. Given the measurement2 ũ ∈ Rn′t×H′×W ′ on a coarse grid (with H ′ ×W ′ resolution) at n′t time steps, we establish a data-driven model that minimizes the misfit error (i.e., mean squared error) between the predictions from the PeRCNN model and the input measurements. In specific, the loss function is given by
L(θ) = MSE ( Û(x̃)− ũ ) + ηMSE ( Û (0) − I(ũ(0)) ) . (5)
Note that the 1st part in RHS of Eq. (5) is adopted for training the network based on the LR data ũ. Since we use x̃ to denote the set of locations (i.e., coarse grid) on which the LR data is collected, Û(x̃) denotes the mapping of HR prediction Û ∈ Rnt×H×W on the coarse grid x̃, where nt is the number of time steps for reconstructed data and H × W the resolution for the fine grid. Besides, the 2nd part in RHS of Eq. (5) is the regularizer term referring to the discrepancy of measurementinterpolated HR initial state I(ũ(0)) and the predicted HR initial state Û (0) from ISG, where the
2ũ herein denotes the LR measurement while u in Eq. (1) represents the symbol of state variable.
superscript “0” indicates the first snapshot of the measurement or prediction. The introduction of the initial condition (IC) regularizer is because of its effectiveness in avoiding the model overfitting. Moreover, η is the weighting coefficient for balancing these two loss components in the model training. We first pretrain ISG to ensure the accuracy of the HR initial state (i.e., Û (0)
). Then the whole network (ISG and recurrent Π-block) is trained using the loss function of Eq. (5).
Sparse regression. With the reconstructed high-fidelity (i.e., HR and denoised) solution, we are able to reliably and accurately perform sparse regression for the explicit form (analytical structure) of PDEs. As described in Section 3.3, given a library of candidate functions Θ(U), sparse regression seeks to find a suitable coefficient vector Ξ that balances model complexity and accuracy. This is realized by solving the optimization problem described by Eq. (4) with the STRidge algorithm. Note that sparse regression is performed on subsampled Θ(U) (i.e., 10% rows are randomly sampled) to avoid a very large library matrix that would lead to extremely slow computation or exceeding computer’s memory limit. Similar strategy is adopted in PDE-FIND (Rudy et al., 2017).
Fine-tuning of coefficients. The obtained coefficients from sparse regression may not fully exploit all the available measurement as the regression is performed on subsampled HR data. To further improve the result of equation discovery, we present a fine-tuning step to produce the final explicit governing equation. In the fine-tuning step, all the spatiotemporal measurements are used to train a recurrent block completely based on the identified PDE structure from the sparse regression (see
Fig. 1(3b)). The coefficient of each term, which is treated as a trainable variable in the network, can be obtained by minimizing the misfit error between the predictions and measurements. In Appendix Section G, we show that fine-tuning can effectively improve the accuracy of the discovered PDE.

heading: 4 EXPERIMENTS

heading: 4.1 DATASETS
text: To examine the effectiveness of the proposed approach, we test on three different datasets (i.e., PDE systems) which cover the 2D Burgers’ equation, 2D Lambda-Omega (λ–Ω) and 2D Gray-Scott (GS) reaction-diffusion (RD) equations. All of them are valuable mathematically and practically, and usually work as the benchmark examples. For instance. 2D Burgers’ equation has wide applications in fluid mechanics, acoustics and traffic flow, while the GS RD equation is extensively used to model the process in chemistry and biochemistry (whose response has very complex patterns). Detailed description of each equation is provided in Appendix Section B.

heading: 4.2 EVALUATION METRICS
text: We adopt three quantitative metrics to examine the performance of our method for PDE discovery:
Relative ℓ2 error of coefficients. The relative ℓ2 error, defined as E = ||Ξid −Ξtrue||2/||Ξtrue||2, measures the relative distance between the identified coefficient vector Ξid and the ground truth Ξtrue. However, when the magnitude of coefficients vary significantly (e.g., the GS RD equation), E cannot reflect the result well as the small coefficients get overwhelmed. Hence, we introduce the following non-dimensional measures to gain a more insightful evaluation of the results.
Precision and recall. The PDE discovery can also be considered as a binary classification problem (e.g., whether a term exists or not) given a candidate library. Therefore, we present the precision and recall to evaluate the performance of the proposed method. The recall measures the percentage of the successfully identified coefficients among the true coefficients, defined as R = ||Ξid ⊙ Ξtrue||0/||Ξtrue||0 where ⊙ denotes element-wise product of two vectors3. Similarly, the precision has the definition of P = ||Ξid ⊙Ξtrue||0/||Ξid||0.

heading: 4.3 EXPERIMENT SETUP
text: Generation of measurements. The noisy and LR measurements are synthetized through numerical simulations with high-order FD methods. 9-point stencil is used for computing the spatial deriva-
3A successful identification occurs when the entries in both identified and true vectors are nonzero.
tives while time marching is performed through Runge-Kutta scheme. For the 2D Burgers’ equation system, we select ν to be 0.005 (Geneva & Zabaras, 2020) and generate the solution on a unit square domain discretized by the 101 × 101 grid. 200 time steps are simulated for a total duration of 0.05 second. The numerical solution is downsampled to measurements as ũ ∈ R40×2×51×51. For the 2D λ–Ω RD system, β, µu and µv are set to be 1, 0.1 and 0.1 respectively. We discretized the 20 × 20 square domain into the 101× 101 grid. 200 time steps are considered for a physical duration of 2.5 seconds. The measurements are downsampled from the solution to be ũ ∈ R40×2×51×51. For the 2D GS RD system, the solution is generated on the 101 × 101 grid for a unit square domain using the parameters of µu = 2 × 10−5, µv = 5 × 10−6, κ = 0.06 and f = 0.04. 800 time steps are simulated for the total duration of 400 seconds. The measurements are obtained by downsampling the numerical solution to ũ ∈ R160×2×26×26. Periodic boundary condition is adopted for all three systems. In addition, to mimic the measurement noise in the real world, we add the Gaussian noise of a given level (e.g., 5% and 10%) to the numerical solution before downsampling for the LR data.
Baselines. In the experiments, we compare the performance of our proposed method with the widely used PDE-FIND (Rudy et al., 2017). Different from the settings in (Rudy et al., 2017), we assume only LR measurements are available for computing the derivatives in Θ(U). Specifically, the FD is used to obtain the derivatives for noiseless data. However, in the case that data is heavily corrupted by noise, the derivatives of the fitted Chebyshev polynomial are used. In addition, we also consider the baselines of the sparse regression coupled with a fully connected neural network (namely FCNN+SR) or PDE-Net4 (namely PDE-Net+SR). In these two baselines, the FCNN and PDE-Net (Long et al., 2018) are used respectively to fit the LR measurement and perform inference for HR data while SR is used to discover the PDE. Notably, automatic differentiation (Baydin et al., 2017) is used in FCNN+SR to compute the partial derivatives in Θ(U). However, like many other methods (Chen et al., 2021; Xu et al., 2021) that utilize the FCNN to fit the sparse data, FCNN+SR lacks the capability of encoding prior physics knowledge. Therefore, the prior knowledge of known diffusion terms would not be utilized in the data reconstruction phase of FCNN+SR.

heading: 4.4 RESULTS
text: 2D Burgers’ equation. We first test the proposed approach on the 2D Burgers’ equation. To reconstruct the high-fidelity data, we use the PeRCNN with 3 parallel Conv layers of 16 channels and filter size of 5. These hyperparameters of network architecture are selected through hold-out validation. The range for selecting the network hyperparameters are given in Appendix Section C. Furthermore, in this and the following examples, we assume the measurement exhibits the ubiquitous diffusion phenomenon. That said, the diffusion layer would be encoded into the network architecture. Also in the sparse regression stage, the coefficient of diffusion term would be exempted from being filtered. We train the network with Adam optimizer for 15,000 iterations. The learning rate is initialized to be 0.002 and decreases to 97% of the previous for every 200 iterations.
Once the training is completed, the HR solution (i.e., 201 snapshots of 101 × 101) can be inferred from the trained model. Figure 2(a) provides the snapshots of the reconstructed HR data under 5% noise from each method. It can be seen that our method and the FCNN has a much smaller reconstruction error, which would give rise to more accurate derivatives for constructing the library matrix. We also observe that PDE-Net struggles in reconstructing the HR data. Based on our experiments and the claim PDE-Net paper (Long et al., 2018) made, this is because the additive operation limits PDE-Net to only learning the response of linear PDEs. Therefore, unlike our model, PDE-Net lacks the capability to exactly express nonlinear terms like uux and u2v. To prepare for the subsequent sparse regression, we established a group of 70 candidate functions that consists of polynomial terms {1, u, v, u2, uv, v2, u3, u2v, uv2, v3}, derivatives {1, ux, uy, vx, vy,∆u,∆v} and their combinations. The selection of the crucial hyperparameter γ in sparse regression is conducted through Pareto analysis, which is detailed in Appendix Section H. In addition to the diffusion terms known a priori, the result of sparse regression indicates the presence of several other terms in F , i.e., Su = {uux, vuy} and Sv = {uvx, vvy} respectively. As the last step, the fine-tuning is performed using the network built completed based on the discovered PDE structure in the second
4As the original PDE-Net does not account for LR measurement data, we replace the Π-block in our model with the δt-block in PDE-Net to make it a HR predictive model.
stage (see Fig. 1(3b)), which gives us the final discovered equation:
ut = 5.0113× 10−3∆u− 1.004uux − 1.004vuy, vt = 4.9953× 10−3∆v − 1.009uvx − 1.002vvy.
(6)
The precision, recall and relative ℓ2 error of the discovered PDE under various noise levels are provided in Table 1. It can be observed that all methods except PDE-Net+SR perform well when the measurement is noise-free. Furthermore, the performance of baselines deteriorate when the noise level increases to 10% while our approach demonstrates better robustness against the noise. Our approach achieves 5.9× 10−3 relative ℓ2 error, 100% recall and precision under the 10% noise level. In the Appendix Sections E and F, we demonstrate that our method could handle even larger noise (up to 30%) and sparser data. We also consider a more challenging case (increasing the Reynolds number to 500 for the Burgers’ equation), in which our method still work well for both response reconstruction and PDE discovery (see Appendix Section I). We also show the interpretability of the trained Π-block network for the 2D Burgers’ case (see Appendix Section J).
2D λ–Ω RD equation. In this example, the recurrent network employed for data reconstruction has 3 parallel Conv layers, 16 channels and 5 × 5 kernels. The same procedure in previous example for training the reconstruction network is adopted here. Once the training is completed, we infer from the network the HR solution at finer spatiotemporal grid. The reconstruction error of our method and the PDE-FIND is shown in Fig. 2(b). The same candidate set in 2D Burgers’ example (with 70 candidate functions) is adopted to establish the library matrix Θ(U). We randomly subsample 10% of the HR points for the sparse regression, whose result (5% noise case) gives the set of terms with nonzero coefficients, i.e., Su = {∆u, u, u3, u2v, uv2, v3} and Sv = {∆v, u, u3, u2v, uv2, v3}. With the set of existing terms in the governing PDE, the finetuning step is performed using all the available measurement. The obtained scalar coefficients from the fine-tuning renders us with the final governing equation, which reads as
ut =0.096∆u+ 1.038u− 1.048u3 + 1.004u2v − 1.050uv2 + 0.998v3, vt =0.100∆v + 1.014v − 0.998u3 − 1.025u2v − 0.999uv2 − 1.015v3.
(7)
The evaluation metrics computed from the discovered PDE under various noise levels are provided in Table 1. Our approach achieves 5.44 × 10−2 relative error, 100% recall and 91.6% precision under 10% noise. Our method is also capable of handling much sparser data (Appendix Section F).
2D GS RD equation. In the last example, we assume the measurements are some extremely LR (i.e., 26× 26) snapshots (see measurements in Fig. 2(c)). To remedy the lack of resolution in spatial dimension, we further introduce the assumption that reaction term is in the form of polynomial after reviewing the existing literatures. As the first step, we employ the network with 3 parallel Conv layers, 8 channels and filter size of 1 to reconstruct the data. From the reconstruction model, we infer the HR (i.e., 101× 101) solution at finer time instances. Snapshots of HR solution in Fig. 2(c) show our method is able to restore the high-fidelity data from the LR and noisy measurement very well. With the HR solution, a library of polynomials up to the third degree Θ(u, v) = [1, u, v, u2, uv, v2, u3, u2v, uv2, v3] is constructed for sparse discovery. We subsample 10% HR data for the sparse regression, which gives us the remaining terms Su = {∆u, 1, u, uv2} and Sv = {∆v, v, uv2}. With the existing terms in the PDE, we perform the fine-tuning using the completely physics-based Conv block. The final obtained PDEs in the case of 5% noise are
ut =2.001× 10−5∆u− 1.003uv2 − 0.04008u+ 0.04008, vt =5.042× 10−6∆v + 1.009uv2 − 0.1007v.
(8)
As presented in Table 1, our proposed approach performs well on discovering the PDE as a result of fully utilizing the prior physics knowledge and the powerful expressiveness of the model, while PDE-FIND , PDE-Net and FCNN+SR struggle due to the extremely LR measurements.

heading: 5 CONCLUSION
text: We propose a hybrid computational framework for discovering nonlinear PDE systems from sparse (LR) and noisy measurement data, which couples the novel network architecture PeRCNN for the data-driven modeling of dynamical systems and the sparse regression for distilling the dominant candidate functions and coefficients. First of all, one major innovation of our framework to scientific discovery comes from the nonlinearity achieved by the element-wise product operation and its powerful expressiveness (see Appendix Section A) for the nonlinear function F(·) in PDEs. Besides, our network also features the capability to encode prior physics knowledge via a highway Conv layer that helps to overcome the challenges brought by the scarcity and noise of the measurements. Secondly, we successfully marry PeRCNN to the sparse regression algorithm to solve the crucial equation discovery issues. The coupled scheme enables us to iteratively optimize the network parameters, and fine-tune the discovered PDE structures and coefficients. Through the numerical validation, we demonstrate the superior reconstruction capability and excellent discovery accuracy of our proposed framework against three baseline approaches. Overall, we provide an effective, interpretable and flexible approach to accurately and reliably discover the underlying physical laws from the imperfect and coarse-meshed measurements.

heading: A Π-BLOCK AS A UNIVERSAL POLYNOMIAL APPROXIMATOR
text: In the proposed PeRCNN, Π-block acts as an universal polynomial approximator to unknown nonlinear functions while the physics-based highway Conv layer (i.e., with FD stencil as the Conv filter) accounts for the prior knowledge on the governing equation. Notably, the Π-block achieves its nonlinearity through the element-wise product operation (see Eq. (2)) instead of what is widely adopted by the traditional Conv network, in which nonlinear layers are sequentially interwined with linear layers. Compared with the additive form representation F(u) =
∑ 0≤i+j≤N fij · (Kij ⊛ u) seen in
related work (Long et al., 2018; Guen & Thome, 2020), this multiplicative representation of Π-block promotes the network expressiveness for nonlinear functions. To support this claim, we need to recognize that Π-block roots on the numerical differentiation (i.e., forward Euler time stepping), which can guarantee the convergence of the solution under two conditions: (i) the RHS F of Eq. (1) can be computed accurately; and (ii) the time spacing δt is sufficiently small. As the second condition can be satisfied by selecting an appropriate δt, we would show that, by construction, our Π-block can approximate any nonlinear polynomial involving spatial derivatives, such as uux + vuy and u∆u. This can be easily argued by proving that each monomial (e.g., uux, vuy , etc.) can be approximated by multiplying the output of given number of Conv layers. To this end, we provided the following lemma and its proof reproduced from (Long et al., 2019).
Lemma 1: By construction, a convolutional filter K with a given size can approximate any differential operator with prescribed order of accuracy.
Proof 1: Without loss of generality, let us consider a bivariate differential operator L(·) w.r.t variables x and y, we have
L(u) = N−1 2∑ k1,k2=−N−12 K[k1, k2] N−1∑ i,j=0 ∂i+ju ∂ix∂jy ∣∣∣∣ (x,y) ki1k j 2 i!j! δxiδyj +O(|δx|N−1 + |δy|N−1)
=
N−1 2∑
k1,k2=−N−12
K[k1, k2]u(x+ k1δx, y + k2δy) +O(|δx|N−1 + |δy|N−1)
= K ⊛ u+O(|δx|N−1 + |δy|N−1)
(9)
where N is the size of the filter. The weight matrix K[·, ·] is indexed by k1 and k2. Letting the filter’s entry K[k1, k2] be the corresponding Taylor series coefficient, we can see the approximation error is bounded by O(|δx|N−1 + |δy|N−1). Therefore, compared with the traditional black-box models (e.g., deep neural networks) for representing nonlinear functions, Π-block possesses two main advantages:
• The nonlinear function F in the form of multivariate polynomial (e.g., u·∇u+u2v) covers a wide range of well-known dynamical systems, such as Navier-Stokes, reaction-diffusion (RD), Lorenz equations, to name only a few. Since the spatial derivatives can be computed by Conv kernels, a Π-block with n parallel Conv layers of appropriate filter size is able to represent a polynomial up to the nth order.
• The Π-block is flexible in representing the nonlinear function F . For example, a Π-block with 2 parallel layers of appropriate filter size ensembles a family of polynomials up to the 2nd order (e.g., u, ∆u, uv, u · ∇u), with no need to explicitly define the basis functions.

heading: B DETAILED DESCRIPTION OF THE CONSIDERED PDES
text: The first 2D Burgers’ equation studied in this paper has the form of ut + u · ∇u = ν∆u (10)
where ∇ is the Nabla operator; ∆ is the Laplacian operator; u = [u, v]T is the fluid velocities and ν is the viscosity coefficient. The second and third RD equation system can both be described by
ut = D∆u+R(u) (11)
where u = [u, v]T is the concentration vector; D = [µu, 0; 0, µv] is the diagonal diffusion coefficient matrix; R(u) denotes the reaction vector. For the λ–Ω RD system, the reaction vector has the form of R(u) = [(1 − u2 − v2)u + β(u2 + v2)v, − β(u2 + v2)u + (1 − u2 − v2)v]T where β is a reaction parameter. For the GS RD system, the nonlinear reaction vector reads R(u) = [−uv2 + f(1 − u), uv2 − (f + κ)v]T where κ and f denote the kill and feed rate respectively.

heading: C HYPERPARAMETER SELECTION OF PERCNN IN DATA RECONSTRUCTION
text: The selection of the hyperparameters of PeRCNN for data reconstruction is performed through the hold-out validation. Among the entire measurement, 10% data is split as the validation dataset. The range of hyperparameters we considered are summarized in Table C.1. The value in the parenthesis denotes the final parameter adopted in the data reconstruction.

heading: D ENCODE PHYSICS KNOWLEDGE TO ENHANCE THE MODEL
text: As introduced in Section 3.2, one salient characteristic of the proposed network is the capability to encode prior physics knowledge. To examine how the encoded physics enhances the accuracy of the data-driven model, here we use the 2D Burgers’ equation system as a testbed. In the experiment, we consider two different cases: (1) no prior knowledge is available, and (2) the diffusion term (∆u) is known to exist in F . The setting in the second case leads to an additional physics-based Conv layer in PeRCNN to account for the diffusion term. Apart from this, other hyperparameters are kept same in the experiment.
Figure D.1 shows the learning curves of these two data-driven models. It can be observed that the network encoded with the diffusion term exhibits faster convergence and better accuracy compared with the competitor. As a result of encoding prior physics knowledge, the PeRCNN does not need to learn a diffusion operator through the back propagation of error. That said, incorporating available physics can facilitate the learning efficiency and accuracy of the data-driven model.

heading: E TOLERANCE TO LARGER NOISE LEVEL
text: To examine the scalability of our method with regard to the Gaussian noise level, we performed additional experiments with the noise level of 20% and 30% for all three cases. The same settings from Section 4 is adopted in these experiments. The performance of our method on various noise levels is summarized in Table E.1. We observed that larger noise affects the accuracy of our method differently. For the 2D GS RD equation, the accuracy of our method deteriorates drastically as the noise level increases. That is because the measurement data used in this example has a relatively low resolution (i.e., 26 × 26). Even worse, the phenomenon in 2D GS RD system is characterized with local patterns and sharp gradients (see Fig. 2c). Increasing the Gaussian noise would further blur the spatial information in the measurement data. However, for the case of 2D λ-Ω RD and Burgers’ equation where the measurement data has a decent resolution (i.e., 51 × 51) and smoother spatial patterns, our method stays competitive even for 30% Gaussian noise. Overall, there is a trade-off between the sparsity and noise level of the measurement data.
Table E.1: Performance of the proposed method on larger noise level.
Cases Noise level /% Precision /% Recall /% Relative ℓ2 error /×10−2
Burgers’
0 100 100 0.50 5 100 100 0.54 10 100 100 0.59 20 75.0 100 7.33 30 75.0 100 9.10
λ–Ω RD
0 100 100 1.18 5 100 100 2.69 10 91.6 100 5.44 20 91.6 100 7.65 30 91.6 100 15.89
GS RD
0 100 100 1.59 5 100 100 2.85 10 85.7 85.7 10.03 20 46.2 85.7 86.05 30 33.3 37.5 238.31

heading: F MEASUREMENT DATA SPARSITY
text: Our method assumes the measurement data is collected on a fixed coarse sensor grid (e.g., 21× 21 sensors), captured either by point-wise sensor units or by imaging techniques. Although the ground truth data is of high fidelity, our measurement data is heavily spatiotemporally-downsampled and, as a result, remains sparse (e.g., a finite number of LR snapshots). As for the scalability with data sparsity, we performed numerical experiments on the 2D GS and λ–Ω RD equations using various spatial resolutions, namely, 51 × 51, 26 × 26, 21 × 21 and 11 × 11. The observation is that the tolerable sparsity of our method depends on the spatial pattern of the system. For example, in the 2D GS RD example, the lowest resolution of data to guarantee the discoverability of the governing PDE is 26 × 26 since the very complex maze-like pattern (see Fig. 2c) is in a relatively fine scale. However, since the solution in the λ–Ω RD example is much smoother and periodic (see Fig. 2b), our method is able to discover a major portion of the PDE (i.e., precision 90.9%, recall 83.3%) with the spatial resolution as low as 11× 11 (a very sparse data scenario).
G IMPORTANCE OF THE COEFFICIENT FINE TUNING
Our proposed approach for discovering PDEs consists of three stages including the data reconstruction, the sparse regression and the fine-tuning of coefficients. In fact, the sparse regression has
already been able to determine the explicit form of the discovered PDE. Then a natural question would be why the last fine-tuning step is necessary. As its name indicates, the fine-tuning step is able to further improve the accuracy of the discovered PDE, for mainly two reasons: (i) The sparse regression step gives an explicit PDE based on the primitive assumption of a sparse Ξ. However, the discovered PDE can be further considered as the latest assumption, based on which a completely physics-based recurrent network can be established. In the fine-tuning step, the functional form of the PDE is unchanged while its scalar coefficients are optimized. That said, the fine-tuning step is an optimization process on a much smaller parameter space. (ii) The sparse regression is performed on subsampled high-dimensional spatiotemporal data while the fine-tuning step could utilize all the available measurement through the recurrent network.
The above argument provides some intuitions of the coefficient fine-tuning step. To further justify it, we compute the relative ℓ2 error of the coefficient vector before/after the fine-tuning step, which is given in Table G.1. It can be seen that fine-tuning step improves the accuracy of the discovered PDE consistently for different noise level. Interestingly, we observe that for the 2D GS and λ–Ω RD case under 10% noise, the improvement brought by the fine tuning step is not as evident as that of the other cases. That is because, in these two cases, the sparse regression does not give the correct form of the PDE. In other words, the fine-tuning is performed based on an assumption that deviates from the fact. The dependency on the sparse regression result also implies one limitation of our approach.
Table G.1: Relative ℓ2 error (unit: 10−2) of the coefficient vector before and after the fine tuning.
Case Stage Noise level
0% 5% 10%
Burgers’ Before 0.93 1.57 2.18
After 0.50 0.54 0.59
λ–Ω RD Before 1.52 3.33 6.65
After 1.18 2.69 5.44
GS RD Before 6.40 9.84 12.33
After 1.59 2.85 10.03

heading: H PARETO ANALYSIS FOR HYPERPARAMETER SELECTION
text: As introduced in Section 3.3, sparse regression is formulated as an optimization problem described by Eq. (4). The aggregated optimization objective consists of two components, ||Ut −Θ(U)Ξ||2 and ||Ξ||0, to ensure the accuracy of the regression while keeping the model parsimonious. This trade-off can be balanced by selecting an appropriate weighting coefficient γ in the aggregated objective function. Pareto analysis is so far considered as the most effective way for selecting this hyperparameter γ. To ensure the scale of two components aligned, we introduce an equivalent hyperparameter κ as the nondimensional weighting coefficient. The κ and γ are interconvertible with the definition γ ≡ κ||Ut − Θ(U)Ξ̃||2 where Ξ̃ is the Least Squares solution to the linear system.
In the Pareto analysis, we would sample some points of κ from [10−2, 20]. For each of the point, we would pass it as the argument for the sparse regression routine, which would return the optimal coefficient vector Ξ∗. The corresponding regression error and ℓ0 penalty can be computed for Pareto analysis. Figure H.1 provides the results of Pareto analysis for 2D λ–Ω equation system (10% noise case). From the subfigure on the left, we can identify the Pareto front that corresponds to the optimal region of κ. Pareto analyses as such are performed on each case to find the suitable weighting coefficient κ (or equivalently γ). The value of κ employed in each case is summarized in Table H.1.
However, readers should note that the Pareto front may not always be so obvious to identify. In many scenarios, trade-off between the model sparsity and accuracy must be done with manual selection, which still appears to be more art than science.
Optimal region
(a)
(b)
Optimal region
𝜅𝜅
𝜅𝜅
*
*
*
*
*
* *
*
Figure H.1: Pareto analysis for selecting the weighting coefficient κ (or γ) in 2D λ–Ω GS system (10% noise level). (a) and (b) correspond to the Fu and Fv respectively. The Pareto front that represents the optimal κ is marked with the red star.
Table H.1: The coefficient κ employed in sparse regression.
Method Burgers’ λ–Ω RD GS RD
0% 5% 10% 0% 5% 10% 0% 5% 10%
Ours 1 1 1 1 1 1 0.1 0.1 1
PDE-FIND 1 1 1 1 1 1 0.1 0.05 0.1

heading: I DISCOVERY OF BURGERS’ EQUATION UNDER HIGHER REYNOLDS NUMBER
text: In this part, we examine the performance of our method on Burgers’ equation when Reynolds number grows. Higher Reynolds number would complicate the equation discovery problem due to the intense convective phenomena (or even chaos). We increase the Reynolds number of the problem to 500 by setting the ν = 0.002 as both the average initial velocity magnitude and characteristic length is 1. 10% Gaussian noise is adopted in this example while the remaining setting follows the Section 4.4. One snapshot of the reconstructed HR data, altogether with the ground truth, is presented in Fig. I.1. It can be observed that the reconstructed HR agree well with the ground truth. Then the reconstructed HR data is used in sparse regression to obtain the structure of the governing PDE, based on which the PeRCNN is built for coefficient fine tuning. The discovered PDEs after the coefficient fine tuning are
ut = 1.9854× 10−3∆u− 1.0072uux − 1.006vuy, vt = 1.9409× 10−3∆v − 0.9953uvx − 1.008vvy.
(12)
which match the ground truth correctly.
J INTERPRET THE LEARNED MODEL
One major drawback of the traditional deep neural network is the lack of interpretability as the network’s output is usually expressed as a prolonged nested function. However, since each channel of
the input to Π-block (i.e., Û (k)
) corresponds to a state variable component (i.e., [u, v]), the multiplicative form of Π-block makes it possible to extract (or interpret) the explicit form of learned F from the learned weights and biases via symbolic computations. This is a notable advantage of the proposed Π-block as an analytical expression would help researcher to disentangle the underlying physics or make inference on different cases (e.g., different I/BCs).
Here a simple experiment is conducted on the 2D Burgers’ case. The network employed in the experiment has two Conv layers with two channels. Two channels of the first Conv layer are associated with ∂/∂x and ∂/∂y respectively, by fixing the filters with the corresponding FD stencils. The remaining settings are kept the same as the 2D Burgers’ case in Section 4.4. The interpreted expression from the recurrent block (i.e., Π-block and the highway diffusion layer) model is
ut = [ 0.0051∆u− 0.95ux(1.07u− 0.0065v − 0.17) + 0.98uy(0.0045u− 1.01v + 0.17) + 0.053 0.0051∆v − 0.82vx(1.22u+ 0.0078v − 0.18)− 0.91vy(0.0063u+ 1.08v − 0.17) + 0.058 ] (13)
It is seen that the equivalent expression of the learned model is close to the genuine governing PDEs, which helps to explain the extraordinary expressiveness of our model. Although the selection of differential operators to be frozen is crucial for identifying the genuine form of the F , this example demonstrates the better interpretability of PeRCNN over “black-box” models. Similarly, with the trained model from Section 4.4 for 2D GS RD equation which has three Conv layers and 1×1 filters, we could extract a third degree polynomial for the reaction term in F .

Review summaries:
Problem: data-driven PDE discovery methods are not robus to low-quality measurement data.  Solution: A novel architecture that can encode prior knowledge (known terms, PDE structure, boundary conditions), and sparse regression procedure that is hypothesized to be more robust to scarce and noisy data scenarios.  Other Contributions: * proof that proposed Pi-block is a universal polynomial approximator. This is why I scored a "4" on technical contributions.

The authors propose a physics-encoded discrete learning framework to improve the robustness of the PDE-FIND algorithm. The physics knowledge such as known terms and initial conditions can be used to reconstruct higher-resolution data. Then, sparse regression is performed with fine-tuning of coefficients to finally determine the discovered PDE.

The focus of this paper is to develop a learning framework to discover spatiotemporal PDEs from scarce and noisy data.  The authors suggest devise a deep convolutional-recurrent network which encodes prior physics knowledge followed by  sparse regression with reconstructed data to identify the analytical form of the governing PDEs. The framework is validated  on three high-dimensional PDE systems where superiority over baselines is demonstrated.

The paper focuses on discovering the underlying physics dynamics, in the form of partial differential equation (PDE) from low-resolution (LR) measurements. It leverages sparse regression-based recovery to identify the constituents of the unknown PDE. The main contribution is a module that constructs high-resolution (HR) observations from LR data, for sparse regression to perform well.