Data-Driven Discovery of PDEs via the Adjoint Method

## Abstract

Abstract In this work, we present an adjoint-based method for discovering the underlying governing partial differential equations (PDEs) given data. The idea is to consider a parameterized PDE in a general form, and formulate the optimization problem that minimizes the error of PDE solution from data.
Using variational calculus, we obtain an evolution equation for the Lagrange multipliers (adjoint equations) allowing us to compute the gradient of the objective function with respect to the parameters of PDEs given data in a straightforward manner. In particular, for a family of parameterized and nonlinear PDEs, we show how the corresponding adjoint equations can be derived.
Here, we show that given smooth data set, the proposed adjoint method can recover the true PDE up to machine accuracy. However, in the presence of noise, the accuracy of the adjoint method becomes comparable to the famous PDE Functional Identification of Nonlinear Dynamics method known as PDE-FIND (Rudy et al., 2017 ) .
Even though the presented adjoint method relies on forward/backward solvers, it outperforms PDE-FIND for large data sets thanks to the analytic expressions for gradients of the cost function with respect to each PDE parameter.

## 1 Introduction

A large portion of data-driven modelling of physical processes in literature is dedicated to deploying Neural Networks to obtain fast prediction given the training data set. The data-driven estimation methods include Physics-Informed Neural Networks (Raissi et al., ), Pseudo-Hamiltonian neural networks (Eidnes & Lye, ), structure preserving (Matsubara et al., ; Sawant et al., ), and reduced order modelling (Duan & Hesthaven, ). These methods often provide efficient and somewhat ”accurate” predictions when tested as an interpolation method in the space of input or boundary parameters. Such fast estimators are beneficial when many predictions of a dynamic system is needed, for example in the shape optimization task in fluid dynamics.
However, the data-driven estimators often fail to provide accurate solution to the dynamical system when tested outside the training space, i.e. for extrapolation. Furthermore, given the regression-based nature of these predictors, often they do not offer any error estimator in prediction. Since we already have access to an arsenal of numerical methods in solving traditional governing equations, it is attractive to learn the underlying governing equation given data instead. Once the governing equation is found, one can use the standard and efficient numerical methods for prediction. This way we guarantee the consistency with observed data, estimator for the numerical approximation, and interoperability. Hence, learning the underlying physics given data has motivated a new branch in the scientific machine learning for discovering the mathematical expression as the governing equation given data.

The wide literature of data-driven discovery of dynamical systems includes equation-free modeling (Kevrekidis et al., ), artificial neural networks (González-García et al., ), nonlinear regression (Voss et al., ), empirical dynamic modeling (Sugihara et al., ; Ye et al., ), modeling emergent behavior (Roberts, ), automated inference of dynamics (Schmidt et al., ; Daniels & Nemenman, , ), normal form identification in climate (Majda et al., ), nonlinear Laplacian spectral analysis (Giannakis & Majda, ) and Koopman analysis (Mezić, ) among others. There has been a significant advancement in this field by combining symbolic regression wtih the evolutionary algorithms (Bongard & Lipson, ; Schmidt & Lipson, ; Tohme et al., ), which enable the direct extraction of nonlinear dynamical system information from data. Furthermore, the concept of sparsity (Tibshirani, ) has recently been employed to efficiently and robustly deduce the underlying principles of dynamical systems (Brunton et al., ; Mangan et al., ).

#### Related work.

Next, we review several relevant works that have shaped the current landscape of discovering PDEs from data:

*PDE-FIND* (Rudy et al., ). This method has been developed to discover underlying partial differential equation by minimizing the $L_{2}^{2}$-norm point-wise error of the parameterized forward model from the data. Estimating all the possible derivatives using finite difference, PDE-FIND constructs a dictionary of possible terms and finds the underlying PDE by performing a sparse search using ridge regression problem with hard thresholding, also known as STRidge optimization method. Several further developments in the literature has been carried out based on this idea, namely (Champion et al., ; Kaheman et al., ). In these methods, as the size (or dimension) of the data set increases, the PDE discovery optimization problem based on point-wise error becomes extremely expensive, forcing the user to arbitrarily reduce the size of data by resampling, or compress the data using proper orthogonal decomposition. Needles to say, in case of non-linear dynamics, such truncation of data can introduce bias in prediction leading to finding a wrong PDE.

*PDE-Net* (Long et al., , ). In this method, the PDE is learned from data using convolution kernels rather than brute-force use of finite differences, and apply neural networks to approximate nonlinear responses. Similar to PDE-FIND, the loss function of PDE-Net is the point-wise error from data which leads to a regression task that does not scale well with the size of data set.

*Hidden Physics Models* (Raissi & Karniadakis, ). This method assumes that the relevant terms of the governing PDE are already identified and finds its unknown parameters using Gaussian process regression (GPR). While GPR is an accurate interpolator which offers an estimate for the uncertainty in prediction, its training scales poorly with the size of the training data set as it requires inversion of the covariance matrix.

*PINN-SR* (Chen et al., ). One of the issues with the PDE-FIND is the use of finite difference in estimating the derivatives. The idea of PINN-SR is to extend PDE-FIND’s optimization problem to also find a PINN fit to the data in order to find smooth estimates for spatial derivatives. In particular, the training of PINN-SR combines the search for weights/biases of PINN approximation of the PDE with the sparse search in the space of possible terms to find the coefficients of the PDE given data. However, similar to PDE-FIND, the point-wise error from data is used as the loss for the regression task which does not scale well with the size of data set.

#### Contributions.

In this paper, we introduce a novel approach for discovering PDEs from data based on the adjoint method.
The idea is to formulate the objective (or cost) functional such that the estimate function $f$ minimizes the $L_{2}^{2}$-norm error from the data points $f^{*}$ with the constraint that $f$ is the solution to a parameterized PDE using the method of Lagrange multipliers. Here, we consider a parameterized PDE in a general form and the task is to find all the parameters including irrelevant ones. By finding the variational extremum of the cost functional with respect to the function $f$, we obtain a backward-in-time evolution equation for the Lagrange multipliers (adjoint equations). Next, we solve the forward parameterized PDE as well as the adjoint equations numerically. Having found estimates of the Lagrange multipliers and solution to the forward model $f$, we can numerically compute the gradient of the objective function with respect to the parameters of PDEs given data in a straightforward manner. In particular, for a family of parameterized and nonlinear PDEs, we show how the corresponding adjoint equations can be elegantly derived.
We note that the adjoint method has been successfully used before as an efficient shape optimization and sensitivity analysis method in fluid mechanics (Jameson, ; Caflisch et al., ) and plasma physics (Antonsen et al., ; Geraldini et al., ).

The remainder of the paper is organized as follows. First in Section [2](#S2), we introduce and derive the proposed adjoint-based method of finding the underlying system of PDEs given data. Next in Section [3](#S3), we present our results on a wide variety of PDEs and compare the solution with the celebrated PDE-FIND in terms of error and computational/training time. In Section [4](#S4), we discuss the limitations for the current version of our approach and provide concluding remarks in Section [5](#S5).

## 2 Adjoint method for finding PDEs

In this section, we introduce the problem and derive the proposed adjoint method for finding governing equations given data.

#### Problem setup.

Assume we are given a data set on a spatial/temporal grid $(\bm{x},t)$ for the vector of functions $\bm{f}^{*}(\bm{x},t)$. Here, $\bm{x}\in\Omega\subset\mathbb{R}^{n}$ is coordinates inside the solution domain $\Omega$, $t\in[0,T]$ denotes time, and output is a map $\bm{f}^{*}:(\Omega,[0,T])\rightarrow\mathbb{R}^{N}$. The goal is to find the governing equations that accurately estimates $\bm{f}^{*}$ at all points. In order to achieve this goal, we formulate the problem using the method of Lagrange multipliers.

#### Adjoint method.

We consider a forward model $\mathcal{L}[\cdot]$ that maps an $N$-dimensional vector of functions $\bm{f}(\bm{x},0)$ into $\bm{f}(\bm{x},t)$ for $t\in[0,T]$ where the $i$-th PDE is given by

$$ $\displaystyle\mathcal{L}_{i}[\bm{f}]:=\partial_{t}f_{i}+\sum_{\bm{d},\bm{p}}\alpha_{i,\bm{d},\bm{p}}\nabla^{(\bm{d})}_{{\bm{x}}}[\bm{f}^{\bm{p}}]=0$ (1) $$

for $i=1,\ldots,N$, resulting in a system of N-PDEs, i.e. the $i$-th PDE $\mathcal{L}_{i}$ predicts $f_{i}$. Here, $\bm{x}=[x_{1},x_{2},\ldots,x_{n}]$ is an n-dimensional (spatial) input vector, and $\bm{f}=[f_{1},f_{2},\ldots,f_{N}]$ is an N-dimensional vector of functions. Here, we use the shorthand $f_{i}=f_{i}(\bm{x},t)$ and $\bm{f}=\bm{f}(\bm{x},t)$. Furthermore, $\bm{p}=[p_{1},\ldots,p_{N}]$ and $\bm{d}=[d_{1},d_{2},\ldots,d_{n}]$ are non-negative index vectors such that $\bm{f}^{\bm{p}}=f_{1}^{p_{1}}f_{2}^{p_{2}}\cdots f_{N}^{p_{N}}$ and

$$ $\displaystyle\nabla^{(\bm{d})}_{{\bm{x}}}$ $\displaystyle=\nabla^{(d_{1})}_{x_{1}}\nabla^{(d_{2})}_{x_{2}}\cdots\nabla^{(d_{n})}_{x_{n}}\leavevmode\nobreak\ ,$ (2) $$

where $\nabla^{(d_{i})}_{x_{i}}$ for $i=1,...,n$ indicates $d_{i}$-th derivative in $x_{i}$ dimension, and $\partial_{t}f_{i}$ denotes the time derivative of the $i$-th function. We denote the vector of unknown parameters by $\bm{\alpha}=[\alpha_{i,\bm{d},\bm{p}}]_{(i,\bm{d},\bm{p})\in\mathcal{D}}$, where $\mathcal{D}$ represents the domain of all valid combinations of $i$, $\bm{d}$, and $\bm{p}$.

Having written the forward model ([1](#S2.E1)) as general as possible, the goal is to find the parameters ${\bm{\alpha}}$ such that $\bm{f}$ approximates the data points of $\bm{f}^{*}$. To this end, we formulate an objective (or cost) functional that minimizes the $L_{2}^{2}$-norm error between what the model predicts and the data $\bm{f}^{*}$ for all $\bm{x}$ and $t$, with the constraint that $\bm{f}$ solves the forward model ([1](#S2.E1)), i.e.

$$ $\displaystyle\mathcal{C}=\sum_{i=1}^{N}\bigg{(}$ $\displaystyle\int\big{(}f_{i}^{*}(\bm{x},t)-f_{i}(\bm{x},t)\big{)}^{2}d\bm{x}dt$ $\displaystyle+$ $\displaystyle\int\big{(}f_{i}^{*}(\bm{x},T)-f_{i}(\bm{x},T)\big{)}^{2}d\bm{x}$ $\displaystyle+$ $\displaystyle\int\lambda_{i}(\bm{x},t)\mathcal{L}_{i}[\bm{f}(\bm{x},t)]d\bm{x}dt\bigg{)}\leavevmode\nobreak\ .$ (3) $$

Clearly, given estimates of $\bm{f}$ and Lagrange multipliers $\bm{\lambda}=(\lambda_{1},\lambda_{2},\ldots,\lambda_{N})$, the gradient of the cost function with respect to model parameters can be simply computed via

$$ $\displaystyle\frac{\partial\mathcal{C}}{\partial\alpha_{i,\bm{d},\bm{p}}}=(-1)^{|\bm{d}|}\int\bm{f}^{\bm{p}}\,\nabla^{(\bm{d})}_{{\bm{x}}}[\lambda_{i}]d\bm{x}dt$ (4) $$

where $i=1,...,N$ and $|\bm{d}|=d_{1}+...+d_{n}$, where $|.|$ denotes $L_{1}$-norm.
Here, we used integration by parts and imposed the condition that $\bm{\lambda}\rightarrow\bm{0}$ on the boundaries of $\Omega$ at all time $t\in[0,T]$. The analytical expression ([4](#S2.E4)) can be used for finding the parameters of PDE using in the gradient descent method with update rule

$$ $\displaystyle\alpha_{i,\bm{d},\bm{p}}\leftarrow\alpha_{i,\bm{d},\bm{p}}-\eta\frac{\partial\mathcal{C}}{\partial\alpha_{i,\bm{d},\bm{p}}}$ (5) $$

for $i=1,\ldots,N$, where $\eta=\beta\min(\Delta\bm{x})^{|\bm{d}|-d_{\mathrm{max}}}$ is the learning rate which includes a free parameter $\beta$ and scaling coefficient for each term of the PDE, and $d_{\mathrm{max}}=\max(|\bm{d}|)$ for all considered $\bm{d}$. Let us also define $p_{\mathrm{max}}=\max(|\bm{p}|)$ as the highest order in the forward PDE model. We note that since the terms of the PDEs may have different scaling, the step size for the corresponding coefficient must be adjusted accordingly.

However, before we can use Eq. ([4](#S2.E4)) and ([5](#S2.E5)), we need to find $\bm{\lambda}$, hence the adjoint equation. This can be achieved by finding the functional extremum of the cost functional $\mathcal{C}$ with respect to $\bm{f}$. First, we note that the total variation of $\mathcal{C}$ can be derived as

$$ $\displaystyle\delta\mathcal{C}=$ $\displaystyle\sum_{i=1}^{N}\bigg{(}-\int 2\big{(}f_{i}^{*}{(\bm{x},t)}-f_{i}{(\bm{x},t)}\big{)}\delta f_{i}d\bm{x}dt$ $\displaystyle-$ $\displaystyle\int 2\big{(}f_{i}^{*}(\bm{x},T)-f_{i}(\bm{x},T)\big{)}\delta f_{i,T}d\bm{x}$ $\displaystyle+$ $\displaystyle\int\bigg{(}-\frac{\partial\lambda_{i}}{\partial t}+\sum_{\bm{d},\bm{p}}(-1)^{|\bm{d}|}\alpha_{i,\bm{d},\bm{p}}\nabla_{f_{i}}[\bm{f}^{\bm{p}}]\nabla^{(\bm{d})}_{{\bm{x}}}[\lambda_{i}]\bigg{)}\delta f_{i}d\bm{x}dt$ $\displaystyle+$ $\displaystyle\int\lambda_{i}(\bm{x},T)\delta f_{i,T}d\bm{x}\bigg{)}$ (6) $$

where $\delta f_{i}$ denotes variation with respect to $f_{i}(\bm{x},t)$, and $\delta f_{i,T}$ variation with respect to $f_{i}(\bm{x},T)$. Here again, we used integration by parts and imposed the condition that $\bm{\lambda}\rightarrow\bm{0}$ on the boundaries of $\Omega$ at all time $t\in[0,T]$ for $i=1,...,N$. Note that $f_{i}(\bm{x},t)$ is the output of $i$-th PDE.

Next, we find the optimums of $\mathcal{C}$ (and hence the adjoint equations) by taking the variational derivatives with respect to $f_{i}$ and $f_{i,T}$, i.e.

$$ $\displaystyle\frac{\delta\mathcal{C}}{\delta f_{i}}$ $\displaystyle=0\implies$ $\displaystyle\frac{\partial\lambda_{i}}{\partial t}$ $\displaystyle=\sum_{\bm{d},\bm{p}}(-1)^{|\bm{d}|}\alpha_{i,\bm{d},\bm{p}}\nabla_{f_{i}}[\bm{f}^{\bm{p}}]\nabla^{(\bm{d})}_{{\bm{x}}}[\lambda_{i}]-2(f^{*}_{i}-f_{i})$ (7) $$

and

$$ $\displaystyle\frac{\delta\mathcal{C}}{\delta f_{i,T}}$ $\displaystyle=0\implies\lambda_{i}(\bm{x},T)=2\big{(}f^{*}_{i}(\bm{x},T)-f_{i}(\bm{x},T)\big{)}$ (8) $$

for $i=1,...,N$. We note that the adjoint equation ([7](#S2.E7)) for the system of PDEs is backward in time with the condition at the final time $t=T$ given by Eq. ([8](#S2.E8)). In order to make the notation clear, we present an examples for deriving the adjoint equations in Appendix [A](#A1).

The training procedure follows the standard gradient descent method, starting by taking an initial guess for parameters $\bm{\alpha}$. Then, the forward model ([1](#S2.E1)) is solved numerically to estimate $\bm{f}(\bm{x},t)$ given the initial condition $\bm{f}^{*}(\bm{x},0)$. Then, the adjoint Eq. ([7](#S2.E7)) is solved backwards in time with the final time condition  ([8](#S2.E8)). Finally, the estimate for parameters of the model is updated using Eq. ([5](#S2.E5)). In Algorithm [1](#alg1), we present a pseudocode for finding the parameters of the system of PDEs using Adjoint method. For simplicity, in this work we use finite difference method to estimate the spatial and Euler for the time derivative, for both forward and backward equations with similar mesh size as the input data $\bm{f}^{*}$ on $\Omega\times[0,T]$.

Figure: Algorithm 1 Finding a system of PDEs using Adjoint method

## 3 Results

We demonstrate the validity of our proposed adjoint-based method in discovering PDEs given time series measurements on a spatial-temporal grid. We consider data collected from a variety of systems. We mainly compare our approach to PDE-FIND in terms of error and time to convergence. All the results are obtained using a single core-thread of a 2.3 GHz Quad-Core Intel Core i7 CPU.

### 3.1 Heat equation

As a first example, let us consider measured data collected from the solution to the heat equation, i.e.

$$ $\displaystyle\frac{\partial f}{\partial t}+D\,\frac{\partial^{2}f}{\partial x^{2}}=0,$ (9) $$

with $D=-1$. The data is constructed using finite difference with initial condition $f(x,0)=5\sin(2\pi x)x(x-L)$ and a mesh with $N_{x}=100$ nodes in $x$ covering the domain $\Omega=[0,L]$ with $L=1$ and $N_{t}=100$ steps in $t$ with final time $T=N_{t}\Delta t$ where $\Delta t=0.05\Delta x^{2}/(1+|D|)$ is the step size and $\Delta x=L/N_{x}$ is the mesh size in $x$.
We consider a system consisting of a single PDE (i.e. $N=1$, $\bm{f}=f$, and $\bm{p}=p$) with one-dimensional input, i.e. $n=1$ and $\bm{x}=x\subset\mathbb{R}$, and $\bm{d}=d\in\mathbb{N}$. In order to construct a general forward model, here we consider derivatives and polynomials with indices $d,p\in\{1,2,3\}$ as the initial guess for the forward model. This leads to $9$ terms with unknown coefficients $\bm{\alpha}$ that we find using the proposed adjoint method (an illustrative derivation of the candidate terms can be found in Appendix [A.1](#A1.SS1)). While we expect to recover the coefficient that corresponds to $D$, we expect all the other coefficients (denoted by $\bm{\alpha}^{*}$) to become negligible. That is what we indeed observe in Fig. [1](#S3.F1) where the error of the coefficient for each term is plotted against the number of epochs.

Figure: Figure 1: The estimated coefficient corresponding to $D$ (left) and the $L_{1}$-norm error of all considered coefficients (right).
Refer to caption: /html/2401.17177/assets/x1.png

Next, we compare the solution obtained via the adjoint method against PDE-FIND with STRidge optimization method. Here, we test both methods in recovering the Heat equation given data on the grid with descritization $(N_{t},N_{x})\in\{(100,100),\ (500,100),\ (1000,100),\ (1000,1000)\}$. As shown in Fig. [2](#S3.F2), the proposed adjoint method provides more accurate results across all data sizes. We also point out that as the size of the data set increases, PDE-FIND with STRidge regression method becomes more expensive, e.g. one order of magnitude more expensive than the adjoint method for the data on a grid size $(N_{t},N_{x})=(1000,1000)$.

Figure: Figure 2: $L_{1}$-norm error of the estimated coefficients (left) and the execution time (right) for discovering the heat equation equation using the proposed Adjoint method (blue) and PDE-FIND method (red), given data on a grid with $N_{t}\in\{100,500,1000\}$ steps in $t$, and $N_{x}\in\{100,1000\}$ nodes in $x$.
Refer to caption: /html/2401.17177/assets/x3.png

### 3.2 Burgers’ equation

As a nonlinear test case, let us consider the data from Burgers’ equation given by

$$ $\displaystyle\frac{\partial f}{\partial t}+\frac{\partial(Af^{2})}{\partial x}=0$ (10) $$

where $A=-1$. The data is obtained with similar simulation setup as for Heat equation (Section [3.1](#S3.SS1)) except for the time step, i.e. $\Delta t=0.05\Delta x/(1+|A|)$.

Similar to Section [3.1](#S3.SS1), we adopt a system of one PDE with one-dimensional input. We also consider derivatives and polynomials with indices $d,p\in\{1,2,3\}$ in the construction of the forward model. This leads to $9$ terms whose coefficients we find using the proposed adjoint method. As shown in Fig. [3](#S3.F3), the proposed adjoint method finds the correct coefficients, i.e. $\alpha_{d=1,p=2}$ that corresponds to $D$ as well as all the irrelevant ones denoted by $\bm{\alpha}^{*}$, up to machine accuracy in $\mathcal{O}(10)$ epochs.

Figure: Figure 3: The estimated coefficient corresponding to $A$ (left) and the $L_{1}$-norm error of all considered coefficients (right) given the discretized data of Burgers’ equation during training.
Refer to caption: /html/2401.17177/assets/x5.png

Next, we compare the solution obtained from the adjoint method to the one from PDE-FIND using STRidge optimization method. Here, we compare the error on coefficients and computational time between the adjoint and PDE-FIND by repeating the task for the data set with increasing size, i.e. $(N_{t},N_{x})\in\{(100,100),\ (1000,100),\ (1000,\ 1000)\}$. As depicted in Fig. [4](#S3.F4), the adjoint method provides us with more accurate solution across the different discretization sizes. Regarding the computational cost, while PDE-FIND seems faster on smaller data sets, as the size of the data grows, it becomes increasingly more expensive than the adjoint method. Similar to the heat equation, for the mesh size $(N_{t},N_{x})=(1000,\ 1000)$ we obtain one order of magnitude speed up compared to PDE-FIND.

Figure: Figure 4: $L_{1}$-norm error of the estimated coefficients (left) and the execution time (right) for discovering the Burgers’ equation equation using the Adjoint method (blue) and PDE-FIND method (red), given data on a grid with $N_{t}\in\{100,1000\}$ steps in $t$, and $N_{x}\in\{100,1000\}$ nodes in $x$.
Refer to caption: /html/2401.17177/assets/x7.png

### 3.3 Kuramoto Sivashinsky equation

As a more challenging test case, let us consider the recovery of the Kuramoto-Sivashinsky (KS) equation given by

$$ $\displaystyle\frac{\partial f}{\partial t}+A\frac{\partial f^{2}}{\partial x}+B\frac{\partial^{2}f}{\partial x^{2}}+C\frac{\partial^{4}f}{\partial x^{4}}=0$ (11) $$

where $A=-1$, $B=0.5$ and
$C=-0.5$. The data is generated similar to previous sections except for the grid $(N_{t},N_{x})=(64,256)$ and the time step size $\Delta t=0.01\Delta x^{4}/(1+|C|)$.

Here again, we adopt a system of one PDE with one-dimensional input. As a guess for the forward model, we consider terms consisting of derivatives with indices $d\in\{1,2,3,4\}$ and polynomials with indices $b\in\{1,2\}$, leading to $8$ terms whose coefficients we find using the proposed adjoint method. As shown in Fig. [5](#S3.F5), the adjoint method finds the coefficient with error of $\mathcal{O}(10^{-5})$, yet achieving machine accuracy seems not possible.

Figure: Figure 5: The estimated coefficients corresponding to $A,B,C$ (left) and the $L_{1}$-norm error of all considered coefficients (right) given the discretized data of the KS equation during training.
Refer to caption: /html/2401.17177/assets/x9.png

Again, in Fig. [6](#S3.F6) we make a comparison between the predicted PDE using the adjoint method against PDE-FIND. In particular, we consider a data set on a temporal/spatial mesh of size $(N_{t},N_{x})=\{(64,256),\ (128,512),\ (256,1024)\}$ and compare how the error and computational cost vary. Similar to previous sections, the error is reported by comparing the obtained coefficients against the coefficients of the exact PDE in $L_{1}$-norm. Interestingly, the PDE-FIND method has $3$ to $4$ orders of magnitude larger error compared to the adjoint method. Also, in terms of cost, the training time for PDE-FIND seems to grow at a higher rate than the adjoint method as the (data) mesh size increases.

Figure: Figure 6: $L_{1}$-norm error of the estimated coefficients (left) and the execution time (right) for discovering the KS equation using the Adjoint method (blue) and PDE-FIND method (red), given data on a grid with $N_{t}\in\{64,128,256\}$ steps in $t$, $N_{x}\in\{256,512,1024\}$ nodes in $x$.
Refer to caption: /html/2401.17177/assets/x11.png

### 3.4 Random Walk

Next, let us consider the recovery of the governing equation on probability density function (PDF) given samples of its underlying stochastic process. As an example, we consider the Itô process

$$ $\displaystyle dX=Adt+\sqrt{2D}dW$ (12) $$

where $A=1$ is drift and $D=0.5$ is the diffusion coefficient, and $W$ denotes the standard Wiener process with $\mathrm{Var}(dW)=\Delta t$. We generate the data set by simulating the random walk using Euler-Maruyama scheme starting from $X(t=0)=0$ for $N_{t}=50$ steps with a time step size of $\Delta t=0.01$. We estimate the PDF using histogram with $N_{x}=100$ bins and $N_{\mathrm{s}}=1000$ samples.
Let us denote the distribution of $X$ by $f$. Itô’s lemma gives us the Fokker-Planck equation

$$ $\displaystyle\frac{\partial f}{\partial t}+A\frac{\partial f}{\partial x}-D\frac{\partial^{2}f}{\partial x^{2}}=0\leavevmode\nobreak\ .$ (13) $$

Given the data set for $f$ on a mesh of size $(N_{t},N_{x})$, we can use finite difference to compute the contributions from derivatives of $f$ in the governing law. Since this is one of the challenging test cases due to noise, here we only consider three three possible terms in the forward model, consisting of derivatives with indices $d\in\{1,2,3\}$ and polynomial power $p=1$. In Fig. [7](#S3.F7), we show how the error of finding the correct coefficients evolves during training for the adjoint method. Clearly, the adjoint method seems to recover the true PDE with $L_{1}$ error of $\mathcal{O}(10^{-2})$ in its coefficients.

Figure: Figure 7: The estimated coefficients corresponding to $A$ and $D$ (left) and the $L_{1}$-norm error of all considered coefficients (right) of the Fokker-Planck equation as the governing law for the PDF associated with the random walk during training.
Refer to caption: /html/2401.17177/assets/x13.png

In Fig. [8](#S3.F8), we make a comparison with PDE-FIND for the same number and order of terms as the initial guess for the PDE. We compare the two methods for a range of grid and sample sizes, i.e. $N_{t}\in\{50,100\}$, $N_{x}=100$, and $N_{s}\in\{10^{3},\ 10^{4}\}$. It turns out that the proposed adjoint method overall provides more accurate estimate of the coefficients than PDE-FIND, though at a higher cost. In Table [1](#S3.T1), we show the discovered PDEs for both methods across the different discretizations.

Figure: Figure 8: $L_{1}$-norm error of the estimated coefficients (left) and the execution time (right) for discovering the Fokker-Planck equation using the proposed Adjoint method (blue) and PDE-FIND method (red), given samples of its underlying stochastic process with $N_{t}\in\{50,100\}$ steps in $t$, $N_{x}=100$ histogram bins, and $N_{s}\in\{10^{3},10^{4}\}$ samples.
Refer to caption: /html/2401.17177/assets/x15.png

**Table 1: Recovery of the Fokker-Planck equation, i.e. $f_{t}+f_{x}-0.5f_{xx}=0$, using the proposed adjoint method against PDE-FIND method given samples of the underlying stochastic process for various discretization parameters.**
| $\bm{N_{t}}$ | $\bm{N_{x}}$ | $\bm{N_{\mathrm{s}}}$ | Method | Recovered PDE |
| --- | --- | --- | --- | --- |
| 50 | 100 | 1000 | Adjoint | $f_{t}+1.025f_{x}-0.465f_{xx}=0$ |
|  |  |  | PDE-FIND | $f_{t}+0.798f_{x}-0.454f_{xx}=0$ |
|  |  | 10000 | Adjoint | $f_{t}+1.022f_{x}-0.495f_{xx}=0$ |
|  |  |  | PDE-FIND | $f_{t}+0.818f_{x}-0.496f_{xx}=0$ |
| 100 | 100 | 1000 | Adjoint | $f_{t}+1.010f_{x}-0.543f_{xx}=0$ |
|  |  |  | PDE-FIND | $f_{t}+0.863f_{x}-0.560f_{xx}=0$ |
|  |  | 10000 | Adjoint | $f_{t}+1.015f_{x}-0.589f_{xx}=0$ |
|  |  |  | PDE-FIND | $f_{t}+0.894f_{x}-0.612f_{xx}=0$ |

### 3.5 Reaction Diffusion System of Equations

In order to show scalability and accuracy of the adjoint method for a system of PDEs in a higher dimensional space, let us consider a system of PDEs given by

$$ $\displaystyle\frac{\partial u}{\partial t}+c_{0}^{u}\nabla^{2}_{{x_{1}}}[u]+c_{1}^{u}\nabla^{2}_{{x_{2}}}[u]+R^{u}(u,v)=0,$ (14) $\displaystyle\frac{\partial v}{\partial t}+c_{0}^{v}\nabla^{2}_{{x_{1}}}[v]+c_{1}^{v}\nabla^{2}_{{x_{2}}}[v]+R^{v}(u,v)=0$ (15) $$

where

$$ $\displaystyle R^{u}(u,v)$ $\displaystyle=c^{u}_{2}u+c^{u}_{3}u^{3}+c^{u}_{4}uv^{2}+c^{u}_{5}u^{2}v+c^{u}_{6}v^{3}$ (16) $\displaystyle R^{v}(u,v)$ $\displaystyle=c^{v}_{2}v+c^{v}_{3}v^{3}+c^{v}_{4}vu^{2}+c^{v}_{5}v^{2}u+c^{v}_{6}u^{3}$ (17) $$

We construct the data set by solving the system of PDEs Eqs. ([14](#S3.E14))-([15](#S3.E15)) using a 2nd order finite difference scheme with initial values

$$ $\displaystyle u_{0}=a\sin\left(\frac{4\pi x_{1}}{L_{1}}\right)\cos\left(\frac{3\pi x_{2}}{L_{2}}\right)\left(L_{1}x_{1}-x_{1}^{2}\right)\left(L_{2}x_{2}-x_{2}^{2}\right)$ $\displaystyle v_{0}=a\cos\left(\frac{4\pi x_{1}}{L_{1}}\right)\sin\left(\frac{3\pi x_{2}}{L_{2}}\right)\left(L_{1}x_{1}-x_{1}^{2}\right)\left(L_{2}x_{2}-x_{2}^{2}\right)$ $$

where $a=100$, and the coefficients

$$ $\displaystyle\bm{c}^{u}=[c^{u}_{i}]_{i=0}^{6}=[-0.1,-0.2,-0.3,-0.4,0.1,0.2,0.3]$ $\displaystyle\bm{c}^{v}=[c^{v}_{i}]_{i=0}^{6}=[-0.4,-0.3,-0.2,-0.1,0.3,0.2,0.1].$ $$

We generate data by solving the system of PDEs Eqs. ([14](#S3.E14))-([15](#S3.E15)) using the finite difference method and forward Euler scheme for $N_{t}=25$ steps with a time step size of $\Delta t=10^{-6}$, and in the domain $\Omega=[0,L_{1}]\times[0,L_{2}]$ where $L_{1}=L_{2}=1$ which is discretized using a uniform grid with $N_{x_{1}}\times N_{x_{2}}=50^{2}$ nodes leading to mesh size $\Delta x_{1}=\Delta x_{2}=0.02$. In Fig. [9](#S3.F9) we show the solution to the system at time $T=N_{t}\Delta t$ for $u$ and $v$.

Figure: Figure 9: Solution to the reaction diffusion system of PDEs at time $t=T$ for $u$ (left) and $v$ (right).
Refer to caption: /html/2401.17177/assets/x17.png

We consider a system consisting of two PDEs, i.e. $N=\dim(\bm{f})=\dim(\bm{p})=2$, with two-dimensional input, i.e. $n=\dim(\bm{x})=\dim(\bm{d})=2$. Here, $\dim(\bm{f})=\dim(\text{Ima}(\bm{f}))$, where Ima$(\cdot)$ denotes the image (or output) of a function.

In order to use the developed adjoint method, we construct a guess forward system of PDEs (or forward model) using derivatives up to 2nd order and polynomials of up to 3rd order. That is, $d_{\mathrm{max}}=2$ and $p_{\mathrm{max}}=3$. This leads to $90$ terms whose coefficients we find using the proposed adjoint method (an illustrative derivation of the candidate terms can be found in Appendix [A.2](#A1.SS2)). The solution to the constructed model $\bm{f}\approx[u,v]$ as well as the adjoint equation for $\bm{\lambda}$ is found using the same discretization as the data set.

As shown in Fig. [10](#S3.F10), the adjoint method finds the correct equations with error up to $\mathcal{O}(10^{-12})$. Furthermore, the coefficients corresponding to the irrelevant terms $\bm{\alpha}^{*}$ tend to zero with error of $\mathcal{O}(10^{-11})$, see Fig. [11](#S3.F11).

Figure: Figure 10: $L_{1}$-norm error in the estimated coefficients of the reaction diffusion system of PDEs during training.
Refer to caption: /html/2401.17177/assets/x18.png

Figure: Figure 11: $L_{1}$-norm error in the estimated coefficients of the irrelevant terms compared to the true reaction diffusion system of PDEs during training, i.e. $||e(\bm{\alpha}^{*})||_{1}=||\bm{\alpha}^{*}||_{1}$.
Refer to caption: /html/2401.17177/assets/x20.png

Furthermore, we have compared the adjoint method against PDE-FIND for a range of grid sizes in Fig. [12](#S3.F12). We observe that the cost of PDE-FIND grows with higher rate than adjoint method as the size of the data set increases.

Figure: Figure 12: Comparing the error and execution time of the adjoint method (blue) to PDE-FIND method (red) against the size of the data set for the tolerance of $10^{-7}$ in the discovered coefficients.
Refer to caption: /html/2401.17177/assets/x21.png

### 3.6 Sensitivity to noise

So far, we have only tested the adjoint method given noise-free data. Here, we investigate how the error increases once noise is added to the data set. In particular, we add noise $\epsilon\sim\mathcal{N}(0,\sigma^{2})$ to each point of the data set for $\bm{f}^{*}(\bm{x},t)$, where $\mathcal{N}(0,\sigma^{2})$ denotes a normal distribution with zero mean and standard deviation $\sigma$. As test cases, we revisit the heat and burgers’ equations with added noise of $\epsilon$ with $\sigma\in\{0.001,\ 0.005,\ 0.01,\ 0.1\}\ \%$.

As shown in Figure [13](#S3.F13), adding noise to the data set deteriorates the accuracy in finding the correct coefficients of the underlying PDE for both the adjoint method and PDE-FIND method. While we observe that the accuracy of the adjoint method is less prone to noise compared to PDE-FIND, we conclude that both methods fail in recovering PDE at high levels of noise. While this motivates us to first consider filtering the noise in the data set, and then use the smoothed out data for the PDE recovery task, we leave out de-noising task here as it is beyond the scope of this paper.

Figure: Figure 13: Error of the adjoint method (blue) and PDE-FIND method (red) in finding the coefficients of the true PDE, i.e. heat equation (left) and Burgers’ equation (right), given noisy data.
Refer to caption: /html/2401.17177/assets/x23.png

## 4 Discussion

Below we highlight and discuss strengths and weaknesses of the proposed adjoint method.

#### Strengths.

The proposed method has several strengths:
i) The proposed adjoint-based method of discovering PDEs can provides coefficients of the true governing equation with significant accuracy.
ii) Since the gradient of the cost function with respect to parameters are derived analytically, the optimization problem converges fast. In particular, the adjoint method becomes cheaper than PDE-FIND as the size of the data set increases.
iii) Since the adjoint method uses a PDE solver to find the underlying governing equation, there’s a guarantee that the found PDE can be solved numerically with the same PDE solver as the one used by the adjoint method.

#### Weaknesses.

Our proposed method has some limitations:
i) In order to use the proposed adjoint method for discovering PDEs, a general solver of PDEs needs to be implemented. Here, we used finite differences which can be replaced with more advance solvers.
ii) Similar to PDE-FIND, the accuracy of the adjoint method reduces enormously when the data is noisy. This could be addressed through the use of noise reduction techniques, such as filtering out high frequencies using Fast Fourier Transform (Li et al., ) or smooth interpolating tools such as PINNs (Raissi et al., ), to increase the regularity of the data before carrying out the PDE-discovery task.
iii) In this work, we used the same discretization as the input data. If the grid of input data is too coarse for the PDE solver, one has to use interpolation to estimate the data on a finer grid that is more appropriate for the PDE solver.

## 5 Conclusion

In this work, we introduce a novel mathematical method for the discovery of partial differential equations given data using the adjoint method. By formulating the optimization problem in the variational form using the method of Lagrange multipliers, we find an analytic expression for the gradient of the cost function with respect to the parameters of the PDE as a function of the Lagrange multipliers and the forward model estimate. Then, using variational calculus, we find a backwards-in-time evolution equation for the Lagrange multipliers which incorporates the error with a source term (the adjoint equation). Hence, we can use the same solver for both forward model and the backward Lagrange equations. Here we used finite differences to estimate the spatial derivatives and forward Euler for the time derivatives, which indeed can be replaced with more stable and advanced solvers.

We compared the proposed adjoint method against PDE-FIND on several test cases. While PDE-FIND seems to be faster for small size problems, we observe that the adjoint method equipped with forward/backward solvers becomes faster than PDE-FIND as the size of the data set increases. Also, the adjoint method can provide machine-accuracy in identifying and finding the coefficients when the data set is noise-free. Furthermore, in the case of discovering PDEs for PDFs given its samples, both methods seem to suffer enormously form noise/bias associated with the finite number of samples and the finite difference on histogram. This motivates the use of smooth and least biased density estimator in these methods such as (Tohme et al., ) in future work. In the future work, we intend to combine the adjoint-based method for the discovery of PDE with PINNs as the solver instead of finite difference method. This would allow us to handle noisy and sparse data as well as deploying larger time steps in estimating the forward and backward solvers.

## Acknowledgements

T. Tohme was supported by the MathWorks Engineering Fellowship. M. Sadr thanks Dr. Andreas Adelmann for his support and Prof. Thomas Antonsen for his constructive discussion on the adjoint method.

## Appendix A Illustrative derivation of the adjoint equations for the considered cases.

Although the proposed method and its algorithm can be and has been computed in an automated fashion,
here we show two detailed illustrative examples for 1-dimensional and 2-dimensional cases presented in Section [3](#S3) for the sake of better understanding the used notation and how the library of candidate terms looks like.

### A.1 Heat and Burger’s Equations

As mentioned in Sections [3.1](#S3.SS1) and [3.2](#S3.SS2), for these two cases, we consider a system consisting of a single PDE, i.e. $N=\dim(\bm{f})=\dim(\bm{p})=1$ where $\bm{f}=f$ and $\bm{p}=p$, in one-dimensional input space, i.e. $n=\dim(\bm{x})=\dim(\bm{d})=1$ where $\bm{x}=x$, and $\bm{d}=d$.
In addition, we consider candidate terms consisting of derivatives with indices $d\in\{1,2,3\}$ and polynomials with indices $p\in\{1,2,3\}$). In other words, $d_{\mathrm{max}}=3$ and $p_{\mathrm{max}}=3$. The resulting forward model in Eq. [1](#S2.E1) takes the form

$$ $\displaystyle\mathcal{L}[f]=\frac{\partial f}{\partial t}$ $\displaystyle+\sum_{d=1}^{3}\sum_{p=1}^{3}\alpha_{d,p}\frac{\partial^{d}\big{(}f^{p}\big{)}}{\partial x^{d}}$ $\displaystyle=\frac{\partial f}{\partial t}$ $\displaystyle+\alpha_{1,1}\frac{\partial f}{\partial x}+\alpha_{1,2}\frac{\partial\big{(}f^{2}\big{)}}{\partial x}+\alpha_{1,3}\frac{\partial\big{(}f^{3}\big{)}}{\partial x}$ $\displaystyle+\alpha_{2,1}\frac{\partial^{2}f}{\partial x^{2}}+\alpha_{2,2}\frac{\partial^{2}\big{(}f^{2}\big{)}}{\partial x^{2}}+\alpha_{2,3}\frac{\partial^{2}\big{(}f^{3}\big{)}}{\partial x^{2}}$ $\displaystyle+\alpha_{3,1}\frac{\partial^{3}f}{\partial x^{3}}+\alpha_{3,2}\frac{\partial^{3}\big{(}f^{2}\big{)}}{\partial x^{3}}+\alpha_{3,3}\frac{\partial^{3}\big{(}f^{3}\big{)}}{\partial x^{3}}$ (18) $$

where $\alpha_{d,p}$ denotes the parameter corresponding to the term with $d$-th derivative and $p$-th polynomial order. As we can observe, we have $9$ terms with unknown coefficients $\bm{\alpha}=[\alpha_{d,p}]_{d\in\{1,2,3\},p\in\{1,2,3\}}$ that we aim to find using the proposed adjoint method.

The cost functional in this case is simply

$$ $\displaystyle\mathcal{C}$ $\displaystyle=\int\big{(}f^{*}(x,t)-f(x,t)\big{)}^{2}dxdt+\int\big{(}f^{*}(x,T)-f(x,T)\big{)}^{2}dx+\int\lambda(x,t)\mathcal{L}[f(x,t)]dxdt\leavevmode\nobreak\ .$ (19) $$

Letting variational derivatives of $\mathcal{C}$ with respect to $f$ to be zero, and using integration by parts, the corresponding adjoint equation can be obtained as

$$ $\displaystyle\frac{\partial\lambda}{\partial t}=-2(f^{*}-f)$ $\displaystyle+\sum_{d=1}^{3}\sum_{p=1}^{3}(-1)^{d}\alpha_{d,p}\frac{\partial\big{(}f^{p}\big{)}}{\partial f}\frac{\partial^{d}\lambda}{\partial x^{d}}$ $\displaystyle=-2(f^{*}-f)$ $\displaystyle-\alpha_{1,1}\frac{\partial\lambda}{\partial x}-\alpha_{1,2}\big{(}2f\big{)}\frac{\partial\lambda}{\partial x}-\alpha_{1,3}\big{(}3f^{2}\big{)}\frac{\partial\lambda}{\partial x}$ $\displaystyle+\alpha_{2,1}\frac{\partial^{2}\lambda}{\partial x^{2}}+\alpha_{2,2}\big{(}2f\big{)}\frac{\partial^{2}\lambda}{\partial x^{2}}+\alpha_{2,3}\big{(}3f^{2}\big{)}\frac{\partial^{2}\lambda}{\partial x^{2}}$ $\displaystyle-\alpha_{3,1}\frac{\partial^{3}\lambda}{\partial x^{3}}-\alpha_{3,2}\big{(}2f\big{)}\frac{\partial^{3}\lambda}{\partial x^{3}}-\alpha_{3,3}\big{(}3f^{2}\big{)}\frac{\partial^{3}\lambda}{\partial x^{3}}$ (20) $$

with final condition $\lambda(x,T)=2\big{(}f^{*}(x,T)-f(x,T)\big{)}$. The parameters $\bm{\alpha}$ are then found using the gradient descent method with update rule

$$ $\displaystyle\alpha_{d,p}\leftarrow\alpha_{d,p}-\eta\,\frac{\partial\mathcal{C}}{\partial\alpha_{d,p}}\qquad\text{where}\quad\eta=\beta\min(\Delta x)^{d-d_{\mathrm{max}}}\quad\text{and}\quad\frac{\partial\mathcal{C}}{\partial\alpha_{d,p}}=(-1)^{d}\int f^{p}\,\frac{\partial^{d}\lambda}{\partial x^{d}}dxdt\leavevmode\nobreak\ .$ (21) $$

This leads to the update rule for each coefficient

$$ $\begin{array}[]{l@{}l@{}l|l@{}l@{}l|l@{}l@{}l}\alpha_{1,1}&\leftarrow&\,\,\alpha_{1,1}+\displaystyle\frac{\beta}{\min(\Delta x)^{2}}\int f\,\frac{\partial\lambda}{\partial x}dxdt&\,\,\alpha_{2,1}&\leftarrow&\,\,\alpha_{2,1}-\displaystyle\frac{\beta}{\min(\Delta x)}\int f\,\frac{\partial^{2}\lambda}{\partial x^{2}}dxdt&\,\,\alpha_{3,1}&\leftarrow&\,\,\alpha_{3,1}+\beta\displaystyle\int f\,\frac{\partial^{3}\lambda}{\partial x^{3}}dxdt\[12.0pt] \alpha_{1,2}&\leftarrow&\,\,\alpha_{1,2}+\displaystyle\frac{\beta}{\min(\Delta x)^{2}}\int f^{2}\,\frac{\partial\lambda}{\partial x}dxdt&\,\,\alpha_{2,2}&\leftarrow&\,\,\alpha_{2,2}-\displaystyle\frac{\beta}{\min(\Delta x)}\int f^{2}\,\frac{\partial^{2}\lambda}{\partial x^{2}}dxdt&\,\,\alpha_{3,2}&\leftarrow&\,\,\alpha_{3,2}+\beta\displaystyle\int f^{2}\,\frac{\partial^{3}\lambda}{\partial x^{3}}dxdt\[12.0pt] \alpha_{1,3}&\leftarrow&\,\,\alpha_{1,3}+\displaystyle\frac{\beta}{\min(\Delta x)^{2}}\int f^{3}\,\frac{\partial\lambda}{\partial x}dxdt&\,\,\alpha_{2,3}&\leftarrow&\,\,\alpha_{2,3}-\displaystyle\frac{\beta}{\min(\Delta x)}\int f^{3}\,\frac{\partial^{2}\lambda}{\partial x^{2}}dxdt&\,\,\alpha_{3,3}&\leftarrow&\,\,\alpha_{3,3}+\beta\displaystyle\int f^{3}\,\frac{\partial^{3}\lambda}{\partial x^{3}}dxdt\leavevmode\nobreak\ .\end{array}$ $$

### A.2 Reaction Diffusion System of Equations

As mentioned in Section [3.5](#S3.SS5), for this case, we consider a system consisting of two PDEs, i.e. $N=\dim(\bm{f})=\dim(\bm{p})=2$ where $\bm{f}=[f_{1},f_{2}]$ and $\bm{p}=[p_{1},p_{2}]$, in a two-dimensional input space, i.e. $n=\dim(\bm{x})=\dim(\bm{d})=2$ where $\bm{x}=[x_{1},x_{2}]$, and $\bm{d}=[d_{1},d_{2}]$. In addition, we consider candidate terms with derivatives such that $\bm{d}\in\mathcal{D}_{\bm{d}}=\{[0,0],[1,0],[0,1],[2,0],[0,2]\}$ and polynomials such that $\bm{p}\in\mathcal{D}_{\bm{p}}=\{[1,0],[0,1],[1,1],[2,0],[0,2],[2,1],[1,2],[3,0],[0,3]\}$. In other words, $d_{\mathrm{max}}=2$ and $p_{\mathrm{max}}=3$. The resulting forward model in Eq. [1](#S2.E1) takes the form

$$ $\displaystyle\mathcal{L}_{i}[\bm{f}]=\partial_{t}f_{i}+\sum_{\bm{d},\bm{p}}\alpha_{i,\bm{d},\bm{p}}\nabla^{(\bm{d})}_{{\bm{x}}}[\bm{f}^{\bm{p}}]$ (22) $$

where $i\in\{1,2\}$, $\bm{f}^{\bm{p}}=f_{1}^{p_{1}}f_{2}^{p_{2}}$ and $\nabla^{(\bm{d})}_{{\bm{x}}}=\nabla^{(d_{1})}_{x_{1}}\nabla^{(d_{2})}_{x_{2}}$. This is equivalent to

$$ $\displaystyle\mathcal{L}_{i}[f_{1},f_{2}]=\frac{\partial f_{i}}{\partial t}$ $\displaystyle+\sum_{[d_{1},d_{2}]\in\mathcal{D}_{\bm{d}}}\sum_{[p_{1},p_{2}]\in\mathcal{D}_{\bm{p}}}\alpha_{i,[d_{1},d_{2}],[p_{1},p_{2}]}\frac{\partial^{\,d_{1}+d_{2}}\big{(}f_{1}^{p_{1}}f_{2}^{p_{2}}\big{)}}{\partial x_{1}^{d_{1}}\partial x_{2}^{d_{2}}}$ $\displaystyle=\frac{\partial f_{i}}{\partial t}$ $\displaystyle+\alpha_{i,[0,0],[1,0]}f_{1}+\alpha_{i,[0,0],[0,1]}f_{2}+\alpha_{i,[0,0],[1,1]}f_{1}f_{2}+\ldots+\alpha_{i,[0,0],[0,3]}f_{2}^{3}$ $\displaystyle+\alpha_{i,[1,0],[1,0]}\frac{\partial f_{1}}{\partial x_{1}}+\alpha_{i,[1,0],[0,1]}\frac{\partial f_{2}}{\partial x_{1}}+\alpha_{i,[1,0],[1,1]}\frac{\partial\big{(}f_{1}f_{2}\big{)}}{\partial x_{1}}+\ldots+\alpha_{i,[1,0],[0,3]}\frac{\partial\big{(}f_{2}^{3}\big{)}}{\partial x_{1}}$ $\displaystyle+\hskip 20.0pt\ldots$ $\displaystyle+\alpha_{i,[0,2],[1,0]}\frac{\partial^{2}f_{1}}{\partial x_{2}^{2}}+\alpha_{i,[0,2],[0,1]}\frac{\partial^{2}f_{2}}{\partial x_{2}^{2}}+\alpha_{i,[0,2],[1,1]}\frac{\partial^{2}\big{(}f_{1}f_{2}\big{)}}{\partial x_{2}^{2}}+\ldots+\alpha_{i,[0,2],[0,3]}\frac{\partial^{2}\big{(}f_{2}^{3}\big{)}}{\partial x_{2}^{2}}$ (23) $$

where $i\in\{1,2\}$. As we can observe, we have $\lvert\mathcal{D}_{\bm{d}}\rvert\times\lvert\mathcal{D}_{\bm{p}}\rvert=5\times 9=45$ terms with unknown coefficients $\bm{\alpha}_{i}=[\alpha_{i,\bm{d},\bm{p}}]_{d\in\mathcal{D}_{\bm{d}},p\in\mathcal{D}_{\bm{p}}}$ for the $i$-th PDE, i.e. a total of $90$ terms for the considered system, that we aim to find using the proposed adjoint method.

The cost functional in this case is simply

$$ $\displaystyle\mathcal{C}=\sum_{i=1}^{2}\bigg{(}$ $\displaystyle\int\big{(}f_{i}^{*}(\bm{x},t)-f_{i}(\bm{x},t)\big{)}^{2}d\bm{x}dt+\int\big{(}f_{i}^{*}(\bm{x},T)-f_{i}(\bm{x},T)\big{)}^{2}d\bm{x}+\int\lambda_{i}(\bm{x},t)\mathcal{L}_{i}[\bm{f}(\bm{x},t)]d\bm{x}dt\bigg{)}\leavevmode\nobreak\ .$ (24) $$

The corresponding adjoint equation is given by

$$ $\displaystyle\frac{\partial\lambda_{i}}{\partial t}=-2(f^{*}_{i}-f_{i})$ $\displaystyle+\sum_{\bm{d},\bm{p}}(-1)^{|\bm{d}|}\alpha_{i,\bm{d},\bm{p}}\nabla_{f_{i}}[\bm{f}^{\bm{p}}]\nabla^{(\bm{d})}_{{\bm{x}}}[\lambda_{i}]$ $\displaystyle=-2(f^{*}_{i}-f_{i})$ $\displaystyle+\sum_{[d_{1},d_{2}]\in\mathcal{D}_{\bm{d}}}\sum_{[p_{1},p_{2}]\in\mathcal{D}_{\bm{p}}}(-1)^{d_{1}+d_{2}}\alpha_{i,[d_{1},d_{2}],[p_{1},p_{2}]}\frac{\partial\big{(}f_{1}^{p_{1}}f_{2}^{p_{2}}\big{)}}{\partial f_{i}}\frac{\partial^{\,d_{1}+d_{2}}\lambda_{i}}{\partial x_{1}^{d_{1}}\partial x_{2}^{d_{2}}}$ (25) $$

and $\lambda_{i}(\bm{x},T)=2\big{(}f^{*}_{i}(\bm{x},T)-f_{i}(\bm{x},T)\big{)}$ where $i\in\{1,2\}$.

Assume, without loss of generality, that $i=1$. Then, we can write

$$ $\displaystyle\frac{\partial\lambda_{1}}{\partial t}=-2(f^{*}_{1}-f_{1})$ $\displaystyle+\alpha_{1,[0,0],[1,0]}\lambda_{1}+\alpha_{1,[0,0],[1,1]}f_{2}\lambda_{1}+\alpha_{1,[0,0],[2,0]}\big{(}2f_{1}\big{)}\lambda_{1}+\ldots+\alpha_{1,[0,0],[3,0]}\big{(}3f_{1}^{2}\big{)}\lambda_{1}$ $\displaystyle-\alpha_{1,[1,0],[1,0]}\frac{\partial\lambda_{1}}{\partial x_{1}}-\alpha_{1,[1,0],[1,1]}f_{2}\frac{\partial\lambda_{1}}{\partial x_{1}}-\alpha_{1,[1,0],[2,0]}\big{(}2f_{1}\big{)}\frac{\partial\lambda_{1}}{\partial x_{1}}-\ldots-\alpha_{1,[1,0],[3,0]}\big{(}3f_{1}^{2}\big{)}\frac{\partial\lambda_{1}}{\partial x_{1}}$ $\displaystyle+\hskip 20.0pt\ldots$ $\displaystyle+\alpha_{1,[0,2],[1,0]}\frac{\partial^{2}\lambda_{1}}{\partial x_{2}^{2}}+\alpha_{1,[0,2],[1,1]}f_{2}\frac{\partial^{2}\lambda_{1}}{\partial x_{2}^{2}}+\alpha_{1,[0,2],[2,0]}\big{(}2f_{1}\big{)}\frac{\partial^{2}\lambda_{1}}{\partial x_{2}^{2}}+\ldots+\alpha_{1,[0,2],[3,0]}\big{(}3f_{1}^{2}\big{)}\frac{\partial^{2}\lambda_{1}}{\partial x_{2}^{2}}$ (26) $$

and $\lambda_{1}(\bm{x},T)=2\big{(}f^{*}_{1}(\bm{x},T)-f_{1}(\bm{x},T)\big{)}$. We can follow the same procedure for $i=2$. The parameters $\bm{\alpha}_{i}$ are then found using the gradient descent method with update rule

$$ $\displaystyle\alpha_{i,\bm{d},\bm{p}}\leftarrow\alpha_{i,\bm{d},\bm{p}}-\eta\frac{\partial\mathcal{C}}{\partial\alpha_{i,\bm{d},\bm{p}}}\qquad\text{where}\quad\eta=\beta\min(\Delta\bm{x})^{|\bm{d}|-d_{\mathrm{max}}}\quad\text{and}\quad\frac{\partial\mathcal{C}}{\partial\alpha_{i,\bm{d},\bm{p}}}=(-1)^{|\bm{d}|}\int\bm{f}^{\bm{p}}\,\nabla^{(\bm{d})}_{{\bm{x}}}[\lambda_{i}]d\bm{x}dt$ (27) $$

leading to the update rule for each coefficient

$$ $\begin{array}[]{l@{}l@{}l|c|l@{}l@{}l}\alpha_{i,[0,0],[1,0]}&\leftarrow&\,\,\alpha_{i,[0,0],[1,0]}-\displaystyle\frac{\beta}{\min(\Delta\bm{x})^{2}}\int f_{1}\lambda_{i}d\bm{x}dt&\qquad\cdots&\alpha_{i,[0,0],[0,3]}&\leftarrow&\,\,\alpha_{i,[0,0],[0,3]}-\displaystyle\frac{\beta}{\min(\Delta\bm{x})^{2}}\int f_{2}^{3}\lambda_{i}d\bm{x}dt\[12.0pt] \alpha_{i,[1,0],[1,0]}&\leftarrow&\,\,\alpha_{i,[1,0],[1,0]}+\displaystyle\frac{\beta}{\min(\Delta\bm{x})}\int f_{1}\frac{\partial\lambda_{i}}{\partial x_{1}}d\bm{x}dt&\qquad\cdots&\alpha_{i,[1,0],[0,3]}&\leftarrow&\,\,\alpha_{i,[1,0],[0,3]}+\displaystyle\frac{\beta}{\min(\Delta\bm{x})}\int f_{2}^{3}\frac{\partial\lambda_{i}}{\partial x_{1}}d\bm{x}dt\[12.0pt] &\hskip 6.0pt\vdots&&\cdots&&\hskip 6.0pt\vdots&\[12.0pt] \alpha_{i,[0,2],[1,0]}&\leftarrow&\,\,\alpha_{i,[0,2],[1,0]}-\displaystyle\beta\int f_{1}\frac{\partial^{2}\lambda_{i}}{\partial x_{2}^{2}}d\bm{x}dt&\qquad\cdots&\alpha_{i,[0,2],[0,3]}&\leftarrow&\,\,\alpha_{i,[0,2],[0,3]}-\displaystyle\beta\int f_{2}^{3}\frac{\partial^{2}\lambda_{i}}{\partial x_{2}^{2}}d\bm{x}dt\leavevmode\nobreak\ .\end{array}$ $$

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: Hopefully, the authors can make your contribution clearly and make your experiment and theories to test your results.
soundness: 2
strengths: The topic looks very interesting.
summary: The authors claim that they can discover some PDE via an adjoint-based method.
weaknesses: The manuscript has not been well-written, so that the reader cannot find their motivation clearly.
The theoretical part has been shown rigorously.
The experiment part: description not clear

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 2
questions: (1) Can you provide more complex examples in high dimensional space?
(2) Can you provide some a posteriori analysis analysis regarding the learning performance with respect to the given data or type of PDEs?
(3) Can you derive some well posed-ness of the parametrization of PDE? i.e, if I provide two similar parameterized PDEs, how about the learning performances regarding interpolation and prediction?
soundness: 2
strengths: The paper is generally well-written, with clear explanations of the methodology and a solid foundation in the adjoint approach. The derivations are technically sound and should be accessible to readers familiar with inverse problems and PDEs.
summary: This paper presents a method for reconstructing a parameterized partial differential equation (PDE) from observed data using the adjoint method. The authors claim that their approach offers a new perspective on PDE reconstruction by leveraging adjoint-based optimization, demonstrating the method on simple examples and show better performance compared to the well known PDE-FIND.
weaknesses: (1) The approach presented here is relatively incremental, as the reconstruction of parameterized PDEs has been extensively studied within the inverse problem community for several decades. Many of the concepts explored, particularly adjoint-based parameter estimation, are already well-established. The paper would benefit from a more explicit discussion of how this work advances or differs from existing methods in the literature on PDE-constrained optimization and inverse problems.
(2) The demonstration example is limited to simple lower dimensional problems, which may not be sufficient to convincingly illustrate the method’s robustness or scalability to higher-dimensional, real-world PDEs. Given the computational efficiency implied by the adjoint method, testing on a more challenging, higher-dimensional example or a PDE with more complex dynamics would strengthen the paper. This would also allow for a more thorough evaluation of the method's effectiveness in a broader range of realistic scenarios.
(3) The paper would benefit from a comparative analysis with other state-of-the-art methods for PDE reconstruction. This comparison could help clarify the unique contributions and limitations of the adjoint-based approach relative to current methods in data-driven PDE reconstruction.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 1
questions: See Weakness above. We include here some additional minor issues:

1. When considering sensitivity to noise (line 411), the noise variance is specified in percentage. With respect to what? Is it an SNR?
1. Where does the normalization of the integral in (3) come from? This doesn’t seem to make sense given that the manuscript considers a regular grid.
1. When discussing the learning rate, it is mentioned that “the gradient of the cost function is most sensitive to the highest order terms of the PDE.” Why is that the case?
1. When defining $f^*$ on line 121, shouldn’t it be $f^*: \Omega \times [0,T] \to \mathbb{R}^N$? As it is, the desired solution is only on $\mathcal{G}$. How would that be consistent with, e.g., coarser time grids?
soundness: 2
strengths: - The proposed method archives good results on the considered PDEs considered compared to PDE-FIND
- The experiments consider relevant settings, including coarse temporal resolution, noisy data, and ill-posed tasks
- Using the adjoint method, it is possible find straightforward algorithms based on gradient descent
summary: This paper considers the problem of identifying the parameters of a system of PDEs based on measurement data. The proposed method fits the solution of the PDE to the data while adapting its parameters. It uses the adjoint method that combines a data goodness-of-fit objective with a PDE satisfaction constraint using a Lagrangian formulation. Using variational techniques, it is then possible to obtain evolution equations for those weights (multipliers), which then allow the gradient objective to be computed with respect to the PDE parameters. The proposed method is compared with PDE-FIND for a number of PDEs in the standard setting as well as when the temporal data is coarse and when the data is noisy.
weaknesses: 1. The algorithms proposed in the work implicitly (i) solve the forward PDE (1), which can be challenging and time consuming (line 242); (ii) solve the adjoint PDE model (7), which again can be challenging and time consuming (line 243); and (iii) compute gradients of the adjoint variables $\lambda$, which are only know implicitly through (7) (line 244). The manuscript does not thoroughly discuss how this affects the computational complexity of the proposed method and whether the choice of solution algorithm could affect the performance of the proposed method. This is particularly important since other methods, such as PDE-FIND, do not appear to require such explicit forward evolutions.
It is not clear from the manuscript how the proposed method differs from classical adjoint methods, proposed originally in the context of design or shape optimization (see, e.g., Jameson 2003). The text only mentions that “unlike the usual use of PDE-constrained adjoint optimization where the governing equation is known, in this paper we are interested in finding the form along with the coefficients of the PDE given data.” Nevertheless, the parameters of the PDE can be mapped directly to design parameters used in traditional adjoint methods. As such, there doesn’t seem to be any challenge requiring substantial contributions in the setting of this paper.

1. The paper is not clearly written. There are several typos (“Eq. equation” in lines 201, 215, 376, 482, 1172) and odd uses of language (“Here, $x(k)$ is coordinates”; ) that make it hard to follow. Many elements are either undefined or defined implicitly/informally. For instance, in (2) are those gradient products or an iterated differentiation (as in multi-index differentials from Sobolev spaces)? What is meant by “semi-discrete total variation of C” (line 176)? The definition is given without a reference or explanation. Additionally, the structure of the experiments is not clear as Section 3 named “Results” is followed by Section 4,5,6 named “Partial observations in time”... all containing experimental results as well. This lack of care in the writing renders the paper confusing for the reader.

1. Many results are provided without derivations or references. Since their definitions are incomplete, it is impossible to judge their validity. For instance, the derivations of (4) and (6) are only informally explained. And since nowhere is it defined in which functional space in which the solution $f$ of (1) must be in, it is impossible to judge the correctness of (4), even if the reader had to assume the smoothness/differentiability of $\lambda$. Additionally, since conditions such as $\lambda(x,t) \to 0$ for $x \in \partial\Omega$ are not justified and without knowing the functional (metric) space in which $\lambda$ is optimized, it is impossible to guarantee its differentiability in (4) and (6). Similar comments apply to the derivation of (7)-(8) from (6).

1. The manuscript only considers very simple PDEs. In particular, the main text only tackles the 1D heat equation with D=1, the 1D Burgers equation, and a simple 1D wave equation. These are low-dimensional, simple first-order equations with very few parameters to fit. In fact, the proposed method does not appear to be able to handle PDEs beyond linear, first order in time (particularly since it must solve (1) explicitly). In particular, it does not consider very overparametrized models either (large $p$ or $d$). Hence, the performance of the method in PDEs whose solutions are challenging or where there is substantial uncertainty on the model (large $p$ or $d$) is not tested.

1. Point (4) is particularly critical given that the proposed method is only compared against one baseline (PDE-FIND) which is already somewhat dated (2017). Many other baselines are available, in particular, those detailed in the introduction. It is unclear why none were used in the experiments.

1. There is only very limited discussion on the use of l2 regularization. Particularly since it is well-known from polynomial fitting that the l2-norm is not well-adapted to find sparse solutions, i.e., where many values are 0. That is why, e.g., Brunton et al. 2016 use the l1 norm. While the thresholding heuristic employed by the paper is reminiscent of iterative hard thresholding (IHT) methods from compressed sensing (see, e.g., Blumensath and Davis, Applied Computational Harmonic Analysis, 2009) it does not inherit the same guarantees or theoretical advantages. In fact, this relation is not even mentioned in the manuscript.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: Please see below questions on the paper content:
1. line 87. PINN-SR uses the Adam optimizer, which can process data in batches. Thus, why do authors feel that PINN-SR will not scale well with the size of the data set?  
2. line 135. The vector $\boldsymbol{f}^p$ may be missing a comma.
3. line 189. Does the framework work only for zero boundary conditions? 
4. line 218. How $\sigma$ is selected. Consider an example of the Navier-Stokes equation, where the viscosity can be in the order of $10^{-4}$ to $10^{-5}$. In such cases, the proposed algorithm will discard the relevant terms.
5. line 282. Is this method applicable to time-independent systems like Darcy and Poisson's equations? If applicable, how averaging of gradients in algorithm 2 will be done in the absence of the time component.
6. Fig. 1(c). It is intriguing to see that the error increases with an increase in data in both methods.
7. Fig. 1(d). Since the proposed method uses a forward solver in the loop, is it not the computational time supposed to increase with increasing training samples?
8. Eq. (13-14). The authors show that even in ill-posed problems, the proposed method can provide an alternate sparse equation. However, estimating the error and visually observing the solution fields before concluding on the accuracy would be best. 

**Limitations:**
The proposed framework does not show discovery from irregular observations (both in time and space). This should be included in the limitations.
soundness: 3
strengths: 1. The idea of using adjoint equations to uncover underlying equations is interesting, and it seems to be new in the literature. 
2. The theoretical derivations and analysis do a good job of explaining the key concepts of the framework. 
3. Even in partial observations, when data at fine mesh is not available, the adjoint equations can discretize the general PDE on a finer mesh and, therefore, outperform other methods in the low data limit.
summary: An adjoint-based equation discovery framework is proposed in this paper to discover the underlying governing partial differential equations (PDEs) of physical systems. The proposed method initializes a general PDE by considering derivatives and polynomials of the derivatives up to a finite degree. A loss function is formed by incorporating the error between the prediction from the general PDE and the training labels and the corresponding adjoint equations. After that, a PDE-constrained optimization problem is solved to estimate the parameters of the initial PDE model, such that the solution of the general PDE will accurately approximate the training labels. The successful optimization yields coefficients of only the relevant derivatives to take significant values, thereby uncovering the underlying governing equation.
weaknesses: Please see below the comments that need to be further clarified. 
1. The proposed framework is motivated by assuming a general PDE, which contains derivatives and their polynomials up to a certain degree. This is equivalent to basis functions in regression-based equation discovery frameworks, which is one of the major limitations of the basis-dependent discovery methods. Since, in an unknown scenario, the knowledge about the underlying physics will be minimal, one may need to consider a large number of basses/derivatives. In most of the examples, the authors consider only up to 3rd-order derivatives and their polynomials. Authors should consider a higher order of derivatives and polynomials to check the performance of their framework, both accuracy and computational efficiency.
2. The proposed framework closely resembles the PINN-SR algorithm. Therefore, a comparison with PINN-SR algorithms is necessary to gauge the effect of the adjoint equations.
3. The motivation for the equation discovery architecture is a little convoluted. The authors state that since data-driven simulators fail to learn the exact physics of dynamical systems, they fail to extrapolate beyond the training regime. Thus, we need an equation discovery framework, which post-discovery can be coupled with any standard numerical methods to obtain accurate predictions. While the second statement is true, it can be argued that after the discovery of a governing equation, one can still train a data-driven or physics-informed machine learning emulator for accelerating computational simulations. Therefore, the motivation needs to be reworked.
4. The Bayesian class of equation discovery algorithms also does a good job in distilling governing equation equations from data, particularly in noisy and low-data limits [1-3]. This paper should have discussed these frameworks. A comparison with a few such frameworks would further benefit the content of this paper.

[1] Nayek, Rajdip, et al. "On spike-and-slab priors for Bayesian equation discovery of nonlinear dynamical systems via sparse linear regression." Mechanical Systems and Signal Processing 161 (2021): 107986.\
[2] Zhang, Zhiming, and Yongming Liu. "Parsimony-enhanced sparse Bayesian learning for robust discovery of partial differential equations." Mechanical Systems and Signal Processing 171 (2022): 108833.\
[3] Hirsh, Seth M., David A. Barajas-Solano, and J. Nathan Kutz. "Sparsifying priors for Bayesian uncertainty quantification in model discovery." Royal Society Open Science 9.2 (2022): 211823.