Year: 2021

Abstract:
Discovering the partial differential equations underlying spatio-temporal datasets from very limited and highly noisy observations is of paramount interest in many scientific fields. However, it remains an open question to know when model discovery algorithms based on sparse regression can actually recover the underlying physical processes. In this work, we show the design matrices used to infer the equations by sparse regression can violate the irrepresentability condition (IRC) of the Lasso, even when derived from analytical PDE solutions (i.e. without additional noise). Sparse regression techniques which can recover the true underlying model under violated IRC conditions are therefore required, leading to the introduction of the randomised adaptive Lasso. We show once the latter is integrated within the deep learning model discovery framework DeepMod1, a wide variety of nonlinear and chaotic canonical PDEs can be recovered: (1) up to O(2) higher noise-to-sample ratios than state-of-the-art algorithms, (2) with a single set of hyperparameters, which paves the road towards truly automated model discovery.

heading: 1 INTRODUCTION
text: Mathematical models are central in modelling complex dynamical processes such as climate change, the spread of an epidemic or in designing aircrafts. To derive such models, conservation laws, physical principles and phenomenological behaviors are key. However, some systems are too complex to model with a purely bottom up approach. In such situations, and when observational data is present, automated model discovery tools are becoming increasingly more useful to derive partial differential equations (PDEs) directly from the data. For example, to model the ocean dynamics in Sanchez-Pi et al. (2020) and for embryo patterning in Maddu et al. (2020). From the mathematical point of view, model discovery of PDEs consists in finding F such that,
ut = F(1, u, ux, uxx, ...), where ut is the temporal derivative of the field u and u, ux, uxx, ... are higher order spatial derivatives. Usually F is identified based on an experiment consisting of n samples of the field u, see Brunton et al. (2016); Rudy et al. (2017); Schaeffer (2017); Raissi et al. (2017). Some recent approaches use symbolic regression to find F , see Maslyaev et al. (2019), but so far the most popular approach to perform model discovery is by linear regression which was first introduced in Rudy et al. (2017) and consists in considering F as a linear combination of some candidate terms,
ut = Θ · ξ, where each column in Θ is a candidate term for the underlying equation, typically a combination of polynomial and spatial derivative functions (e.g. u, ux, uux). In order to obtain a parsimonious PDE, many of the coefficients ξ must be zero, motivating the use of a sparse regression method to infer the equation,
ξ̂ = arg min ξ ‖ut −Θ · ξ‖22 + λ ∑ i ‖ξi‖ρ.
The best subset selection is obtained for ρ = 0, Hastie et al. (2015). However it requires solving a nonconvex and combinatorial optimisation problem. The Lasso, Tibshirani (1996), is special in the sense that ρ = 1 is the smallest value of ρ which leads to a convex constraint region and hence a
1Data, code and results shared on: https://anonymous.4open.science/r/sparsistent_ model_disco-56F8/
convex optimisation problem - for which very efficient solvers exist. Furthermore, there is a large corpus of theoretical work available for the Lasso from which the model discovery community could benefit. It is anecdotaly known in the model discovery community that the Lasso does not perform very well compared to other relaxations of the original problem, see Rudy et al. (2017), Li et al. (2019), Rudy et al. (2019) and Maddu et al. (2019), however it has never been studied why. We trace back this lack of performance to the potential variable selection inconsistency of the Lasso.
Let us focus on variable selection consistency or sparsistency, a property defined as: when the number of samples n −→∞, the estimated vector ξ̂ contains the same nonzero terms as the true vector ξ. By first revisiting a geometric interpretation of the irrepresentability condition (IRC) of the Lasso in the context of model discovery, we provide some theoretical insights on Lasso’s inconsistency but also on the design of the libraries generated for sparse regression based model discovery. In addition, these libraries have to be estimated using methods that introduce some deterministic noise into the sparse regression problem. We show when and why the adaptive Lasso Zou (2006) design matrices might have more chances of satisfying the IRC and introduce the randomised adaptive Lasso to perform variable selection under violated IRC. The latter is tuned by stability selection Meinshausen & Bühlmann (2010), which has been used in the past for model discovery in pure sparse regression based approaches Li et al. (2019); Maddu et al. (2019) but without any variable selection error control.
Furthermore, purely relying on sparse regression such as Rudy et al. (2017), Li et al. (2019) and Maddu et al. (2019) heavily limits to low noise and dense data sets, due to the differentiation method used to build the library (typically numerical differentiation or splines). Purely relying on deep learning, see Raissi et al. (2017), to discover F will result in hardly interpretable equations. In Both et al. (2021a) and Chen et al. (2020), the two problems are tackled by concurrently learning a solution of the PDE using a physics informed neural network and inferring an explicit equation by performing sparse regression on a library built by automatic differentiation. However, deep learning model discovery frameworks typically require manually tuning many hyperparameters which are sensitive to the input data: (1) the original DeepMod Both et al. (2021b) requires to tune a threshold to prune coefficients with small magnitudes, (2) PiDL Chen et al. (2020) introduces a couple of multipliers to parameter the amount of physics informed regularisation and the amount of regularisation in the sparsity estimator. Our results show how once the randomised adaptive Lasso with stability selection is integrated within the deep learning model discovery framework DeepMod, a single set of hyperparameters can be used to recover a wide variety of PDEs.
Contributions
• We show the design matrices used to infer the equations by sparse regression can violate the irrepresentability condition (IRC) of the Lasso, even when derived from analytical PDE solutions, i.e. without additional noise. This implies any sparse regression based model discovery framework needs to deal with highly correlated irrelevant variables with relevant ones, no matter the differentiation method used to compute the library.
• We introduce a randomised adaptive Lasso (rAdaLasso) with stability selection and error control algorithm, to recover the true underlying PDE in the presence of design matrices that are highly correlated and violate the IRC.
• By integrating rAdaLasso within the deep learning model discovery framework DeepMod, we show a wide variety of nonlinear and chaotic canonical PDEs can be recovered: (1) at higher noise-to-sample ratios than state-of-the-art algorithms, (2) with a single set of hyperparameters, paving the road towards truly automated model discovery.

heading: 2 THEORY
text: Sparse regression based model discovery sometimes fails to discover the correct underlying PDE from a data set, even when the model is present in the library and contains little noise. To illustrate this, we present a two-soliton solution2 of the Korteweg-de-Vries (KdV) equation in figure 1(a). Considering a library Θ of p =12 terms, the Lasso fails to select the correct terms of the underlying
2Obtained from an analytical solution, see details in Appendix D.
PDE, see figure 1(b), even in the absence of noise. In this section, we explain why it occurs and introduce the irrepresentable condition to identify the cause.
2.1 ON LASSO’S INCONSISTENCY
The Lasso is known to be variable selection consistent under the irrepresentable condition (IRC), see3 Hastie et al. (2015), which requires the existence of η > 0 such that,
max j∈N ||(ΘTTΘT )−1ΘTTΘF,j ||1 < 1− η, (1)
where ΘT is the subset of the design matrix Θ that contains the true model and ΘF,j a column of the subset of the design matrix that contains the rest of vectors. Conversely, the Lasso will not consistently select the true variables for any design matrix, even if they are present in the library.
Geometric interpretation it can be remarked that (ΘTTΘT ) −1ΘTTΘF,j is the least squares solution of,
ΘF,j = ΘT ξ + , (2)
with Gaussian noise, see figure 2(a). The IRC will be violated if any of the projections of ΘF,j onto the column space of ΘT is larger than 1, meaning they lie outside the unit sphere see figure 2(b). The IRC will be trivially satisfied when ΘF,j is orthogonal to the column space of ΘT , in which case η = 1: meaning irrelevant vectors are orthogonal to the relevant vectors Hastie et al. (2015). Furthermore, it becomes clear from the example on figure 2, that if the relevant vectors are orthogonal to each other then the projections of irrelevant terms onto the space of relevant vectors will be exactly the correlations and the IRC will be less likely to be violated. While it is obvious the IRC will more likely be violated in the case the data is highly correlated: correlations among irrelevant vectors will not lead to a violation of the IRC.
3Many formulations of the IRC can be found in literature such as the seminal work from Zhao & Yu (2006) and Meinshausen & Bühlmann (2010) - we do not pretend to be exhaustive on the matter but rather choose the one that better suits our purpose.
A diagnostics metric ∆ based on the IRC we introduce,
∆(Θ, T ) = max j∈N ||(ΘTTΘT )−1ΘTTΘF,j ||1, (3)
where T are the vectors of the support set S and F on its complementary set Sc. In practice, the true support S is unknown. Its estimation Ŝ can be obtained for example by taking for granted the result from a first variable selection. ∆(Θ, T̂ ) can help us determining if our library Θ is sufficiently well designed: if ∆ < 1 we know a Lasso can distinguish Ŝ from Ŝc, otherwise the variable selection result should be taken with caution. It is worth pointing out, that by no means if ∆(Θ, T̂ ) < 1 we can claim T̂ = T .

heading: 2.2 ON THE ADAPTIVE LASSO’S SPARSISTENCY
text: We have seen in the previous section that the Lasso might not be variable selection consistent when ∆(Θ, T ) > 1, meaning that even if the true model is present in the library it might not be selected. Instead, we propose to use the adaptive Lasso which is known to preserve the consistency, see Zou (2006). In figure 1 we illustrate that applying it, results in the correct underlying KdV equation to be recovered. Let us give insights of why it might perform better on model discovery problems. The adaptive Lasso is a two-step estimation procedure where first an initial estimation of the coefficients is obtained to derive a weight vector ŵi = 1/|ξ̂i|γ . We fix γ = 2 throughout this work. ξ̂i is preferably obtained using a Ridge regression to handle multicollinearity. In the second step, a Lasso is applied on the weighted coefficient, ŵi, penalising the terms with their respective weights, i.e.,
ˆ̃ ξ = arg min
ξ̃
( 1 2n ||∂tu− Θ̃ξ̃||22 + λ p∑ i=1 ||ξ̃i||1 ) , (4)
where Θ̃i = Θi/ŵi and ξ̃i = ŵiξi. Let us consider the impact of the transformation on the design matrix by the adaptive weights. Proposition 1. by assuming all relevant coefficients are larger than the irrelevant ones in magnitude and γ ≥ 1, then as a result of the transformation, the projection of the jthirrelevant vector onto the kth relevant vector will shrink by a factor |ξ̂F,j/ξ̂T,k|γ and,
∆(Θ̃, T ) ≤ ∆(Θ, T ). (5) The design matrix Θ̃ will therefore have more chances to verify the IRC than Θ.
See proof on appendix A. By coming back to the example presented on figure 1, ∆(Θ, T ) = 1.69 while ∆(Θ̃, T ) = 2e−14, which illustrates inequality 5. It is worth insisting the correlation matrices of Θ and Θ̃ are identical. However, the projection of the irrelevant vector ux into the space of relevant vectors (uux, uxx) becomes very small, and small enough in this case for the IRC not to be violated anymore (by replacing Θ with Θ̃). This gives insights why the adaptive Lasso manages to identify the true underlying model when the Lasso would not.

heading: 2.3 WHAT HAPPENS WITH MORE REALISTIC LIBRARIES
text: In practice, the design matrix Θ cannot be observed, nor measured, nor derived from analytical solutions, as we did on figure 1, and has to be estimated. In classic model discovery the library Θ is typically built using splines and/or numerical differentiation. In neural network based model discovery, there will also be a non-random misfit between the output of the neural network and the data. As a result, these methods introduce non-random approximations of the higher order derivatives which can be seen as additional deterministic noise δ on top of the ground truth,
ut = (Θ + δ) · ξ (6)
We see in practice that if the correlations introduced by δ are not too large, the true PDE can still be discovered by choosing the correct amount of regularisation for the randomised adaptive Lasso. To see the effect of the interpolation and differentiation methods, we compute ∆ from 2 libraries based on a numerical solution of the chaotic Kuramoto-Sivashinsky (KS) equation (see Rudy et al. (2017))
with varying noise levels. Both libraries contain with p = 36 potential terms and are derived by polynomial interpolations and numerical differentiation. When no noise is added, ∆(Θ0%, T ) = 1.33 and ∆(Θ̃0%, T ) = 7e−3 meaning the adaptive weights used by the adaptive Lasso help casting better the design matrix. However as soon as 1% noise is added, ∆(Θ1%, T ) = 1.38 and ∆(Θ̃1%, T ) = 1.77: the Lasso nor the adaptive Lasso would be able to select the true model. In this case, the interpolation and differentiation methods introduce enough correlations for any design matrix to violate the IRC. On figure 3, the stability plots show the adaptive Lasso would not be able to select the true model even if present in the library, no matter the amount of regularisation; this motivates the introduction of a supplementary ingredient.

heading: 2.4 RANDOMISED ADAPTIVE LASSO (RADALASSO)
text: To work under violated IRC due to correlations that are due to the underlying physical process itself and/or the interpolation and/or the differentiation method, we introduce a randomised adaptive Lasso. It is inspired by the random Lasso presented in Meinshausen & Bühlmann (2010),
ˆ̃ ξ = arg min
ξ̃
( 1 2n ||∂tu− Θ̃ξ̃||22 + λ p∑ i=1 ||ξ̃i||1 Wi ) (7)
where Wi is randomly selected from a beta distribution, w ∼ β(1, 2) to promote weights close to 0. Such randomisation of the regularisation is equivalent to a random rescaling of each column of the design matrix followed by an adaptive Lasso, it is therefore straightforward to solve4. Empirically, we can see the randomisation breaks the correlations and allows to disentangle the group of relevant from the group of irrelevant variables in a stability selection loop: see the stability plots without and with randomisation on figure 3.
Controlling the selection error The determination of λ can be done using stability selection, see Meinshausen & Bühlmann (2010), which has been used in the past for model discovery Li et al. (2019); Maddu et al. (2019). However, error control has not been used and is in our opinion underutilised: as we will see in the experimental results its hyperparameters can become data insensitive. For the sake of completeness we revisit the details of stability selection on Appendix B and recall
4It is (slightly) more computationally expensive than the adaptive Lasso with stability selection as the randomness prevents from the benefit of warm starts.
main results in the following lines. The set of variables to be selected SΛ ∗
stable is determined using two hyperparameters: πthr is the minimum probability of being selected and an upper bound on the expected number of false positives EVmax,
SΛ ∗ stable = { k such that max Π̂λk ≥ πthr for λ ∈ Λ∗ } (8)
where the regularisation path is restricted by an upper bound,
Λ∗ = { λ ∈ Λ such that, E(V ) ≤ q 2 Λ
(2πthr − 1)p ≤ EVmax
} (9)
where qΛ is the average of selected variables. On figure 3(b), the upper bounds on the selection errors are represented using contours on top of the stability plot. The set of variables to be selected SΛ ∗
stable is given by the terms inside a given contour line while being above the threshold πthr.
To conclude this section, in the context of model discovery, sparse regression is usually performed on highly correlated data, due to the data itself and/or to the differentiation method used to estimate the library, which will tend to violate the IRC. This means that even if the true model is present in the library it might not be selected by the Lasso. As a mitigation, we introduce a randomised adaptive Lasso and show once in a stability selection loop with error control, the underlying true model can still be recovered.

heading: 2.5 DEEPMOD INTEGRATION
text: Neural network based model discovery improves the quality of the library with respect to numerical differentiation based methods, see Both et al. (2021b) . We can therefore expect the deterministic noise δ to be much smaller. To leverage such capability, we implement the randomised adaptive Lasso with stability selection and error control in the deep learning model discovery framework DeepMod5, Both & Kusters (2020). The framework combines a function approximator of u, typically a deep neural network which is trained with the following loss,
L = 1 n ||u− û||22︸ ︷︷ ︸ Lmse + 1 n ||∂tû−Θ(ξ̂ ·M)||22︸ ︷︷ ︸ Lreg
(10)
The first term Lmse learns the data mapping (x, t) → û, while the second term Lreg constrains the function approximator to solutions of the partial differential equation given by ∂tu,Θ and (ξ̂ ·M). The terms to be selected in the PDEs are determined using a mask M derived from the result of the randomised adaptive Lasso with stability selection and error control,
Mi =
{ 1 if ξ̃i ∈ SΛ ∗
stable 0 otherwise (11)
where i ∈ [1, p] is the index of a potential term and SΛ∗stable is determined by equation (8). The coefficients ξ̂ in front of the potential terms are computed using a Ridge regression on the masked library (Θ ·M). During training, if Lmse on the test set does not vary anymore or if it increases, the sparsity estimator is triggered periodically. As a result, the PDE terms are selected iteratively by the dynamic udpate of the mask M during the training. In practice, this promotes the discovery of parsimonious PDEs.

heading: 3 EXPERIMENTS
text: In this section, we first show how the randomised adaptive Lasso compares with state-of-the-art sparsity estimators. Second, once within DeepMod, we compare it to the original DeepMod framework.
5The randomised adaptive Lasso promoted here, uses the Ridge and Lasso implementations from scikitlearn, Pedregosa et al. (2011). DeepMod is implemented in JAX, Bradbury et al. (2018)
Comparing with state-of-the art sparsity estimators In order to get an idea of the performance of the randomised adaptive Lasso with stability selection and error control, we compare it to two pure sparse regression based model discovery approaches: PDE-FIND Rudy et al. (2017) and PDESTRIDE Maddu et al. (2019). While the first is a heuristic, the latter solves a relaxation of the best subset selection (l0 regularisation) using an Iterative Hard Thresholding algorithm. To make sure the comparison is fair, we compare our approach with the ones from literature using the data from the original authors of those approaches. Furthermore, we restrict ourselves to cases where the original authors have tuned their algorithms and present the cases as being hard ones, see table 1. In these cases, ∆(Θ, T ) > 1, meaning they violate the IRC, see table 1. The results from the benchmark are presented in table 2. For case 1, ∆(Θ̃, T ) ≈ 1.77 and for case 2, ∆(Θ̃, T ) ≈ 19 explaining why the adaptive Lasso alone will not work in those cases. The result for case 1 is presented on figure 3. From figure 56, with proper tuning both the randomised adaptive Lasso as well as the Iterative Hard Thresholding (IHT) algorithm can recover the true underlying PDE of case 2. However, the computational cost of the IHT is much higher (×100) than the one of the randomised adaptive Lasso (rAdaLasso), which solves a convex optimisation problem.
Impact of rAdaLasso in DeepMod To quantify the impact of the proposed sparsity estimator within DeepMod we compare DeepMod with rAdaLasso and a baseline (the original DeepMod). The latter leverages a thresholded Lasso with a preset threshold of 0.1 (to cut-off small terms) and λ found by cross validation on 5 folds. We simulate model discoveries for the Burgers, KuramotoSivashinsky (KS) and two additional PDEs that introduce different nonlinearities and derivative orders: Kortweg-de-Vries (KdV), ut = −6uux−uxxx and Newell-Whitehead (NW), ut = 10uxx+ u(1 − u2) − 0.4. A single set of hyperparameters is used in all cases see Appendix C. The results are reported on figure 47. Our approach allows to recover all 4 PDEs without overfitting while the original DeepMod would for all, except for the KdV equation. The stability plot obtained on figure 4(b) for the KS equation can be compared to the one presented on figure 3(b): the combination
6The computational cost reported here is obtained by running the code with both the data and hyperparameters from the authors of the original work.
7In terms of computational cost, an epoch takes in average around 0.04s (with 2k samples) on a GeForce RTX 2070 GPU from NVIDIA: discovering the KS equation takes around 90k epochs and around 1 hour.
of rAdaLasso and DeepMod allow to recover the chaotic equation with greater confidence as the probability of selecting irrelevant terms is null.
Benchmarking noise-to-sample ratios We compare the ratios Ψ = α/n where α is the Gaussian white noise expressed in percent of successful model discoveries of additional frameworks for the two most investigated cases in literature: Burgers and KS equations. The compared frameworks are PiDL(deep learning based thats uses the sparsity estimator of PDE-FIND, Chen et al. (2020)), S3d (sparse bayesian learning, Yuan et al. (2019)), SNAPE (basis function approximations based, Bhowmick & Nagarajaiah (2021)) and R-DLGA (symbolic regression based, Xu & Zhang (2021)). The results can be found on table 3. Symbolic regression based approaches such as R-DLGA are more general than sparse regression based approaches in the sense they do need a predefined library, but do not perform as well. For the Burgers equation, the approach proposed in this work can perform the discovery at the same noise-to-sample ratio than PiDL but with a larger library (×2) and a much higher noise-to-sample ratio (×100) for the KS equation. Some limitations of our approach are presented on Appendix F.

heading: 4 CONCLUSION
text: In this paper, we show that no matter the method used to compute the derivatives (numerical or automatic differentiation), the design matrices used for PDE discovery can violate the irrepresentability conditions (IRC) of the Lasso. This means irrelevant variables might be highly correlated with relevant ones. To perform model variable selection under violated IRC, we introduce a randomised adaptive Lasso (rAdaLasso). In addition, it allows to preserve a convex optimization problem and experimental results show it can select the true model in challenging cases at much lower computational cost than state-of-the-art approaches. Furthermore, once integrated in the deep learning model discovery framework DeepMod, the rAdaLasso allows to recover a wide variety of nonlinear and chaotic canonical PDEs up to O(2) higher noise-to-sample ratios than state-of-the-art algorithms. Finally, the hyperparameters used to perform the model discoveries are identical across experiments. These contributions pave the road towards truly automated model discovery. Future work will focus on (1) performing discoveries by leveraging the data from several datasets using multitask learning and (2) including an evolutionary approach in the proposed framework to build the library Θ automatically.

heading: ACKNOWLEDGMENTS
text: Acknowledgements will be included in the final version.
REPRODUCIBILITY STATEMENT
As reported in the abstract of the paper, we share the data and code on a public repository to reproduce our results. In addition,
• for figure 1, Appendix D.1 recalls the expression of the analytical solution.
• for figure 3, the repository allows to reproduce it.
• for figure 4, Appendix D.3 gives more details about the input data for DeepMod. The hyperparameters are identical across experiment and shared in Appendix C. Our repository also contains files with the requirements to reproduce the experiments using python.
The notebooks shared on the repository contain supplementary figures where one can see how well or not the modified version of DeepMod interpolates the data (especially for the chaotic KuramotoSivashinsky equation).

heading: A PROOF OF PROPOSITION 1
text: Proof. by definition for i ∈ [1, p], ŵi = 1/|ξ̂i|γ . We denote Θ̃T,k a vector of the subset of the design matrix Θ̃ that contains the relevant vectors and Θ̃F,j a vector of the subset of the design matrix that contains the irrelevant vectors, then, Θ̃T,k = ΘT,k/ŵk = |ξ̂T,k|γ · ΘT,k and Θ̃F,j = ΘF,j/ŵj = |ξ̂F,j |γ · ΘF,j where k is some column index of Θ̃T and j of Θ̃F . The projection of an irrelevant vector onto the space of relevant vectors in the least-squares sense ΘF,j = ΘT ξ + can be decomposed into the space of relevant vectors, ΘF,j = ∑ k ξF,j,kΘT,k + , where ξF,j,k are the projections of a given irrelevant vector onto the relevant space. Using such decomposition the irrepresentability condition becomes, maxj ∑ k |ξF,j,k| < 1 − η. By applying the transformation
due to the adaptive weights we get, Θ̃F,j = |ξ̂F,j |γ ∑ k ξF,j,k · 1|ξ̂T,k|γ · Θ̃T,k + · |ξ̂F,j | γ , which
can be reduced to Θ̃F,j = ∑ k ∣∣∣ ξ̂F,j ξ̂T,k ∣∣∣γ ξF,j,k · Θ̃T,k + · |ξ̂F,j |γ . Now if we assume that since ξ̂T,k is a relevant term and ξ̂F,j is an irrelevant one: |ξ̂T,k| > |ξ̂F,j |. If in addition, γ ≥ 1, then, 0 ≤ |ξ̂F,j/ξ̂T,k|γ < 1, which leads to,
∑ k ∣∣∣∣∣ ξ̂F,jξ̂T,k ∣∣∣∣∣ γ |ξF,j,k| ≤ ∑ k |ξF,j,k| (12)
which is equivalent to the inequality 5.

heading: B STABILITY SELECTION WITH ERROR CONTROL
text: With stability selection, variables are chosen according to their probabilities of being selected with a warranty on the selection error, Meinshausen & Bühlmann (2010). The first step consists in finding the probability of a variable k of being selected under a data perturbation: let Ib be one of B random
sub-samples of half the size of the training data drawn without replacement. For a given λ an estimation of the probability of k being selected is given by,
Π̂λk = 1
B B∑ b=1 { 1 if fλk (Ib) > 0 0 otherwise (13)
where fλk (Ib) = ||ξ̂λk (Ib)||1 for the randomised adaptive Lasso. Second, by computing the probabilities of being selected over a given range of λ’s the stability paths can be obtained, see figure 3. That range is denoted Λ = [ λmax;λmax], where is the path length and λmax is the regularisation parameter where all coefficients ξ̂ are null. In Meinshausen & Bühlmann (2010) derive an upper bound on the expected number of false positives E(V ), that can help determining a smaller Λ region where a control on the selection error can be warrantied. By fixing this bound to EVmax, the regularisation region becomes, see figure 3,
Λ∗ = { λ ∈ Λ such that, E(V ) ≤ q 2 Λ
(2πthr − 1)p ≤ EVmax
} (14)
where πthr is the minimum probability threshold to be selected and qΛ is the average of selected variables. We propose here to approximate qΛ by q̂Λ = ∑ b |Sb| B = ∑ k Π̂ λ k . Finally, the set of stable variables with an upper bound on the expected number of false positives is,
SΛ ∗ stable = { k such that max Π̂λk ≥ πthr for λ ∈ Λ∗ } (15)

heading: C A SINGLE SET OF HYPERPARAMETERS
text: Stability selection expected number of false positives upper bound EVmax = 3, number of resamples B = 40 and the minimum probability to be selected πthr = 0.9.
Library consists of polynomials and partial derivatives up to the fifth order leading to a library size of p =36 potential terms: {1, ux, uxx, ..., uxxxxx, u, uux, ..., uuxxxxx, ..., u5, u5ux, ..., u5uxxxxx}.
Neural network architecture & optimiser NNs are 4 layers deep with 65 neurons per layer and sinus activation functions with a specific initialisation strategy, see Sitzmann et al. (2020). The NNs are trained by an Adam optimiser with a learning rate of 5 · 10−5 and β = (0.99, 0.99).
Randomness Seeds are identical across datasets: (1) to initialise the NNs and (2) generate the noise vectors.

heading: D ABOUT THE DATA
text: This appendix provides data source details to reproduce the examples of this paper.
D.1 ANALYTICAL LIBRARY OF THE KDV EQUATION
The Kortweg-de-Vries (KdV) ut = −6uux−uxxx PDE analytical solution for 2 travelling solitons is,
u(x, t) = 2(c1−c2)· c1 cosh(
√ c2ξ2/2) 2 + c2 sinh( √ c1ξ1/2)
2( ( √ c1 − √ c2) cosh[( √ c1ξ1 + √ c2ξ2)/2] + ( √ c1 + √ c2) cosh[( √ c1ξ1 − √ c2ξ2)/2] )2 where c1 > c2 > 0, ξ1 = x − c1t and ξ2 = x − c2t. 40 points equally distributed such that x ∈ [−5, 12], 50 points equally distributed such that t ∈ [−1, 2] and c1 = 5, c2 = 2. By automatic differentiation we obtain the 12 terms library: {1, ux, uxx, uxxx, u, uux, uuxx, uuxxx, u2, u2ux, u2uxx, u2uxxx}.
D.2 LIBRARIES FROM SPLINES/NUMERICAL DIFFERENTIATION
Burgers, ut = νuxx − uux, shared on the github repository mentionned in Maddu et al. (2019). The solution here is very similar to the one obtained using the analytical expression below using Dirac delta initial conditions.
Kuramoto-Sivashinky (KS), ut = −uux − uxx − uxxxx, shared on the github repository mentionned in Rudy et al. (2017).
D.3 INPUT DATA FOR DEEP LEARNING EXPERIMENTS
We generate numerical solutions from several equations, on top of which we add α Gaussian white noise,
unoisy = u+ α · σ(u) · Z (16) where Z ∼ N(0, 1). The following PDEs are considered:
Burgers, initial condition: Dirac delta, analytical solution,
u(x, t) =
√ ν
πt · (e
A 2ν − 1)e−x
2 4νt
1 + 12 (e A 2ν − 1)erfc( x√ 4νt )
where A is a constant and ν is the viscosity, ν = 0.1, A = 1 and 40 points equally distributed such that x ∈ [−2, 3], 50 points equally distributed such that t ∈ [0.5, 5].
Kortweg-de-Vries (KdV), see subsection D.1.
Newell-Whitehead (NW), ut = 10uxx + u(1 − u2) − 0.4, numerical solution using a finite differences solver and the following initial condition:
u(x, 0) = 3∑ i=1 αi sin(βiπx)
where αi and βi are constants. 40 points equally distributed such that x ∈ [0, 39], 50 points equally distributed such that t ∈ [0, 1.96] and α1 = 0.2, α2 = 0.8, α3 = 0.4, β1 = 12, β2 = 5, β3 = 10.
Kuramoto-Sivashinky (KS), see subsection D.2. 2000 samples are randomly drawn from a subset of the dataset, details can be found on our github repository, see note 1.

heading: E ADDITIONAL RESULTS
text: Stability plots for case 2 comparison In this case the performance of PDE-STRIDE and rAdaLasso are compared on figure 5.
DeepMod interpolations for the experiments see figure 6.

heading: F LIMITATIONS OF THE APPROACH
text: Incomplete library An obvious limitation of the approach is that the library of potential PDE terms must contain the true underlying PDE terms. Typically this can be diagnosed as we report (in TensorboardX) the mean square error on a test set which allows to verify if the discovered PDE generalizes well or not.
Non-unique solutions If the discovery problem has non unique solutions, our approach will propose one that might not be the one we were looking for. So far, we have devised two cases in which the solutions are not unique and we present them in the next paragraphs. It is known that the solution of a single soliton from the KdV equation is also a solution of the simpler travelling wave equation: ut = −cux, see Rudy et al. (2017). So trying to discover the underlying equation from an analytical of KdV with a single soliton like,
u(x, t) = c
2 cosh( √ cx−ct2 )
(17)
where c > 0, will result in the discovery of ut = −cux. We obtained a similar result with our approach while trying to find the underlying equation from data generated by a solution of Fishers equation (ut = uxx + u(1− u)) from Ma & Fuchssteiner (1996),
u(x, t) = ( 1 + ( √ 2− 1)e−σ(x+2λt)2 )−2
(18)
where σ = λ− √ λ2 − 1 and λ = 5
2 √ 6 . For both analytical solutions x and t play a symmetric role
explaining why ut = k · ux, where k is some constant.
Coefficient bias in chaotic systems for the successful discoveries of non chaotic PDEs, the coefficient mean errors are typically below (< 2%). For the chaotic Kuramoto-Sivashinsky PDE, we show that with 50% noise and a fraction of the data we can recover the terms of the PDE with half the mean error on the coefficients magnitudes with respect to PDE-FIND. However, chaotic systems are sensitive to their initial conditions and the coefficients errors (20%) of the discovered PDE are very large with respect to the ground truth. As a result, the discovered PDE cannot be used for predictions but can be used to identify the underlying physical processes.

Review summaries:
The paper discusses the unidentifiability of lasso-based differential equation discovery, and proposes to use a more stable lasso variant (randomized adaptive lasso), which can more accurately identify PDE coefficients. The method is well-motivated and adresses a major bottleneck in differential discovery. The method is correctly presented, and the initial results are excellent. 

This paper describes an improvement to previous methods for discovering PDEs from data. They show that even for clean data, the design matrices for discovering PDEs using sparse regression can violate irrepresentability condition (IRC), which motivates their approach. They introduce randomised adaptive Lasso (rAdaLasso) and incorporate stability selection and error control. They demonstrate this approach for discovering PDEs from data that was numerically differentiated, comparing to two previous methods on two tricky PDE test cases. They also incorporate rAdaLasso with stability selection and error control into DeepMod (an existing model discovery framework that uses deep learning and avoids numerical differentiation) and show with four PDEs that this is an improvement over the original DeepMod method. Then they compare their version of DeepMod to six previous papers on two PDEs, showing the ability to tolerate more noise and larger libraries. 

This paper deals with the problem of model discovery based on sparse regression. Authors propose a new randomised adaptive lasso for that selects the terms to be used in the PDEs selection that is used in conjunction with Deep learning framework. Authors then experimentally verify on 4 PDEs that their method perform better for model selection and at least as well when looking at MSE loss.   For a detailed summary, authors first verify the violation of IRC which is defined as in eq (1) for Lasso on an analytic PDE (Korteweg-de-Vries (KdV) equation). As their first contribution, authors develop a diagnostic metric that indicates the violation of Lasso when one know the true support of the design matrix. In proposition 1, authors show that under certain assumptions on weights, the adaptive Lasso, where one reweighs the design matrix has better diagnostic score than Lasso.  Since the design matrix is typically approximated, if we incur a deterministic error term, authors empirically show that adaptive Lasso is also susceptible to the violation of the IRC. In lieu of this experiment, authors propose randomised adaptive Lasso, where each column of the design matrix of adaptive Lasso is randomly scaled with $\beta(1,2)$. Author integrate it in deep learning model discovery framework and empirically verify the performance of of rAdaLasso for 1) MSE loss and 2) model selection error.   

This manuscript studies the variable selection consistency in model discovery of PDEs. A randomised adaptive Lasso (rAdaLasso) with stability selection and error control algorithm is proposed to recover the true underlying PDE in the presence of design matrices that are highly correlated and violate the IRC. Better model recovery performance is obtained by integrating rAdaLasso within the deep learning model discovery framework DeepMod.