Feature engineering and symbolic regression methods for detecting hidden physics from sparse sensor observation data

Since the dawn of mathematical modeling of complex physical processes, scientists have been attempting to formulate predictive models to infer current and future states. These first principle models are generally conceptualized from conservation laws, sound physical arguments, and empirical heuristics drawn from either conducting experiments or hypotheses made by an insightful researcher. However, there are many complex systems (some being climate science, weather forecasting, and disease control modeling) with their governing equations known partially, and their hidden physics wait to be modeled. In the last decade, there have been rapid advances in machine learning^1,2 and easy access to rich data, thanks to the plummeting costs of sensors and high performance computers.

This paradigm shift in data driven techniques can be readily exploited to distill new or improved physical models for nonlinear dynamical systems. Extracting predictive models based on observing complex patterns from vast multimodal data can be loosely termed reverse engineering nature. This approach is not particularly new; for example, Kepler used planets’ positional data to approximate their elliptic orbits. The reverse engineering approach is most appropriate in the modern age as we can leverage computers to directly infer physical laws from data collected from omnipresent sensors that otherwise might not be comprehensible to humans. Symbolic regression methods are a class of data driven algorithms that aim to find a mathematical model that can describe and predict hidden physics from observed input-response data. Some of the popular machine learning techniques that are adapted for the task of symbolic regression are neural networks,^3,4 compressive sensing/sparse optimization,^5,6 and evolutionary algorithms.^7,8

Symbolic regression (SR) approaches based on evolutionary computation^7,9 are a class of frameworks that are capable of finding analytically tractable functions. Traditional deterministic regression algorithms assume a mathematical form and only find parameters that best fit the data. On the other hand, evolutionary SR approaches aim to simultaneously find parameters and also learn the best-fit functional form of the model from input-response data. Evolutionary algorithms search for functional abstractions with a preselected set of mathematical operators and operands while minimizing the error metrics. Furthermore, the optimal model is selected from Pareto front analysis with respect to minimizing accuracy vs model complexity. Genetic programming (GP)⁷ is a popular choice leveraged by most of the SR frameworks. GP is an extended and improved version of a genetic algorithm (GA),^10,11 which is inspired by Darwin’s theory of natural evolution. Seminal work was done in identifying hidden physical laws^12,13 from the input-output response using the GP approach. GP has been applied in the context of the SR approach in digital signal processing,¹⁴ nonlinear system identification,¹⁵ and aerodynamic parametric estimation.¹⁶ Furthermore, GP as an SR tool was applied to identify complex closed-loop feedback control laws for turbulent separated flows.^17–20 Hidden physical laws of the evolution of a harmonic oscillator based on sensor measurements and the real world prediction of solar power production at a site were identified using GP as an SR approach.²¹ 

Improved versions of GP focus on better representation of the chromosome, which helps in the free evolution of the chromosome with constraints on the complexity of its growth, and faster searches for the best chromosome. Some of these improved versions of GP are gene expression programming (GEP),⁸ parse matrix evolution (PME),²² and linear genetic programming (LGP).²³ GEP takes advantage of the linear coded chromosome approach from GA and the parse tree evolution of GP to alleviate the disadvantages of both GA and GP. GEP was applied to diverse applications as an SR tool to recover nonlinear dynamical systems.^24–27 Recently, GEP was modified for tensor regression, termed as multi-GEP, and has been applied to recover functional models approximating the nonlinear behavior of the stress tensor in the Reynolds-averaged Navier-Stokes (RANS) equations.²⁸ Furthermore, this novel algorithm was extended to identify closure models in a combustion setting for large eddy simulations (LES).²⁹ Similarly, a new damping function has been discovered using the GEP algorithm for the hybrid RANS/LES methodology.³⁰ Generally, evolutionary based SR approaches can identify models with complex nonlinear compositions given enough computational time.

Compressive sensing (CS)^5,6 is predominately applied to signal processing in seeking the sparsest solution (i.e., a solution with the fewest number of features). Basis pursuit algorithms,³¹ also identified as sparsity promoting optimization techniques,^32,33 play a fundamental role in CS. Ordinary least squares (OLS) optimization generally results in identifying models with large complexity, which are prone to overfitting. In sparse optimization, the OLS objective function is regularized by an additional constraint on the coefficient vector. This regularization helps in taming and shrinking large coefficients and thereby promoting sparsity in feature selection and avoiding overfitted solutions. The least absolute shrinkage and selection operator (LASSO)^32,34 is one of the most popular regularized least squares (LS) regression methods. In LASSO, an L₁ penalty is added to the LS objective function to recover sparse solutions.³⁵ In Bayesian terms, LASSO is a maximum a posteriori estimate (MAP) of LS with Laplacian priors. LASSO performs feature selection and simultaneously shrinks large coefficients, which may manifest to overfit the training data. Ridge regression³⁶ is another regularized variant where an L₂ penalty is added to the LS objective function. Ridge regression is also defined as a MAP estimate of LS with a Gaussian prior. The L₂ penalty helps in grouping multiple correlated basis functions and increases robustness and convergence stability for ill-conditioned systems. The elastic net approach^37,38 is a hybrid of the LASSO and ridge approaches combining the strengths of both algorithms.

Derived from these advances, a seminal work was done in employing sparse regression to identify the physical laws of nonlinear dynamical systems.³⁹ This work leverages the structure of sparse physical laws, i.e., only a few terms represent the dynamics. The authors have constructed a large feature library of potential basis functions that has the expressive power to define the dynamics and then seek to find a sparse feature set from this overdetermined system. To achieve this, a sequential threshold least squares (STLS) algorithm³⁹ has been introduced in such a way that a hard threshold on OLS coefficients is performed recursively to obtain sparse solutions. This algorithm was leveraged to form a framework called sparse identification of nonlinear dynamics (SINDy)³⁹ to extract the physical laws of nonlinear dynamical systems represented by ordinary differential equations (ODEs). This work re-envisioned model discovery from the perspective of sparse optimization and compressive sensing. The SINDy framework recovered various benchmark dynamical systems such as the chaotic Lorenz system and vortex shedding behind a cylinder. However, STLS regression finds it challenging to discover physical laws that are represented by spatiotemporal data or high-dimensional measurements and have highly correlated features in the basis library. This limitation was addressed using a regularized variant of STLS called the sequential threshold ridge regression (STRidge) algorithm.⁴⁰ This algorithm was intended to discover unknown governing equations that are represented by partial differential equations (PDEs), hence forming a framework termed as PDE-functional identification of nonlinear dynamics (PDE-FIND).⁴⁰ PDE-FIND was applied to recover canonical PDEs representing various nonlinear dynamics. This framework also performs reasonably well under the addition of noise to data and measurements. These sparse optimization frameworks generally have a free parameter associated with the regularization term that is tuned by the user to recover models ranging from highly complex to parsimonious.

In a similar direction of discovering governing equations using sparse regression techniques, L₁ regularized LS minimization was used to recover various nonlinear PDEs^41,42 using both high fidelity and distorted (noisy) data. Additionally, limited and distorted data samples were used to recover chaotic and high-dimensional nonlinear dynamical systems.^43,44 To automatically filter models with respect to model complexity (number of terms in the model) vs test accuracy, Bayes information criteria were used to rank the most informative models.⁴⁵ Furthermore, SINDy coupled with model information criteria is used to infer canonical biological models⁴⁶ and introduce a reduced order modeling (ROM) framework.⁴⁷ STRidge⁴⁰ was applied as a deterministic SR method to derive algebraic Reynolds-stress models for the RANS equations.⁴⁸ Recently, various sparse regression algorithms such as LASSO,³² STRidge,⁴⁰ sparse relaxed regularized regression,⁴⁹ and the forward-backward greedy algorithm⁵⁰ were investigated to recover truncation error terms of various modified differential equations (MDEs) coming from canonical PDEs.⁵¹ The frameworks discussed above assume that the structure of the model to be recovered is sparse in nature; that is, only a small number of terms govern the dynamics of the system. This assumption holds for many physical systems in science and engineering.

Fast function extraction (FFX)⁵² is another deterministic SR approach based on pathwise regularized learning that is also called the elastic net algorithm.³⁷ The resulting models of FFX are selected through nondominated filtering concerning accuracy and model complexity, similar to evolutionary computations. FFX is influenced by both GP and CS to better distill physical models from data. FFX has been applied to recover hidden physical laws,²¹ canonical governing equations,⁵³ and Reynolds stress models for the RANS equations.⁵⁴ Some other potential algorithms for deterministic SR are elite bases regression (EBR)⁵⁵ and prioritized grammar enumeration (PGE).⁵⁶ EBR uses only elite features in the search space selected by measuring correlation coefficients of features for the target model. PGE is another deterministic approach that aims for the substantial reduction of the search space where the genetic operators and random numbers from GP are replaced with grammar production rules and systematic choices.

An artificial neural network (ANN), also referred to as deep learning if multiple hidden layers are used, is a machine learning technique that transforms input features through nonlinear interactions and maps to output target features.^3,4 ANNs attracted attention in recent times due to their exemplary performance in modeling complex nonlinear interactions across a wide range of applications including image processing,⁵⁷ video classification,⁵⁸ and autonomous driving.⁵⁹ ANNs produce black-box models that are not quite open to physical inference or interpretability. Recently, physics-informed neural networks (PINNs)⁶⁰ were proposed in the flavor of SR that is capable of identifying scalar parameters for known physical models. PINNs use a loss function in symbolic form to help ANNs adhere to the physical structure of the system. Along similar directions, a Gaussian process regression (GPR) has been also investigated for the discovery of coefficients by recasting unknown coefficients as GPR kernel hyperparameters for various time dependent PDEs.^61,62 As a nonlinear system identification tool, the GPR approach provides a powerful framework to model dynamical systems.^63,64 State calibration with the four dimensional variational data assimilation (4D VAR)⁶⁵ and deep learning techniques such as long short-term memory (LSTM)⁶⁶ have been used for model identification in ROM settings. Convolutional neural networks (CNNs) are constructed to produce hidden physical laws from using the insight of establishing direct connections between filters and finite difference approximations of differential operators.^67,68 This approach has been demonstrated to discover underlying PDEs from learning the filters by minimizing the loss functions.^69,70

In this paper, we have exploited the use of SR in three different applications: equation discovery, truncation error analysis, and hidden physics discovery. We demonstrate the use of the evolutionary computation algorithm, GEP, and the sparse regression algorithm, STRidge, in the context of the SR approach to discover various physical laws represented by linear and nonlinear PDEs from observing input-response data. We begin by demonstrating the identification of canonical linear and nonlinear PDEs that are up to the fifth order in space. For identifying one particular PDE, we demonstrate the natural feature extraction ability of GEP and the limits in the expressive and predictive power of using a feature library when dealing with STRidge in discovering physical laws. We then demonstrate the discovery of highly nonlinear truncation error terms of the Burgers MDE using both GEP and STRidge. We highlight that the analysis of truncation errors is very important in the implicit large eddy simulation as a way to determine inherent turbulence models. This analysis is usually very tedious and elaborate, and our study provides a clear example of how SR tools are suitable in such research. Following truncation error terms identification, we apply GEP using sparse data to recover hidden source terms represented by complex function compositions for a one-dimensional (1D) advection-diffusion process and a two-dimensional (2D) vortex-merger problem. Furthermore, both GEP and STRidge are used to demonstrate the identification of the eddy viscosity kernel along with its ad hoc modeling coefficient closing LES equations simulating the 2D decaying turbulence problem. An important result is the ability of the proposed methodology to distill the Smagorinsky model from an array of tailored features in solving the Kraichnan turbulence problem.

The rest of the paper is organized as follows: Section II gives a brief description of the GEP and STRidge algorithms. In Sec. III, GEP, and STRidge are tested on identifying different canonical PDEs. Section IV deals with the identification of nonlinear truncation terms of the Burgers MDE using both STRidge and GEP. In Sec. V we exploit GEP for identification of hidden source terms in a 1D advection-diffusion process and a 2D vortex-merger problem. We additionally demonstrate recovery of the eddy viscosity kernel and its modeling coefficient by both GEP and STRidge for closing the LES equations simulating the 2D decaying turbulence problem in the same section. Finally, Sec. VI draws our conclusions and highlights some ideas for future extensions of this work.

We recover various physical models from data using two symbolic regression tools, namely, GEP, an evolutionary computing algorithm, and STRidge, which is a deterministic algorithm that draws its influences from compressive sensing and sparse optimization. We take the example of the equation discovery problem that is discussed in Sec. III to elaborate on the methodology of applying GEP and STRidge for recovering various physical models. We restrict the PDEs to be recovered to quadratic nonlinear and up to the fifth order in space. The general nonlinear PDE to be recovered is in the form of

where subscripts denote the order of partial differentiation and σ is an arbitrary parameter. For example, consider the problem of identifying the viscous Burgers equation

where $u(x,t)∈Rm×n$ is the velocity field and ν is the kinematic viscosity. In our study, m is the number of time snapshots and n is the number of spatial locations. The solution field u(x, t) is generally obtained by solving Eq. (2) analytically or numerically. The solution field might also be obtained from sensor measurements that can be arranged as shown follows:

For recovering PDEs, we need to construct a library of basis functions called feature library that contains higher order derivatives of the solution field u(x, t). Higher order spatial and temporal partial derivative terms can be approximated using any numerical scheme once the recording of the discrete data set given by Eq. (3) is available. In our current setup, we use the leapfrog scheme for approximating the temporal derivatives and central difference schemes for spatial derivatives as follows:

where temporal and spatial steps are given by dt and dx, respectively. Within the expressions presented in Eq. (4), the spatial location is denoted using subscript index j, and the temporal instant using superscript index p.

We note that other approaches such as automatic differentiation or spectral differentiation for periodic domains can easily be adopted within our study. Both GEP and STRidge take the input library consisting of features (basis functions) that are built using Eqs. (2) and (3). This core library, used for the equation discovery problem in Sec. III, is shown as follows:

The solution u(x, t) and its spatial and temporal derivatives are arranged with size m · n × 1 in each column of Eq. (5). For example, the features (basis functions) U and U_2x are arranged as follows:

where subscript j denotes the spatial location and subscript p denotes the time snapshot. The features (basis functions) in the core library $Θ̃(U)$ are expanded to include interacting features limited to quadratic nonlinearity and also a constant term. The final expanded library is given as

where the size of the library is $Θ(U)∈Rm⋅n×Nβ$ and N_β is the number of features (basis functions, i.e., N_β = 28 for our setup). For example, if we have 501 spatial points and 101 time snapshots with 28 bases, then Θ(U) [Eq. (7)] contains 501 × 101 rows and 28 columns.

Note that the core feature library $Θ̃(U)$ in Eq. (5) is given as an input to GEP to recover PDEs and the algorithm extracts higher degree nonlinear interactions of core features in $Θ̃(U)$ automatically. However, for sparse optimization techniques such as STRidge, explicit input of all possible combinations of core features in Eq. (5) is required. Therefore, Θ(U) in Eq. (7) forms the input to the STRidge algorithm for equation identification. This forms the fundamental difference in terms of feature building for both algorithms. Subsection II A gives a brief introduction to GEP and its specific hyperparameters that control the efficacy of the algorithm in identifying physical models from observing data. Furthermore, Subsection II B describes how to form linear system representations in terms of V(t) and Θ(U) and briefly describes the STRidge optimization approach to identifying sparse features and thereby building parsimonious models using spatiotemporal data.

Gene expression programming (GEP)^8,71 is a genotype-phenotype evolutionary optimization algorithm which takes advantage of the simple chromosome representation of genetic algorithm (GA)¹⁰ and the free expansion of complex chromosomes of genetic programming (GP).⁷ As in most evolutionary algorithms, this technique also starts with generating initial random populations, iteratively selecting candidate solutions according to a fitness function, and improving candidate solutions by modifying through genetic variations using one or more genetic operators. The main difference between GP and GEP is how both techniques define the nature of their individuals. In GP, the individuals are nonlinear entities of different sizes and shapes represented as parse trees, and in GEP, the individuals are encoded as linear strings of fixed length called genomes and chromosomes, similar to GA representation of individuals, and later expressed as nonlinear entities of different size and shape called phenotypes or expression trees (ET). GEP is used for a very broad range of applications, but here it is introduced as a symbolic regression tool to extract constraint-free solutions from input-response data.

The arrangement of a typical gene/chromosome in GEP is shown in Fig. 1. The GEP gene is composed of head and tail regions as illustrated in Fig. 1. The head of a gene consists of both symbolic terms from functions (elements from a function set F) and terminals (elements from a terminal set T) whereas the tail consists of only terminals. The function set F may contain arithmetic mathematical operators (e.g., +, ×, −, /), nonlinear functions (e.g., sin, cos, tan, arctan, sqrt, and exp), or Boolean operators (e.g., not, nor, or, and and) and the terminal set contains the symbolic variables. The gene always starts with a randomly generated mathematical operator from the function set F. The head length is one of the important hyperparameters of GEP, and it is determined using trial and error as there is no definite method to assign it. Once the head length is determined, the size of the tail is computed as a function of the head length and the maximum arity of a mathematical operator in the function set F.⁹ It can be calculated by the following equation:

where a_max is the maximum argument of a function in F. The single gene can be extended to multigenic chromosomes where individual genes are linked using a linking function (e.g., +, ×, /, −). The general rule of thumb is to have a larger head and higher number of genes when dealing with complex problems.⁹ 

The structural organization of the GEP gene is arranged in terms of open reading frames (ORFs) inspired from biology where the coding sequence of a gene equivalent to an ORF begins with a start codon, continues with an amino acid codon, and ends with a termination codon. In contrast to a gene in biology, the start site is always the first position of a gene in GEP, but the termination point does not always coincide with the last position of a gene. These regions of the gene are termed noncoding regions downstream of the termination point. Only the ORF region is expressed in the ET and can be clearly seen in Fig. 1.

Even though the noncoding regions in GEP genes do not participate in the final solution, the power of GEP evolvability lies in this region. The syntactically correct genes in GEP evolve after modification through diverse genetic operators due to this region’s chromosome. This is the paramount difference between GEP and GP implementations, where in the latter, many syntactically invalid individuals are produced and need to be discarded while evolving the solutions and additional special constraints are imposed on the depth/complexity of the candidate solution to be evolved to avoid bloating problem.¹⁹ 

Figure 2 displays the typical flowchart of the GEP algorithm. The process is described briefly below:

1. The optimization procedure starts with a random generation of chromosomes built upon combinations of functions and terminals. The size of the random population is a hyperparameter and the larger the population size, better the probability of finding the best candidate solution.

2. After the population is generated, the chromosomes are expressed as ETs, and then this population is converted to a numerical expression. This expression is then evaluated using a fitness function. In our setup, we employ the mean squared error between the best predicted model f^* and the true model f as the fitness function given by

where $flk*$ is the value predicted by the chromosome k for the fitness case l (out of N samples cases) and f_l is the true or measurement value for the l^th fitness case.

1. The termination criteria are checked after all fitness evaluations, to continue evolving or to save the best fitness chromosome as our final predicted model. In our current setup, we terminate after a specified number of generations.

2. The evolvability/reproduction of the chromosome through genetic operators, which is the core part of the GEP evolutionary algorithm, executes if termination criteria are not met. Before the genetic operations on the chromosome begin, the best chromosome according to the fitness function is cloned to the next generations using a selection method. Popular selection methods include tournament selection with elitism and roulette-wheel selection with elitism. In our current setup, we use tournament selection with elitism.

3. The four genetic operators that introduce variation in populations are mutation, inversion, transposition, and recombination. The GEP transposition operator is applied to the elements of the chromosome in three ways: insertion sequence (IS), root insertion sequence (RIS), and gene insertion sequence, and similarly, three kinds of recombination are applied, namely, one point, two point, and gene recombination.

4. The process is continued until the termination criteria are met, which is the number of generations in our current setup.

Numerical constants occur in most mathematical models, and therefore, it is important to any symbolic regression tools to effectively integrate floating point constants in their optimization search. GP⁷ handles numerical constants by introducing random numerical constants in a specified range to its parse trees. The random constants are moved around the parse trees using the crossover operator. GEP handles the creation of random numerical constants (RNCs) by using an extra terminal “?” and a separate domain Dc composed of symbols chosen to represent random numerical constants.⁹ This Dc specific domain starts from the end of the tail of the gene.

For each gene, RNCs are generated during the creation of a random initial population and kept in an array. To maintain the genetic variations in the pool of RNCs, additional genetic operators are introduced to take effect on Dc specific regions. Hence, in addition to the usual genetic operators such as mutation, inversion, transposition, and recombination, the GEP-RNC algorithm has Dc specific inversion, transposition, and random constant mutation operators. Hence, with these modifications to the algorithm, an appropriate diversity of random constants can be generated and evolved through operations of genetic operators. The values for each genetic operator selected for this study are listed in Table I. These values are selected from various examples given by Ferreira⁹ combined with the trial and error approach. Additionally, to simplify our study, we use the same parameters for all the test cases even though they may not be the best values for the test case under investigation.

Once decent values of genetic operators that can explore the search space are selected, the size of the head length, population, and the number of genes form the most important hyperparameters for GEP. Generally, larger head length and a greater number of genes are selected for identifying complex expressions. Larger population size helps in a diverse set of initial candidates, which may help GEP in finding the best chromosome in a lower number of generations. However, computational overhead increases with an increase in the size of the population. Furthermore, the best chromosome can be identified in fewer generations with the right selection of the linking function between genes. The GEP algorithm inherently performs poorly in predicting the numerical constants that are ubiquitous in physical laws. Hence, the GEP-RNC algorithm is used where a range of random constants are predefined to help GEP to find numerical constants. This also becomes important in GEP for identifying the underlying expression in fewer generations. Finally, we note that due to the heuristic nature of evolutionary algorithms, any other combinations of hyperparameters might work perfectly in identifying the symbolic expressions. In this study, we use geppy,⁷² an open source library for symbolic regression using GEP, which is built as an extension to the distributed evolutionary algorithms in Python (DEAP) package.⁷³ All codes used in this study are made available on Github (https://github.com/sayin/SR).

Compressive sensing/sparse optimization^5,74 has been exploited for sparse feature selection from a large library of potential candidate features and recovering dynamical systems represented by ODEs and PDEs^39,40,45 in a highly efficient computational manner. In our setup, we use this STRidge⁴⁰ algorithm to recover various hidden physical models from observed data. In continuation with Sec. II where we define feature library Θ(U) and target/output data V(t), this subsection briefly explains the formation of an overdetermined linear system for STRidge optimization to identify various physical models from data.

The Burgers PDE given in Eq. (2) or any other PDE under consideration can be written in the form of linear system representation in terms of Θ(U) and V(t),

where β = [$β1,β2,…,βNβ$] is the coefficient vector of size $RNβ$⁠, where N_β is the number of features (basis functions) in library Θ(U). Note that Θ(U) is an over-complete library (the number of measurements is greater than the number of features) and has rich feature (column) space to represent the dynamics under consideration. Thus, we form an overdetermined linear system in Eq. (10). The goal of STRidge is to find a sparse coefficient vector β that only consists of active features, which best represent the dynamics. The rest of the features are hard thresholded to zero. For example, in the Burgers equation given by Eq. (2), STRidge ideally has to find the coefficient vector β that corresponds to the features uu_x and u_2x, and simultaneously, it should set all other feature coefficients to zero.

The linear system defined in Eq. (10) can be solved for β using the ordinary least squares (OLS) problem. However, OLS minimization tries to form a functional relationship with all the features in Θ(U) resulting in all nonzero values in the coefficient vector β. Thus, solving Eq. (10) using OLS infers a radically complex functional form to represent the underlying PDE and generally results in overfitted models. Regularized least square minimization can be applied to constrain the coefficients and avoid overfitting. Hence, regularized LS optimization is preferred to identify the sparse features (basis functions) along with their coefficient estimation. Typical estimation of a sparse coefficient vector with P nonzero entries in β is shown in Fig. 3. A general sparse regression objective function to approximate the solution of the coefficient vector β is given by

where λ is the regularizing weight and ∥β∥₀ corresponds to the L₀ penalty, which makes the problem np-hard. Hence, to arrive at the convex optimization problem of Eq. (12), L₁ and L₂ penalties are generally used to approximate the solution of the coefficient vector β.

The addition of the L₁ penalty to the LS objective function, which corresponds to a maximum a posteriori estimate (MAP) of a Laplacian prior and is termed the least absolute shrinkage and selection operator (LASSO) in compressive sensing, is defined by

However, the performance of LASSO deteriorates when the feature space is correlated.⁴⁰ The sequential threshold least squares (STLS) algorithm was proposed to identify dynamical systems represented by ODEs.³⁹ In STLS, a hard threshold is performed on least square estimates of regression coefficients and hard thresholding is recursively performed on remaining nonzero coefficients. However, the efficacy of STLS reduces when dealing with the identification of systems containing multiple correlated columns in Θ. Hence, L₂ regularized least squares termed ridge regression,³⁶ which corresponds to the maximum a posteriori estimate using a Gaussian prior, is proposed to handle the identification of PDEs. Ridge regression is defined by

Ridge regression is substituted for ordinary least squares in STLS and the resulting algorithm is termed sequential threshold ridge regression (STRidge).⁴⁰ The STRidge framework⁴⁰ is illustrated in Algorithm 1 for the sake of completeness. Note that, if λ = 0, STRidge becomes the STLS procedure. For more details on updating tolerance (tol) to perform hard thresholding in Algorithm 1, readers are encouraged to refer to the supplementary document of Rudy et al.⁴⁰ 

We use the framework provided by Rudy et al.⁴⁰ in our current study. The hyperparameters in STRidge include the regularization weight λ and tolerance level tol, which are to be tuned to identify appropriate physical models. In the present study, the sensitivity of feature coefficients for various values of λ and the final value of λ where the best model is identified is shown. The following sections deal with various numerical experiments to test the GEP and STRidge frameworks.

Partial differential equations (PDEs) play a prominent role in all branches of science and engineering. They are generally derived from conservation laws, sound physical arguments, and empirical heuristic from an insightful researcher. The recent explosion of machine learning algorithms provides new ways to identify hidden physical laws represented by PDEs using only data. In this section, we demonstrate the identification of various linear and nonlinear canonical PDEs using the GEP and STRidge algorithms from data alone. Analytical solutions of PDEs are used to form the data. Table II summarizes various PDEs along with their analytical solutions u(t, x) and domain discretization. Building a feature library and corresponding response data to identify PDEs is discussed in detail in Sec. II.

We reiterate the methodology for PDE identification in Sec. II. The analytical solution u(t, x) is solved at discrete spatial and temporal locations resulting from the discretization of the space and time domains as given in Table II. The discrete analytical solution is used as input data for calculating higher order spatial and temporal data using the finite difference approximations listed in Eq. (4). Furthermore, the feature library is built using discrete solution u(t, x) and higher order derivative, which is discussed in Sec. II. As GEP is a natural feature extractor, core feature library $Θ̃(U)$ given in Eq. (5) is enough to form input data, i.e., GEP terminal set. Table III shows the function set and terminal set used for equation identification, and Table I lists the hyperparameter values for various genetic operators. However, the extended core feature library Θ(U), which contains a higher degree of interactions of features, is used as input for STRidge as the expressive power of STRidge depends on exhaustive combinations of features in the input library. The temporal derivative of u(t, x) is target or response data V(t) given in Eq. (5) for both GEP and STRidge.

Our first test case is the wave equation, which is a first order linear PDE. The PDE and its analytical solution are listed in Table II. We choose the constant wave speed a = 1.0 for propagation of the solution u(t, x). Figure 4 shows the analytical solution u(t, x) of the wave equation. The GEP hyperparameters used for identification of the wave equation are listed in Table IV. We use a smaller head length and a single gene for simple cases like a linear wave PDE. We note that any other combinations of hyperparameters may identify the underlying PDE. Figure 5 illustrates the identified PDE in the ET form. When the ET form is simplified, we can show that the resulting equation is the correct wave PDE, identified with its wave propagation speed parameter a.

The regularization weight (λ) in STRidge is swept across various values, as shown in Fig. 6. The yellow line in Fig. 6 represents the value of λ at which the best identified PDE is selected. Note that in this simple case, STRidge was able to find the wave equation for almost all the values of λ’s that are selected. Table V shows the wave PDE recovered by both GEP and STRidge.

We use the heat equation, which is a second order linear PDE to test both SR approaches. The PDE and its analytical solution are listed in Table II. The physical parameter α = 1.0 may represent thermal conductivity. Figure 7 displays the analytical solution u(t, x) of the heat equation. Table IV lists the GEP hyperparameters used for the identification of the heat equation. Figure 8 shows the identified PDE in the form of an ET. When the ET form is simplified, we can show that the resulting model is the heat equation identified with its coefficient α.

The regularization weight (λ) in STRidge is swept across various values as shown in Fig. 9. The yellow line in Fig. 9 represents the value of λ selected at which STRidge finds the heat equation accurately. Note that STRidge was able to find the heat equation for low values of the regularization weight λ as shown in Fig. 9. Table VI shows the heat equation recovered by both GEP and STRidge. STRidge was able to find a more accurate coefficient (α) value than GEP. Furthermore, a small constant value is also identified along with the heat equation by GEP.

The Burgers equation is a fundamental nonlinear PDE occurring in various areas such as fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow.^75,76 The interest in the Burgers equation arises due to the nonlinear term uu_x and presents a challenge to both GEP and STRidge in the identification of its PDE using data. The form of the Burgers PDE and its analytical solution⁷⁷ is listed in Table II. The physical parameter ν = 0.01 can be considered as the kinematic viscosity in fluid flows. Figure 10 shows the analytical solution u(t, x) of the Burgers equation. Table IV shows the GEP hyperparameters used for the identification of the Burgers equation. Figure 11 shows the identified PDE in the form of the ET. When the ET form is simplified, we can show that the resulting model is the Burgers equation identified along with the coefficient of the nonlinear term and the kinematic viscosity. GEP uses more generations for identifying the Burgers PDE due to its nonlinear behavior along with the identification of feature interaction term uu_x.

The regularization weight (λ) in STRidge is swept across various values as shown in Fig. 12. The yellow line in Fig. 12 represents the value of λ at which the best identified PDE is selected. Note that the STRidge algorithm was able to find the Burgers equation at multiple values of regularization weights λ. Table VII shows the Burgers PDE recovered by both GEP and STRidge. There is an additional constant coefficient term recovered by GEP. Furthermore, the recovery of the nonlinear term using a limited set of input features shows the usefulness of GEP.

The Burgers PDE with a different analytical solution is used to test the effectiveness of GEP and STRidge as the input data are changed but represented by the same physical law. The analytical solution of the Burgers equation (ii) is listed in Table II. The physical parameter ν = 0.01 is used to generate the data. Figure 13 shows the alternate analytical solution u(t, x) of the Burgers equation. Table IV shows the GEP hyperparameters used for the identification of the Burgers equation (ii). Figure 14 shows the identified PDE in the form of ET. When the ET form is simplified, we can show that the resulting model is the Burgers equation identified along with the coefficient of nonlinear term and kinematic viscosity. With an alternate solution, GEP uses a larger head length, more genes, and a larger population for identifying the same Burgers PDE.

The regularization weight (λ) in STRidge is swept across various values as shown in Fig. 15. The yellow line in Fig. 15 represents the value of λ at which the best identified PDE is selected. Note that STRidge was able to find the Burgers equation at various values of regularization weight λ. Table VIII shows the Burgers PDE recovered by both GEP and STRidge.

Korteweg and de Vries derived the KdV equation to model Russell’s phenomenon of solitons.^78,79 The KdV equation also appears when modeling the behavior of magnetohydrodynamic waves in warm plasmas, acoustic waves in an inharmonic crystal, and ion-acoustic waves.⁸⁰ Many different forms of the KdV equation are available in the literature, but we use the form given in Table II. Figure 16 shows the analytical solution u(t, x) of the KdV equation.⁸¹ It can be seen that this analytical solution refers to two solutions colliding together, which forms a good test case for SR techniques such as GEP and STRidge. Table IX shows the GEP hyperparameters used for the identification of the KdV equation. Due to the higher order nonlinear dynamics represented by a higher order PDE, GEP requires a large head length and genes compared to other test cases in equation discovery. Figure 17 shows the identified PDE in the form of the ET. When the ET form is simplified, we can observe that the resulting model is the KdV equation identified along with its coefficients.

The regularization weight (λ) in STRidge is swept across various values as shown in Fig. 18. The yellow line in Fig. 18 represents the value of λ at which the best identified PDE is selected. Note that STRidge was able to find the KdV equation at various values of the regularization weights (λ). Table X shows the KdV equation recovered by both GEP and STRidge. The physical model identified by STRidge is more accurate to the true PDE than the model identified by GEP.

We consider the Kawahara equation, which is a fifth-order nonlinear PDE⁸² shown in Table II. This equation is sometimes also referred to as a fifth-order KdV equation or singularly perturbed KdV equation. The fifth-order KdV equation is one of the most well known nonlinear evolution equation, which is used in the theory of magnetoacoustic waves in a plasma,⁸² capillary-gravity waves,⁸³ and the theory of shallow water waves.⁸⁴ This test case is intended to test GEP and STRidge for identifying higher order derivatives from observing data. We use an analytical solution,⁸⁵ which is a traveling wave solution given in Table II. This analytical solution also satisfies the linear wave equation and hence both GEP and STRidge may recover a wave PDE (not shown here) as this is the sparsest model represented by observed data (Fig. 19). For simplifying the analysis, we remove the potential basis u_x from the feature library⁴² [Θ(U)] for STRidge and additionally include uu_x basis in core feature library [$Θ̃(U)$] for GEP.

Table IX shows the GEP hyperparameters used for the identification of the Kawahara equation. Due to simplifying the feature library, GEP requires smaller head length and single gene. Figure 20 shows the identified PDE in the form of ET. When the ET form is simplified, we can show that the resulting model is the Kawahara equation identified correctly along with its coefficients. For STRidge, the regularization weight (λ) is swept across various values as shown in Fig. 21. The yellow line in Fig. 21 represents the value of λ at which the best identified PDE is selected. Note that STRidge was able to find the Kawahara equation at various values of regularization weights (λ). Table XI shows the Kawahara equation identified by both GEP and STRidge.

Newell-Whitehead-Segel (NWS) equation is a special case of the Nagumo equation.⁸⁶ Nagumo equation is a nonlinear reaction-diffusion equation that models pulse transmission line simulating a nerve axon,⁸⁷ population genetics,⁸⁸ and circuit theory.⁸⁹ The NWS equation and its analytical solution are shown in Table II. We use a traveling wave solution⁹⁰ that satisfies both wave and NWS equations (Fig. 22). We carry similar changes to the feature library that was applied to discovering the Kawahara equation.

Table IX shows the GEP hyperparameters used for the identification of the NWS equation. However, in contrast to identifying the Kawahara equation with smaller head length and single gene from simplifying the feature library, for NWS case, GEP requires larger head length and more genes for identifying PDE as shown in Table IX. This is due to the identification of nonlinear interaction feature u² that appears in the NWS equation. Figure 23 shows the identified PDE in the form of ET. When the ET form is simplified, we can show that the resulting model is the NWS equation identified along with its coefficients. For STRidge, the regularization weight (λ) is swept across various values, as shown in Fig. 24. The yellow line in Fig. 24 represents the value of λ at which the best identified PDE is selected. Note that STRidge was able to find the NWS equation at various values of regularization weights (λ). Table XII shows the NWS equation identified by both GEP and STRidge.

The Sine-Gordon equation is a nonlinear PDE that appears in propagating fluxions in Josephson junctions,⁹¹ dislocation in crystals,⁹² and nonlinear optics.⁷⁶ The Sine-Gordon equation has a sine term that needs to be identified by GEP and STRidge by observing data (Fig. 25). This test case is straightforward for GEP as the function set includes trigonometric operators that help to identify the equation. However, the application of STRidge is suitable if the feature library is limited to basic interactions and does not contain a basis with trigonometric dependencies. STRidge may recover infinite series approximations if higher degree basic feature interactions are included in the feature library.³⁹ Note that the output or target data for the Sine-Gordon equation consists of a second order temporal derivative of velocity field u(t, x). Hence, V(t) consists of u_2t measurements instead of u_t.

Table IX shows the GEP hyperparameters used for identifying the Sine-Gordon equation. For our analysis, GEP was found to be the best model when the larger population size was used. Figure 26 shows the identified PDE in the form of ET. When the ET form is simplified, we can show that the resulting model is the Sine-Gordon equation identified along with its coefficients. Table XIII shows the equation identified by GEP. This test case demonstrates the usefulness of GEP in identifying models with complex function composition and the limitations of the expressive and predictive power of the feature library in STRidge.

This section deals with constructing a modified differential equation (MDE) for the Burgers equation. We aim at demonstrating both GEP and STRidge techniques as SR tools in the identification of truncation errors resulting from an MDE of the Burgers nonlinear PDE. MDEs provide valuable insights into discretization schemes along with their temporal and spatial truncation errors. Initially, MDE analysis was developed to connect the stability nonlinear difference equations with the form of the truncation errors.⁹³ In continuation, the symbolic form of MDEs were developed and a key insight was proposed that only the first few terms of the MDE dominate the properties of the numerical discretization.⁹⁴ These developments of MDE analysis lead to increasing accuracy by eliminating leading order truncation error terms,⁹⁵ improving stability of schemes by adding artificial viscosity terms,⁹⁶ preserving symmetries,^97,98 and ultimately sparse identification of truncation errors.⁵¹ Therefore, MDE analysis plays a prominent role in implicit large eddy simulations (ILES)⁹⁹ as truncation errors are shown to have inherent turbulence modeling capabilities.¹⁰⁰ Discretization schemes are tuned in the ILES approach to model the subgrid scale tensor using truncation errors. As the construction of MDEs becomes cumbersome and intractable for complex flow configurations, data driven SR tools such as GEP and STRidge can be exploited for the identification of MDEs by observing the data.

For demonstration purposes, we begin by constructing an MDE of the Burgers equation,

and discretizing Eq. (14) using first order schemes (i.e., forward in time and backward in space approximations for the spatial and temporal derivatives, respectively) and a second order accurate central difference approximation for the second order spatial derivatives. The resulting discretized Burgers PDE is shown below,

where the temporal and spatial steps are given by dt and dx, respectively. Within the expressions presented in Eq. (15), the spatial location is denoted using subscript index j and the temporal snapshot using superscript index p.

To derive the modified differential equation (MDE) of the Burgers PDE, we substitute the Taylor approximations for each term,

When we substitute these approximations into Eq. (15), we obtain the Burgers MDE as follows:

where R represents the truncation error terms of the Burgers MDE given as

Furthermore, the temporal derivative in Eq. (18) is substituted with the spatial derivatives resulting in

The truncation error or residual of the discretized equation considering u(t, x) as exact solution to the Burgers PDE is equal to the difference between the numerical scheme [Eq. (15)] and differential equation [Eq. (14)].¹⁰¹ This results in discretized equation with the residual as shown as

We follow the same methodology for constructing the output data and feature library as discussed in Sec. II for the equation discovery. However, the output or target data V(t) are stored with the left hand side of Eq. (20) denoted from now as U_er. The resulting output and core feature library are shown as follows:

The computation of the output data V(t) in Eq. (21) can be obtained using the analytical solution of the Burgers PDE. Furthermore, the derivatives in the core feature library $Θ̃(U)$ are calculated using the finite difference approximations given by Eq. (4). We use both analytical solutions listed in Table II for the Burgers equation (i) and the Burgers equation (ii) to test GEP and STRidge for recovering truncation error terms.

We use the same extended feature library $Θ̃(U)$ as input to STRidge given in Eq. (7), but without the fifth order derivative. However, we add an additional third degree interaction of features to $Θ̃(U)$ to recover the truncation error terms containing third degree nonlinearities. The extra nonlinear features that are added to $Θ̃(U)$ are given as follows:

In contrast, GEP uses the core feature $Θ̃(U)$ as input as it identifies the higher order nonlinear feature interactions automatically. This test case shows the natural feature extraction capability of GEP and the need to modify the feature library to increase the expressive power of STRidge.

The functional and terminal sets used for the identification of the truncation error are listed in Table XIV. First, we test the recovery of truncation errors using the analytical solution of the Burgers equation (i) with the same spatial and temporal domain listed in Table II. However, we set spatial discretization to dx = 0.005 and temporal discretization to dt = 0.005 for storing the analytical solution u(t, x). This test case needs a large population size, bigger head length, more genes, and more iterations as given in Table XV, as the truncation error terms consist of nonlinear combinations of features and the coefficients of error terms that are generally difficult for GEP to identify. Figure 27 shows the ET form of the identified truncation error terms. The regularization weight λ for STRidge is swept across a range of values, as shown in Fig. 28. The vertical yellow line in Fig. 28 is the value of λ where STRidge identifies the best truncation error model. Table XVI shows the recovered error terms by GEP and STRidge along with their coefficients. Both GEP and STRidge perform well in identifying the nonlinear spatial error terms with STRidge predicting the error coefficient better than GEP.

In the second case, we test the recovery of truncation errors using an analytical solution of the Burgers equation (ii) with the same the spatial and the temporal domain listed in Table II. We select spatial discretization dx = 0.005 and the temporal discretization dt = 0.1 for propagating the analytical solution u(t, x). This test case also follows the previous case where a large population size, bigger head length, more genes, and more iterations are needed, as shown in Table XV. Figure 29 shows the ET form of the identified truncation error terms. The regularization weight λ for STRidge is swept across a range of values as shown in Fig. 30. In this test case, the coefficients change rapidly with respect to λ, and the best model is recovered only at the value of λ shown by the vertical yellow line in Fig. 30. Table XVII shows the recovered error terms by GEP and STRidge along with their coefficients. Similar to the previous test case, STRidge predicts the truncation error coefficients better than GEP.

In this section, we demonstrate the identification of hidden physical laws from sparse data mimicking sensor measurements using GEP and STRidge. Furthermore, we demonstrate the usefulness of GEP as a natural feature extractor that is capable of identifying complex functional compositions. However, STRidge in its current form is limited by its expressive power, which depends on its input feature library. Many governing equations of complex systems in the modern world are only partially known or in some cases still awaiting first principle equations. For example, atmospheric radiation models or chemical reaction models might not be fully known in governing equations of environmental systems.^102,103 These unknown models are generally manifested in the right hand side of the known governing equations (i.e., dynamical core) behaving as a source or forcing term. The recent explosion of rapid data gathering using smart sensors¹⁰⁴ has enabled researchers to collect data that represent the true physics of complex systems but their governing equations are only known partially. To this end, SR approaches might be able to recover these unknown physical models when exposed to data representing full physics.

To demonstrate the proof of concept for identification of unknown physics, we formulate a 1D advection-diffusion PDE and a 2D vortex-merger problem. These problems include a source term that represents the hidden physical law. We generate synthetic data that contains true physics and substitute this data set into the known governing equations. This results in an unknown physical model left as a residual that must be recovered by GEP when exposed to a target or output containing the known part of the underlying processes. Furthermore, both GEP and STRidge are tested to recover eddy viscosity kernels for the 2D Kraichnan turbulence problem. These eddy viscosity kernels are manifested as source terms in the LES equations that model unresolved small scales. Additionally, the value of the ad hoc free modeling parameter that controls the dissipation in eddy viscosity models is also recovered using GEP and STRidge.

In the first test case, we consider a 1D nonhomogeneous advection-diffusion PDE, which appears in many areas such as fluid dynamics,¹⁰⁵ heat transfer,¹⁰⁶ and mass transfer.¹⁰⁷ The nonhomogeneous PDE takes the form

where c = $13π$⁠, $α=14$ and S(t, x) is the source term.

We use an analytical solution u(t, x) for solving Eq. (22). The exact solution for this nonhomogeneous PDE is as follows:

where the spatial domain x ∈ [0, 1] and the temporal domain t ∈ [0, 1]. We discretize the space and time domains with n = 501 and m = 1001, respectively. Figure 31 shows the corresponding analytical solution u(t, x).

The source term S(t, x), which satisfies Eq. (22) for the analytical solution provided by Eq. (23), is given as