Improved identification accuracy in equation learning via comprehensive $\boldsymbol{R^2}$-elimination and Bayesian model selection
Authors: Daniel Nickelsen
Bubacarr Bah
Venue: Accepted by TMLR ()

Abstract: In the field of equation learning, exhaustively considering all possible combinations derived from a basis function dictionary is infeasible. Sparse regression and greedy algorithms have emerged as popular approaches to tackle this challenge. However, the presence of strong collinearities poses difficulties for sparse regression techniques, and greedy steps may inadvertently exclude important components of the true equation, leading to reduced identification accuracy. In this article, we present a novel algorithm that strikes a balance between comprehensiveness and efficiency in equation learning. Inspired by stepwise regression, our approach combines the coefficient of determination, $R^2$, and the Bayesian model evidence, $p(y|\mathcal{M})$, in a novel way. Through three extensive numerical experiments involving random polynomials and dynamical systems, we compare our method against two standard approaches, four state-of-the-art methods, and bidirectional stepwise regression incorporating $p(y|\mathcal{M})$. The results demonstrate that our less greedy algorithm surpasses all other methods in terms of identification accuracy. Furthermore, we discover a heuristic approach to mitigate the overfitting penalty associated with $R^2$ and propose an equation learning procedure solely based on $R^2$, which achieves high rates of exact equation recovery.

Published in Transactions on Machine Learning Research (11/2023)
Improved identification accuracy in equation learning via comprehensive R'-elimination and Bayesian model selection
danieln Gaims Computational Statistics and Data Analysis (CSDA), Institute for Mathematics, University of Augsburg, Germany Africon Institute for Mathematical Sciences (AIMS), Cape Town, South Africa
Daniel Nickelsen
Bubacarr Bah
bubacarr. bohl Wlshtm. ac. uk
Medical Research Council Unit The Combia, London School of Hygiene & Tropical Medicine, The Gambia African Institute for Mathematical Sciences (AIMS), Cope Town, South Africa
Reviewed on OpenReview: https: // openrevvee.......... net/ pet//orrmm id= Ock
Abstract
In the field of equation learning, exhaustively considering all possible equations derived from a basis function dictionary is infeasible. Sparse regression and greedy algorithms have emerged as popular approaches to tackle this challenge. However, the presence of multicollinearity poses difficulties for sparse regression techniques, and greedy steps may inadvertently exclude terms of the truc equation, leading reduced identification accuracy. In this article, we present an approach that strikes balance between comprehensiveness and efficiency in equation learning. Inspired by stepwise regression, our approach combines the coefficient of determination, R², and the Bayesian model evidence, P(y in novel way, Our procedure is characterized by an comprehensive search with just a minor reduction of the model space at each iteration step. With two flavors of our approach and the adoption of P(yl.M) for bi-directional stepwise regression, we present a total of three new avenues for equation learning. Through three extensive numerical experiments involving random polynomials and dynamical systems, we compare our approach against four state-of-the-art methods and two standard approaches. The results demonstrate that our comprehensive search approach surpasses all other methods in terms of identification accuracy, In particular, the second flavor of our approach establishes an efficient overfitting penalty solely based on R², which achieves highest rates of exact equation recovery.
1 Introduction
Uncovering the underlying laws governing system is essential for understanding its behavior, and mathemat- ical equations serve as concise representation ion of these laws. Equipped with such equations, predictions can be made, and valuable insights can be derived analytically. The pursuit of inferring these governing equations directly from observations has a long history, dating back to Johannes Kepler's deduction of planetary motion laws in 1609. In modern times, significant progress has been made automated equation inference using machine learning techniques, a discipline commonly referred to as equation learning or symbolic regression.
In contrast to deep learning, equation learning focuses on maximizing model expressivity while minimizing complexity to ensure interpretability. This delicate balance between interpretability and expressivity consti- tutes the central challenge of equation learning. Other challenges include stable feature selection, solvability of learnt equations, computational feasibility, and small data. Greedy algorithms and regularization techniques applied to regression models built. from basis functions are often employed to address these challenges.
Greody algorithms, like stepwise regression and itcrative thresholding, build or reduce the model space by iterative inclusion or exclusion of terms without reconsideration at later steps. Although these methods are computationally efficient, they come at the price of potentially excluding the true model from consideration.

---

Published in Transactions on Machine Learning Research (11/2023)
On the other hand, regularization transforms the search into optimization problem, where all candidate models serve as points in the objective function landscape. However, finding the one global minimum that corresponds to the true model is not guaranteed, and local minima may only accidentally lead to the true model. Multicollinearity among basis functions further complicates optimization algorithms due to the jagged nature of the optimization space.
In this work, we tackle the aforementioned challenges of equation learning by enhancing the stepwise regression approach by stage wise comprehensive model search with only a minor reduction of model space. Our method begins with an elimination process that employs the computationally inexpensive coefficient of determination, R², in order to assess almost all individual candidate models belonging to complexity class. Selecting the best model from that process enables subsequent model selection based on the Bayesian model evidence. The model evidence reflects the probability of a model being the true model for the given data, making an idcal criterion for penalizing overfitting and addressing the central challenge of equation learning: achieving the optimal balance between interpretability and expressivity. By employing specifically tuned conjugate priors within an empirical Bayes framework, we circumvent the computationally demanding estimation of the evidence and instead derive it analytically.
To evaluate our approach, we employ artificially generated data from known models. This choice of evaluation allows us to assess whether our strategy can uncover the ground truth model, serving as a testament to its ability to strike the delicate balance between expressivity and interpretability in realistic scenarios. We compare our method against two standard techniques, least absolute shrinkage and selection operator (LASSO) and least-angle regression (LARS), as well as four ofte--of-theeart methods available in the PySINDy package: sparse relaxed regression (SR3), forward regression orthogonal least squares (FROLS), sequentially thresholded least squares (STLSQ), and best subset selection via mixed-integer optimized sparse regression (MIOSR). Additionally, we propose bi-directional stepwise regression equipped with the Bayesian model evidence (BSR). We evaluate the identification accuracy of approximately 80 scenarios for each system and assess the forccasting accuracy based on 100 initial values for cach scenario and system, amounting to a total of around 8,000 tests. Our findings remonstrate that our proposed methods outperform other approaches in terms of identifying the correct model and can achieve competitive forecasting accuracy.
Our paper is organized as follows. We begin with short overview of existing method in Sec. 2 and introduce our approach and the methods we compare to in Sec. 3 In Sec. 4 we present our nunerical results, which we discuss in Sec. 5 and conclude in Sec. [ A largely self-contained description of our approach with more details and background is provided in the appendix. Our code is available in the Supplemental Material.
2 Related work
The literature on equation learning can be roughly divided into three approaches: evolutionary algorithus,, tree and neural network representations, and regression. A well-known example of evolutionary algorithms is EUREQA, a commercial software used for symbolic regression (Dubčáková 2011 Stoutemyyr 2013). A recent open source alternative to EUREQA also using evolutionary algorithms is the python package PYSR [Crammer] 2023). Tree represent atinss construct equations by combining basic operations and are then employed to minimize regularized objective functions (Vaddireddy et al. 2020), using posterior sampling with sparsily-promooing priors (Jin et al. 2019), or through mixed-integer linear programming (Neumann et al.) 2020 Similarly, neural network architectures have been used to represent equations, where network nodes are replaced by expression building blocks (Martius & Lampert] 2016 Sahoo et al. 2018 Werner al., 2021] ). These neural networks are trained with a regularized minimizer applied objective functions (Rackauckas et al. 2020 Another approach, which allows for the incorporation of physics-infone inffrmee properties such symmetries, has also been developed (Udrescu & Tegmark 2020
These methods find applications in complementing deep neural networks to enhance generalization (Arabshahi et al. 2018) and reducing the data requirement. for training (Yang et. al. 2021 They are also utilized in uncovering complex ecosystem dynamics (Chen et al., 2019).
The above approaches are highly non-linear and typically involve complex optimization algorithms. A simpler approach is based on linear regression models, where features are replaced by basis functions derived from
2

---

Published in Transactions on Machine Learning Research (11/2023)
observed data. Regularized regression ensures sparse weight estimates on the basis functions (Hastie et al. 2009), resulting in concise mathematical expressions. One prominent and widely used method for sparse regression is the LASSO with fr-regularization (Tibshirani) 1996). The advantage of e-regularization is the convexity of the objective functions, which allows for efficient optimization. However, as we will discuss later, while sparse regression performs well with independent features in reconstruction (asks, the correlations introduced by basis function expansion can lead to detrimental instability in equation learning (Su al.) 2017; Rudy et al. 2017). Recent advancements of LASSO, such as SR3, can be found in Zheng et al. 2011 Tibshirani & Friedman (2020).
Greedy algorithms like stepwise regression, FROLS and STLSQ, as well as sparse relaxed regression (i.c. SR3) have proven to be more successful than LASSO in equation learning (Champion et al. 2020 Kaiser et al.)
2018). The latter three algorithms (FROLS, STLSQ, SR3) are implemented in the PYTHON package PySINDy, which is designed for the sparse identification of nonlinear dynamics (SINDy) (Brunton et al.
2016). SINDy,
initially using STLSQ, has seen numerous extensions, including applications to partial differential equations (Rudy et al. 2017], improved noise robust.ness through automated disferentiation (Kaheman al. 2020) and
the weak form of differential equations (Messenger & Bortz 2021a b). supplemented with a-posteriori (MAP)
estimates Niven ct al. 2020), re- re-weighted l1 éi-rogularization (Corticlla ct al. 2021 relaxed rogularization
(Champion et al. 2020 and most recently, best subset selection using MIOSR (Bertsimas & Gurnee 2023
In Nardini et al.
2020), spatio-temporal biological data is first denoised using artificial neural network before
partial differential equations are learnt from the processed data using STLSQ.
In Bayesian linear regression, sparsity-promoting promoting priors are used instead of regularization for equation learning NNyek et al., 2021). Thresholded sparse regression has been translated into Bayesian formulation with implied error bars in Zhang & Lin 20018 2021). By restricting basis functions to quadratic order, the linear structure of the models allows for deterministic results even in the case of the more challenging fo-regularization Schaafffer et al. 2018
3 Equation learning as a regression problem
The goal of equation learning is to find a function f(x) that accurately represents the relationship between the imput variables a and the output variable y. In the context of regression models, we consider a dataset consisting of N observations, where the inputs are organized in a design matrix X and the corresponding outputs are modeled by a random variable Y. We can express the relationship between æ and y as follows:
Y f(x) aZ,
(1)
where the feature vector x is a row of X, Z is a standard normal random variable, and o² is the variance of the noise term. Our goal is to represent the unknown function f(x) using a basis function expansion, iee. as a linear combination of p basis functions kn (x),
P f(z) Σ (x),
where Wn denotes the weights associated with each basis function. Common choices for the basis functions (z) involve products and powers of the individual features Ty. For example, we can have ky (x) = x? k2(xr) = 2122, kz(x)) = 23 and build f(xx) = wa + W2F1122 + 133
Constructing basis function matrix K, with elements Kin kn (æ:), the ordinary least squares (OLS)
estimates for the weights wn are given by Montgomery et al. 2012
y, where y represents a vector consisting of N samples of the target variable Y. With the estimated weights w we can make predictions ÿ = KW.
(3)
Equation learning involves finding a small subset of (x) for which the corresponding weights wn are estimated to fit the data, while the remaining weights can be set to zero without. compromising expressivity.
3

---

Published in Transactions on Machine Learning Research (11/2023)
The challenging part lies indeed in the selection of the appropriate basis functions ka(x), after that, in view of Eq. (3,, the actual learning of the model is straightforward, On one hand, we require the model to have enough flexibility (expressivity) to minimize bias, ensuring that it captures the underlying relationship between x and Y. On the other hand, we want to avoid overfitting, which can lead to high variance in predictions. This trade-off between bias and variance guides the determination of the number non-zero weights, denoted as m, or oquivalently, the model size. addition to finding the appropriate model size, we also aim to identify the "correct" set of basis functions hi (x), which corresponds to recovering the true f()) from data generated by Eq. (1).
3.1 Model class
Apart from the equations that directly be written in the form of Eq. (2) prominent application of equation learning is the learning of dynamical systems,
(tt) = f(x())),
(4)
where (t) denotes the time derivative of (t). To map this problem to Eq. (2). the response variable can be computed from finite differences, i.c. 1/1 #:41-2; for each component of and fixed time step At.
At
It must be noted that the assumption of normally distributed Y in Eq. 1 might not be justified in the context of dynamical system learning due to temporal correlations. On the other hand, if the error structure in (t) is normally distributed, then the differences Ti+1 22 will again be normally distributed, even in the presence of correlations. However, a time varying noise structure would cause heteroscedasticity which could impede statistical methods. In a different approach, the weak form of the differential equation (4) can be utilized; we refer the reader to Messenger Bortz (2021a) for details.
A restriction for the regression models to stay linear in its parameters is that parameters of basis functions may only enter weights w. Basis functions like ea In(a + x), cos as, x², (a+x)m ... with internal parameter a are not suitable.
This restriction might seem quite limiting. On the other hand, the function s(x((), w) defining a dynamical system typically is linear in its parameters w. The reason for that is that these functions often reproduce when diffcrentiated, which can be used to climinate these functions, retrieving the standard linear form in Eq. (2). Many special functions like Bessel, Hankel, Struve and Meijer functions are even defined as solutions of differential equations lincar their coefficients. In general, by considering the differentiated response variable, y + de ~ =((z)))-((z)) if necessary to higher order, we can learn a surprisingly broad class of equations relating y and æ, even relations that do not exist in closed form.
dy
3.2 Regularized regression
The OLS estimates 33 would always involve all terms k, (x) of the basis function dictionary. Equation learning hence is directly related to sparse regression in this context. A common way to mitigate overfitting and thus promote sparsity in regression is to use regularization,
à arg min
||| | -1122 + > Allwolgg]
(5)
E R"
where || 1I.. [r,, 192.199/94 and À is the Lagrange parameter that sets the strength of the e, penalty. The standard, sparsity promoting choice is g == 1 which is known as the LASSO [Tibshirani 1996). The choice
0, for which l/wl|| = Σ, ōw,, is the number of non-zero weight estimates, is often called best subset
selection and requires specialized optimization algorithms (Zhu et al. 2020 Hastie et al. 2020 like MIOSR
Bertsimas & Gurnee 2023
A relaxed regularization like SR3 can be introduced by letting the penalty act on an auxiliary variable u, where distance of to the actual weights controlled by an additional regularization term, e.g. || -uu-a twull
Zheng et al. 2019).
4

---

Published in Transactions on Machine Learning Research (11/2023)
3.3 Model selection
Taking a different roule, sparse solutions of Eq. 3 may be realized by model selection, where basis functions are not discarded by shrinkage of their weights, instead criterion is used to select a model from a set of candidate models. Hore, a candidate model M is defined via a subset of m basis functions kn (1,)....., (x)), which mathematically translates into building a N m submatrix K' from the full N x basis function matrix K. To cope with the huge candidate model space, equation learning utilizing model selection typically involves greedy iterations in which the model space is reduced drastically each iteration.
Selection based equation learning differ in the strategy to trawl through model space and by the choice of selection criteria. A computationally cheap criterion would be the coefficient of determination,
n² y² K' (K'T K y
(6)
yTy
here in a simplified form for standardizod Y (Mongoomrry et al. 2012). However, R² known for its lack of overfitting penally as it would just increase with decreasing sparsily.
Alternatives like adjusted versions of R,2 exist to incorporate an overfitting penalty, but here we make use of the Bayesian model evidence
p(M) x plulM) - do f dw pn pilylue, a, M, Ppr (w,
(7)
where M defines a normal likelihood function (y w. ar M) via Eq. ! substituting K' for K in Eq. ? For a conjugate prior Ppr (w, o), the marginalization above can be done analytically and p(y (M) is known exactly, as detailed in App.
The use of M(yl.M) is often motivated by its excellent overfitting penalizing properties (Occam's razor), which can be ascribed to p(y] MM) being proportional to the probability p(M) of the model M being the true model for the data (y, X) Murphy 2012). The downside of using p(y(M) is the imperative to specify Per (w, o)) even in cases of scarce prior knowledge, which can have a significant impact on p(y]M). In an empirical Bayes approach, we fix Ppe σ) by exploiting the normality property of lincar regression which implics that estimates v are normally distributed with the mean given by the true values of w and the variance given by
T) D)
N p
The details of this procedure included in App.
3.4 Stepwise regression
A standard greedy algorithm to build K' slepwise regression, where terms are added or removed from the model step by step based on criterion (Montgomery al. 2012). Here, we promote using the Bayesian model evidence P(y].M) in Eq. (7) as criterion, which is rarely used due to its intricacics in terms of prior selection and computational cost (Hastie et al.) 2009). Our Bayesian stepwise regression (BSR) procedure is bi-directional, that is, we start with an empty model and add the (x)) that maximizes p(y),M, add a second (x) to K' maximizing p(y]MM, and so on (forward selection). Once p(y)M) cannot be increased further by adding more (z) to K', we start removing columns from K' in the same fashion until again P(y[M) is maximized (backward selection). We continue forward and backward selection until p(y]M) cannot. be increased either selection direction. From including the backward direction and due to the overfitting penalty of M(y] we expect our procedure to be more parsimonious in terms of model size.
Other stepwise regression procedures are FROLS, LARS and STLSQ. In FROLS, forward selection is applied with the correlation to the target variable as criterion such that Eq. (5) with q 2 (Ridge regression) is minimized Billings 2013). A related stepwise algorithm is least-angle regression (LARS) where the correlation between k (x) and residuals is used to build the model in a forward procedure. Following a backward thresholding procodure, STLSQ starts with the full model and alternates between Ridge regression
and removing terms
(z) with weights wn below a pre-defined threshold (Brunton et al. 2016 The STLSQ procedure is the standard method in PySINDY.
5

---

Published in Transactions on Machine Learning Research (11/2023)
3.5 Comprehensive search (CS)
(P) 11 combinations of terms in the model equation which obviously quickly becomes computationally infeasible. However, since for equation learning we are interested in parsimonious models, we may choose small candidate model size m and explore all (m) possible models M defined by K' within that budget. For a fixed model size m, overfitting penalty is not important, and we may use R which is particularly cheap to compute for standardized data, c.f. Eq. (6) In this way, we can single out the best models in terms of R2 for different budgets m.
An exhaustive search algorithm would consider all Σ
To find the best model size m. we utilize the overfitting penalty property of the model evidence p(y]M). Starting with m= we compute M) for fixed number continue to do so for increasing m until a stopping criterion is reached. After that, select the model that maximizes p(y (M) out of the top models selected by R across all considered model sizes. We refer to this method as CS-p
p/2 of best models selected by R², and
Algorithm Generate dictionary of basis functions to produce matrix K
1: Function FUNCDICT(X)
2: Input: data X with N datapoints, features
3: Parameters: maximum degree M1 for individual features, maximum degree M2 for term
4: build all powers "", mj (1,... M1)
pre-computed for speed, limited to power M2
5: initialize counter p p=
6: for all unique I-tuples (ma,...,m do
7:
if Σ, m, E M2 then
D ensure that collective power limited Mr
8:
p: p==++1
9:
P = 11, x
b basis function matrix, candidate models in columns
end if 11: end for
12: Return: basis function matrix K, shape N x
Algorithm Model ranking using R²
Function TopRsq(y, (y. K,m)
2: Input: response data y, basis function matrix K, number of terms for candidate models
3: Parameters: number t= 25 of top models
4: initialize criterion e
b flexible length, to store n² for candidate models D p rows indicating terus part models, flexible number of columns
5: initialize model indices M
6: initialize model number i ==0 0
7: for all n= (m1,..., plp) E 1J,, E,a)=
b possible selections of terms for fixed size m.
8: 1+
9: reduce K' :='==K K.n
b extract column of indicates model with R² valuee o
define candidate model via K'
determine R² for K' and y from Eq. 6
store 6 : R2 and M.
end for
sort columns of and by e (descending)
first columns M now indicate top models in terms of R²
14: Return: top models M.. (shape pxt), criterion c: (length t)
> only return the best models
As a stopping criterion, we may use the first decrease p(y M), but we propose a criterion purely based on R² which proved superior in our study. The R² criterion we propose builds on the observation that true terms (x) have a tendency to be consistently selected in the top R² models. We therefore keep incrementing m until no new kn, (x) is selected consistently, and then build the inferred model from those consistently selected terms. While this is empirically developed method, it is motivated by the reasonable assumption that the maximization of R² is dominated by the true (x)) Therefore, as soon as we look at the R² values of models one term larger than the true model, the procedure is forced to randomly select one extra term in addition to the consistently selected true terms. We refer to this met solely based on R² as CS-R²
Since (P) still becomes very large for larger model sizes m, we add an additional pruning step, in which all (z) that have consistently not been selected for two subsequent iterations are removed from the basis function dictionary for all following iterations. The stopping criterion and the pruning step may classify
6

---

Published in Transactions on Machine Learning Research (11/2023)
Algorithm Comprehensive search with R and M)
1: Function CHS(y, x)
2: Input: data y, X
3: Parameters: maximum number mmax 8 of terms, number a s= of top models, threshold Cmin = 0.75 for term selection
4: K, FUNCDICT(X)
5: initialize feature rating F of shape x mmax
b rate importance of terms for different model sizes mn
6: initialize search := True and number of terms m := 0
7: while search and m S Mmax do
> increment model sixe mn in each iteration until stopping criterion reached
8: m := 111
9:
M, c := TooRssqyy K, m)
obtain list of models and with their It values for fixed model size in append models to index matrix keeping models for all m counts how often terms are selected in top models, weighted by R² criterion b to have ratings between 0 and
10:
append M to Mall
11:
F.m Ej M
12:
normalize For := F.r/ max P..))
13:
r 2 then
14:
index 8o= == (F.r + F..--1 =)
> True if terms not selected for two succesaive model sizes reduce model equation space by those terms
15:
remove K, from K
16:
index in = Frr F
Gmin)
indexes terms with signifieant rating across best models of size
17:
index i2 == r1 Cmin)
b same for provinusly considered model size r-1
18:
if i i2 then
19:
search Falsse
term is selected twice in row like this, conclude search b defines model MM learnt by CS-It" It?
20:
M' == i1, K" K.M
21: end if
22: end
23:
end while
24:
Compute M(y]MM for models Mau using Eq. 7 store (1) maximizing p(aLM)
b defines model M) learnt by CS-p(M)
25:
Obtain OLS estimates vn for K and (1) for K(4) = (MM using Eq. (3)
26: Return: K with v" (model inferred by CS-R*), K()) with wi()) (model inferred by CS-p(y M))
our procedure as a greedy algorithm. However, due to the comprehensive search for each considered m, we consider a drastically larger model space than other greedy algorithms. It is also worth noting that our procodure can be extended to cases where P(y) needs to be estimated by computationally costly evidence estimators, as we effectively reduce the pool of candidate models to just few.
The near-exponential increase of (#) for m < naturally raises the questions of how our model scales with candidate model size m. On the other scales the computation of '22 almost linearly in N and m. Given the fact that in equation learning the goal is to find a minimal model (in fact, models considered in the literature hardly ever excoed m terms), the runtime per itcration in m stays well below 1 min on a standard laptop. A detailed analysis is provided in App. [and Fig. 7
We summarized details of our CS approach in Alg. 3
Illustrations and more details can be found in App.
4 Numerical experiments
We compare all introduced methods (LASSO, LARS, CS-R CS-P(M), BSR, SR3, FROLS, STLSQ, MIOSR) in three mumerical experiment. For LASSO and LARS, we use the PYTHON package scikit-lears (Pedregosa) et al. 2011), for CS-R22 CS-p(M) and BSR we use our own implementation, and for SR3, FROLS, STLSQ, MIOSR we use the PYTHON package PYSINDY (de Silva et al. 2020 Kaptanoglu et al. 2022). As mixed-integer optimizer for MIOSR we used GUROBI with an academic license (Gurobi Optimization ion, L.C., 2023). For all mothods from scikit-learn and PySINDY we used 5-fold cross validation and 3 refinement steps to determine optimal hyperparameters for each application separately. Our own methods, CS-R² and CS-p(M), are not very sensitive to hyperparameters and we worked out the universally best values which, in contrast to the methods we compare to, were used for all applications. The values of the hyperparameters are included in Alg. 3 The BSR method we implemented comes without hyperparameters. Tab. I gives an overview over all considered methods.
In all experiments, the data artificially generated, with the advantage that we know the true model. In the first experiment, we assess the identification accuracy of 100 random polynomials. Since PYSINDY is tailored to learning dynamical systems and cannot directly be applied to learning polynomials, we omitted those
7

---

Published in Transactions on Machine Learning Research (11/2023)
Table 1: Equation learning techniques used in this study building on regression models as in Eqs.
Acronym
Full name
Description
Reference
LASSO
Least absolute shrink- age selection operator
Regularized regression as in Eq. for q=1.
Tibssirani
1996).
LARS
Least augle regression
Forward stepwise regression using correlations between basis functions h2 and residuals.
Efron et al. 2004b.
CS-R2
Comprehensive search R² elimination
For an increasing model size m, the R²² score is calculated for all possible models. Based on the number of times basis functions (x) contribute to the models with largest R², terms (xr) are rated and the lowest rated (x) are excluded. Once no new highly rated k, (x)) are found, the iteration in terminates and the model with largest R is returned.
This work, App.
CS-p(M)
Comprehensive search Bayesian selection
Of models M with highest R² from CS-R² the model M with maximal model evidence p(y(M) returned.
This work, App.
BSR
Bayesian stepwise gression
Bi-directional stepwise selection using the Bayesian model evidence p(y(M) as score.
Hastie et al. 2009
this work for p pyy M),
App. B
SR3
Sparse relaxed regres- sion
The regularization is put an auxiliary vari- able and extra distance term like ||w-ulll is added to Eq. (5)
Zheng al.
2019
FROLS
Forward regression or- thogonal least-squares
Forward stepwise regression selecting terms most correlated with target minimizing Eq. (5) for 2=2.
Billings 2013
STLSQ
Sequentially thresh- olded least squares
Backward selection stepwise regression using results of regularized regression for q = 2 in Eq.
Brunton ct al. 2016
MTOSR
Best subset selection via mixed-integer opti- mized sparse regression
Formulation of regularized regression for q 0 as a mixed-integer lincar program and specialized algorithms.
Bertsimas &
Gurnee 2023)
methods in this experiment. In the second and third experiment the methods are applied to two (chaotic) dynamical systems, where, for the sake readability, we omitted LARS due to its similar performance compared to LASSO.
In all experiments, we have three features polynomials where the power of individual factors was limited to maximum of M1 and the combined power of terms k,,z)) to maximum of M2 = 6. For instance, k(x)) x x2 x would be valid basis function, while k(x) x Ø2 would be excluded for exceeding the individual power limit M1, as would k(Ø) = x1 22 4 for exceeding the combined power limit M2. The resulting basis function dimension of K is = 72.
(x1 22), and we used a basis function expansion of
The results of cach experiment are illustrated statistically by boxplots, where the box indicates the interquartile range, the whiskers extend by factor 1.5 beyond the box, and outcomes exceeding the whiskers are shown individual crosses. The median is shown darker horizontal line, the mean as closed circle.
8

---

Published in Transactions on Machine Learning Research (11/2023)
4.1 Random polynomials
For the artificial data, randomly generated 100 polynomials with 2, 3, and 4 non-zero weights respectively. We restricted the terms of the polynomials to have maximum collective power M2 4, where individual features are restricted to maximum power M1 = 2. The non-zero weights are randomly selected with equal probabilities and their values are uniformly sampled from the set [-4,-1]us [1,4]. We generated artificial data X for each polynomial by sampling from normal distributions with means randomly selected from the interval [-20, 20] and standard deviations such that 5% of their probability mass overlap respoctively. For the polynomials sized 2, 3, and
generated N 20 and N = 65 and N 95 datapoints, respectively. Plugging X into Eq. (1) with f(x;) given by the random polynomial, and corrupting the output with normal noise with standard deviation 0.01, generate data for the response variable Y.
1%
98%
52%
75%
0%
0%
50%
47%
56%
0%
16%
7%
0
Man
Nar
-11
-11
-14
-15
-16
LA550
CS-R!
(S-p(M)
BSS
CS-R2
(3-P(M)
(S-p(M)
Figure 1: Identification accuracy of learning 100 random polynomials measured in terms of the difference Nair between the number of true terms found and wrong terms found. We generated polynomials with 2, 3 and 4 terms (left to right), which also constitutes the highest possible value for Nair and is indicated by a gray horizontal line. The percentages above this line indicate how often all exact terms and no wrong terms have been recovered.
The statistics of the identification accuracy for the 100 polynomials are shown in Fig.
4.2 Lorenz system
As a first dynamical system lest our method, we use the chaotic Lorenz system defined as
=(t) = e(((t) - z(t)),
j(t) = =(p =(t)) y.
2(t) x(t)x(t) -Bz(t).
(9)
The parameters are fixed to its standard values = 10, p 28 and 8 = 8/3. As initial condition,
(in -8, v/ = 8, =0 = 27). We obtain between N = 200 and N = 10000 data points by solving Eq. ( numerically for timestep widths between At = 0.001 and = = 0.1. We corrupt the solutions with normal noise levels between = 0.001 and a = 0.1. Forcing the simulation time to be larger than T == 2, we thus have 80 different scenarios with varying N, At, a and T. We apply the equation learning methods to all scenarios to obtain some statistics on identification accuracy and mean-absolute error (MAE).
We measure the identification accuracy the number of correctly identified equations comprising the model (9). The MAE is obtained by solving the learnt dynamical system equations numerically for randomly selected initial values, and comparing the result with the solution we get solving the true model Eq. D for the same initial condition.
To further increase the sample size for the statistical evaluation the methods, and to exclude the possibility have a particular initial condition that suits a certain method by chance, we sample 100 initial values from the Lorenz attractor and solve each learnt model from all scenarios and equation learning methods for the same 100 initial values.
9

---

Published in Transactions on Machine Learning Research (11/2023)
3.0
3.0
3.0
2.5
(dentfiiatioo accuracy
2.5
identiication
1.5
1.5
ddndificaioon
1.5
1.0
1.0
1.0
0.5
200
500
1500
3500
10000
0.001
0.003
0.01
0.03
24
12
SEAB
10
MIAE
08
8daa
200
3500
10000
0.003
0.01
1000
1000
1000
800
800
Nouke MAL
Novtc:
200
200
200
200
500
1500
3500
10000
0.001
0.003
0.01
0.03
0.2
27
LARO 05-8" CRru5 MRR
SR3 PRITA SITURQ AITOR
JASSO GEAT CRUCAA BER
503
PRONE STEPQ MI0SR
LASSO CENR CRYORE
SR1 FROES ST19Q 3I05R
Figure Statistical results for models learnt from data created by solving the Lorenz system 9 In cach of the nine plots different metric display: the number of equations identified correctly, the MAE to solutions of the true model, the MAE multiplied by the size of the learnt model. The first column uses the statistics across all scenarios, the second column splits it up in terms of number of datapoints N, and the third column in lerms of timestep width AA.
In Fig. 2 we show the statistical results of our experiment. To emphasize the goal to parsimoniously learn models, we also plot the statistics of the MAE multiplied by the size Nmadel of the learnt model. In addition to the overall statistics, we also show the dependency on N and At (splitting the data up in terms of a or T did not reveal any insights).
4.3 Rabinovich- equations
As second dynamical system, use the Rabinovich-Pabrikann equations
((f(()))++((())))))))))))))))))) 7r(t). *(t) (()) -1 1+
-
((t) = (3)(() + 1- - ry(t),
=(1) = ((0())((+(())))))
(10)
We choose the parameter values = 0.14 and 7 = 0.1, and the initial conditions (my = -1.5, 3o = 0, 2g = 1).
We consider 84 scenarios comprising values N 1000 to = = A = 0.001 to A/ = 0.1, and a = 0.0001
and a = 0.01, where we enforce > 5. We solve Eq. 10 numerically and corrupt the solutions with noise as for the Lorenz system.
The results shown in Fig. in the fashion for the Lorenz system.
10

---

Published in Transactions on Machine Learning Research (11/2023)
2.0
2.0
2.0
idtetiiicaioon accuracy
identiiication accuracy
deréification accursyy
1.5
1.5
2.0
1.0
1.0
0.5-
0.5
10000
2000
4000
8000
16000
0.001
0.003
0.01
0.03
0.1
24
0.4
0.6
0.6
0.5
MAB
0.5
0.3
0.4
0.4
0.3-
0.3
1000
2000
4000
8000
16000
0.001
40
50-
35
25
25
20
Novt::
20
15
10:
10
1000
2000
4000
8000
16000
0.001
0.003
0.01
0.03
0.2
27
TAS8O 13-R7
PROTA VTLAQ METER
LARRO CHR DS-AAM
SPA
FROLS
SOOTE
LARO (5-P"
580
FROCS ST1RQ
2000
Figure 3: Statistical results from models learnt from data created by solving the Rabinovich-Fabrikant 10]. The plots have the same structure in Fig. 2 each show the overall statistics split up in terms of N or At.
Fabrikant system shows different metric, whereas the columns
5
Discussion
In this section, we discuss the results of the numerical experiments shown in the previous section.
We begin with the learning of the random polynomials and stress the difficulty of this task, which consists various aspects: i) Randomly generated polynomials are not hand-picked examples where some fine-tuning is possible, ii) The artificially generated data was not tested represent the polynomial uniquely. iii) Even if inferred polynomial deviating from the true polynomial would describe the data well, it does not contribute positively to the identification accuracy.
With these aspects in mind, it is quite remarkable that overall many of the polynomials could be recovered, as shown in Fig. I Particularly striking is the success rate CS-R² solely based on R whereas the LASSO and LARS clearly overfit terms of model size. However, the overall accuracy clearly declines with incrcasing number of non-zero terms. This can be explained with higher chance of terms being selocted spuriously leading
early or too late stopping of exploring the required nunber of terms. Also the risk of erroneously removal of true terms from K increases. We are confident that more refined stopping and removal criteria can overcome these inaccuracies. In combining several criteria, we see an opportunity to improve the accuracy even further.
The Lorenz system has between 2 and terms per equation, and as such was learnt quite successfully, as the first row of Fig. shows. The highost identification accuracy was again achieved by CS-R?, the LASSO was not able to identify any equation, and of the SINDy methods the relaxed regularized regression SR3 and the best subset-selectin using mixed integer optimization MTOSR performed best. In terms of MAF, the
11

---

Published in Transactions on Machine Learning Research (11/2023)
LASSO is among the best, but requires significantly larger models is shown by Mondel MAE. Most of the methods are relatively robust against N and At. Interestingly, Nmodel MAF worsens with more datapoinis in the case of MIOSR but improves for the LASSO, signifying that &o regression "sees" more in the data than there is if given enough data, while & requires more data to produce smaller models (c.f. Eq. 5). Regarding robustness against larger timesteps, it turns out that MIOSR particularly sensitive with deteriorating performance for smaller At, as well as FROLS and the LASSO to a lesser extent.
The Jobinovich-Fabrikann system turned out to be particularly difficult to learn. The identification accuracy for the two CS methods, the BSR. and the relaxed regularized SR3 was on average just below 1 equation, while the other methods failed to identify any equation in all scenarios. No method was able to learn the complete system of equations correctly. Interestingly, the MAE was still reasonable, in particular the LASSO performed well.
A particular problem in learning dynamical systems is that there is no guarantce that the learnt models are solvable. In cases where the numerical solutions failed, we excluded the results from the statistics and kept a record of how often this happened for the various equation learning methods. For the Lorenz system, LASSO, SR3 and STLSQ failed at about 1%, MIOSR at about 20% and the other methods always produced solvable systems. This is little different for the Rabinovich-Fabrieen models bet ween about 1% and 15,, with STLSQ the most problematic followed by SR3 and CS-R2
system, where all methods produeed unsolvable
6 Conclusions
We extensively tested our methods CS-R?, CS-p(MM and the BSR adaption against the standard methods LASSO and LARS, as well as state-of-the-art methods SR3, FROLS, STLSQ and MTOSR. To our knowledge, we are the first to explore equation learning based on a comprehensive model search outperforming existing state-of-the-art methods in terms of identification accuracy, and at least on equal terms regarding forecast. quality.
A disadvantage of our comprehensive search approach is the higher computational cost. However, a direct advanlage is that the model evaluation is trivially parallelizable. Also the little amount of data needed for high success rates is striking - tests on two sets of random polynomials where even done with less datapoints than number of basis functions, N < p. Testing candidate models individually also allows for great ficxibility when it comes to constraints or conditions on models such as solvability, as well as eliminating the risk of getting stuck in local minima of an objective function. It also allows combining several criteria for model and feature selection, in particular complementing existing methods in an independent way with the potential of synergy effects.
The observation that learnt chaotic model comprising all true terms plus few extra terms can decrease the MAE has an interesting implication worth exploring in a future work: It seems possible to learn correction terms from data that lead to a better forecast horizon than the true model itself.
Finally, we contributed towards the utilization of the Bayesian model evidence p(ylM) in equation learning. Here, we benefit from using a conjugate prior for which P(y]MM can be computed analytically and showed in our numerical experiments that our choice of empirical prior is well suited for the tasks considered. However, in general, one would like to have the freedom to select any prior which may entail particularly computationally costly evidence estimators. Performing the comprehensive search with R² can boil down the number of candidate models to a feasible nunber, an approach planned to be explored in a future work.
Considering these possibilities and the promising identification accuracy achieved, we hope to open new avenues of equation learning.
Acknowledgements
DN gratefully acknowledges support by the African Institute for Mathematical Sciences (AIMS) and the Oppenheimer Memorial Trust (OMP). BB has been supported by the BMBF through the German Research Chair at ATMS administered by the Humboldt Foundation.
12

---

Published in Transactions on Machine Learning Research (11/2023)
References
Forough Arabshahi, Sameer Singh, and Anima Anandkumar. Combining symbolic expressions and black-box function evaluations in neural programs. In ICLR, 2018.
Dimitris Bertsimas and Wes Gurnee. Learning sparse nonlinear dynamics via mixed-integer optimization. Nonlinear Dynamics, 1 2023. ISSN 0924-090X. doi: 10.0007/s11077--22--811--- URL https://link wttpinger.om//..100000 wwt.:///0.0..00/111111000000-000
Stephen Billings. Nonlinear system identification: Narmax methods in the time, frequency, and spatio- temporal domains. Nonlinear System Idenlification: NARMAX Methods in the Time, Frequency, and Spatio- Temporal Domains, 08 2013. doi:
Steven L. Brunton, Joshua L.. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113 (15):3932-3937, apr 2016. ISSN 0027-8424. doi: 10.1073/ 20.107/ppnas.le pttp:/pass....1 1517384113 URL Stttp://wwi http://wwwwpppppaaa pnas. org/
tookup//o1111..1033/paaaaa www..//ppppa....m... .ttpppamna...1..1...1.1111 1517384113
Kathleen Champion, Peng Zheng, Aleksandr Y. Aravkin, Steven L. Brunton, and J. Nathan Kutz. A Unified Sparse Optimization Framework to Learn Parsimonious Physics-Informed Models From Data.
IEEE Access, 8: 8:169259 169271, 2020. ISSN 2169-3536. doi: 10.100/ACCEESS.20200 BCCCCCSSSS00.003365533 2020.002362.. URL
https://Leeexppore. iese. org/ htt////tttttnn////////
Yize Chen, Marco Tulio Angulo, and Yang Yu Liu. Revealing Complex Ecological Dynamics via Symbolic Regression. BioEssays, 41(12):1-9, 2019. ISSN 15211878. doi: 10.1002/bies. 201900069.
Alexandre Cortiella, Kwang Chun Park, and Alireza Doostan. Sparse identification of nonlinear dynam- ical systems via reweighted 11-regularized least. squares. Compuler Methods in Applied Mechanics and Engineering, 376:1-38 33, 2021. ISSN 00457825. doi:
cmma20200111622)
Miles Cranmer. Interpretable Machine Learning for Science with PySR and SymbolicReerresiio...] 2023.
URL http://arxxv.orr/aa//33... 01582
Brian de Silva, Kathleen Champion, Markus Quade, Jean- -hrrstsose Loiseau, Kutz, and Steven Brunton. Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49):2104, 2020. doi: 10.21 105//00s..22004.. URL https://ddii.org/10 21101/
joss. 02104
Renáta Dubčáková. Eureqa: software review. Genetic Programmning and Evolvable Machines, 12(2):173 178, jun 2011. ISSN 1389-2576 2576. doi:
URL http:://inn..prrinnee.. com/10.
s 810710-010-9124-2
Bradley Efron, Trevor Hastie, Tain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of Statistics, 32(2):407 499, 2004. doi: 10.1114/000533040000000000 URL https:://bii.rr//100
0090536040000000677
Vincent Fortuin. Priors in Bayesian Deep Learning: A Review. (1):1 28, 2021. URL http://arxiv.org:
abs/2105. 06868
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023. URL https://www.gurobb
Trevor Hastic, Robert Tibshirani, and Jerome Friodman. The Elements of Statistical Learning. Springer
Series Statistics. Springer New York, New York, NY, 2nd edition, 2009. ISBN 978-0-387--888577--000 doi: 10.1007/978-0--887--46555 URT. http://hii..pppieeer..mo//o./... /ttp://///////000/00000000000000
Trevor Hastie, Robert Tibshirani, and Ryan Tibshirani. Best subset, forward stepwise or lasso? analysis and
recommendations based on extensive comparisons. Statistical Science, 35, 11 2020. ISSN 0883-4237. doi: 10.1214/19-STS733. URL Stwerdd--pwwiiiiieieeee Lasso
Sasto-4/Best-Suust-Forrree 1214/19-S79733.
Lasso- Lwsso-analllsssaaaaaeeea la
13

---

Published in Transactions on Machine Learning Research (11/2023)
Ying Jin, Weilin Fu, Jian Kang, Jiadong Guo, and Jian Guo. Bayesian Symbolic Regression. 2019. URL
http::/arxiv.orgg arggabb//9110.06 08892
Kadierdan Kaheman, Steven L. Brunton, and J. Nathan Kutz. Automatic Differentiation to Simultaneously
Identify Nonlinear Dynamics and Extract Noise Probability Distributions from Data. pp. 30, 2020. URL
http: /tappv///axiiv. wwrxiv..rg//aee ab0// 08810
E. Kaiser, J. N. Kutz, and S. L. Brunton. Sparse identification nonlinear dynamics for model predictive control in the low-data limit. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2219), 2018. ISSN 14712946. doi: 10.1098/rssa.20153 htt.sppp...11.33.......
Alan A. Kaptanoglu, Brian M. de Silva, Urban Fasel, Kadierdan Kaheman, Andy J. Goldschmidt, Jared Callaham, Charles B. Delahunt, Zachary G. Nicolaou, Kathleen Champion, Jean- Cean-Christophe Loiseau, J. Nathan Kutz, and Steven L. Brunton. Pysindy: A comprchensive python package for robust sparse system identification. Journal of Open Source Software, 39694)3994,, 2022. doi: 10.21105/joss.03900 URL https: https://doi.org org/ 10. 21805/joss.... 03994
Kevin H. Knuth, Michael Habeck, Nabin K. Malakar, Asim M. Mubeen, and Ben Placek. Bayesian evidence and model selection. Digital Signal Processing: A Review Journal, 47:50 67, 2015. ISSN 10512004. doi:
10.1016/j.ds.2005..06
URL http://áx dww.d.i..ooooo.................. htt..orr/...../////..... dap. 201..06.00000000
Georg Martius and Christoph H. Lampert. Extrapolation and learning equations. 2016. URL http:
//axxiv.com wwa.iiv.rg/abss///.. org/abs//110/ 029955
Daniel A. Messenger and David M. Bortz. Weak SINDy: Galerkin- Data-Driven Model Selection.
Multiscale Model. Simul., 19(3):1474-1497 Jamuary 2021a. ISSN 1540-3459, 1540-3467.
doi: 10.1137/
20M1343166. URL https://epubs. siam. www..oogg///oliie htt.////....111//0011111100 117722011/111
Daniel A. Messenger and David M. Bortz. Weak SINDy for partial differential equations. Journal of
Computational Physics, 443:110525, October 2021b. ISSN 00219991. doi: 20.1016/......22...11 URL
elsevier. httpa//aww........eeeeeeae/////a
Montgomery, E. Peck, G. Vining, and an O'Reilly Media Company Safari. Introduction to Linear Regression Anolysis, 5th Edition. John Wiley & Sons, 2012. ISBN 978-0-470-54281-1. URL https
7/books. google. co. aaabookss7ddDD442EEA.....
K.P. Murphy. Machine Learning: A Probobilistic Perspective. Adaptive Computation and Machine Learn- ing series. MIT Press, 2012. ISBN 9780262018028 URL https://books. google. co. za/books?id=
NZP6AQAAQBAJ
John T. Nardini, John H. Lagergren, Andrea Hawkins- Daarud, Lee Curtin, Bethan Morris, Erica M. Rutter, Kristin R. Swanson, and Kevin B. Flores. Learning Equations from Biological Data with Limited Time Samples. Bulletin of Mathematicnl Biology, 82(9), 2020. ISSN 15229602. doi: 10. 10.1007/s1 10007/11138---1 17388020-00094-2
Nayek, R. Fuentes, K. Worden, and E. J. Cross. On spike-and-slab priors for Bayesian equation discovery of nonlinear dynamical systems via sparse linear regression. Mechanical Systems and Signol Processing, 161:1-22, 2021. ISSN 10961216. doi: 10. 10.1016/j.ymssp.2200 2021.00986.
Pascal Neumann, Liwei Cao, Danilo Russo, Vassilios S. Vassiliadis, and Alexei A. Lapkin. A new formulation for symbolic regression to identify physico-chemical laws from experimental data. Chemical Engineering Journal, 387:123412, may 2020. ISSN 13858947. doi: 10.1016/j..c..20111122/ 12341111 URL https://doi. org/ 10. 10.1016/j.ce) 2019. 2011..123121111....
elsevier httpa///aaaa.aaaaaaaaaaaaaaaaaaa
Robert K. Niven, Ali Mohammad-Djaafar, Laurent Cordier, Markus Abel, and Markus Quade. Dynamical System Identification by Bayesian Inference. In Proceedings of the 22nd Australasian Fluid Mechanics Conference A AFMC2020, dec 2020. doi: 10.14264/692fc55 wi9ffU
URL httpss//espacee library..... .ug........... uq. edu, au/
14

---

Published in Transactions on Machine Learning Research (11/2023)
A. O'Hagan and M.G. Kendall. Advanced Theory of Statistics: Bayesian inference. Volume 2B. Number 2, pt. 2 in KENDALL, MAURICE CFOOGG///FFDATTTSS ADVANCED THEORY OF STATISTICS 6TH ED. Edward Arnold, 1994. ISBN 9780340529225.. URL https://books. co www.z.......otdott DIrEMgEACAAJ.
F. Podregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 2830, 2011.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal Differential Equations for Scientific Machine Learning. 2020. ISSN 2331-8422. URL http://arxiv. http:/axxxvv..rg/bbbbe/e abs///b////// 2001b04485
Samuel H. Rudy, Steven L. Brumton, Joshua L.. Proctor, and J. Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4):1-7, apr 2017. ISSN 2375-2548. doi: 2..116//ciadv.16000
URL atttps://rt http://wwww..wnc.....///o
science...... scwence.org//oo///...
RUser4512.
Computational complexity of machine learning algorithms.
https:/ httpt://wwwtttht www.theeernelttrp.ccomm com baachine/eerrinng//oee httpan///ooootttttoaaoaappttttaa alaatr//nnnngggggttttsss///////
April 2018. Accessod: Nov 2023.
Subham S. Sahoo, Christoph H. Lantpert, and Georg Martius. Learning equations for extrapolation and control. 35th International Conference on Machine Learning, 1CML 2018, 10:7053 7061, 2018.
Hayden Schaeffer, Giang Tran, and Rachel Ward. Extracting Sparse High-Dimensiona Dynamics from Limited Data. SIAM Journal on Applied Mathematics, 78(6) 3295, jan 2018. ISSN 0036-1399. doi: 10.1137/18M116698X, URL htttps://epbbs..... http://pppbs.ssaaaoorrr/ooa///// /ttp////0000011111181111616689
David R. Stoutemyer. Can the Eureqa Symbolic Regression Program, Computer Algebra, and Numerical Analysis Help Each Other? Notices of the American Mathematical Society, 60(06):713, jan 2013. ISSN 0002-9920. doi: 10. 10.1090/noti1000. URL htttp://wwx ams www.ar/......................... http:////ooggg//ooooore/ee/aaaaa
Weijie Su, Malgorzata Bogdan, and Emmanuel Candès.
False discoveries occur early on the
Lasso path. The Annals of Statistics, 45(5):2133 2150, oct 2017. ISSN 0090-5364. doi: 10.
1214/16-A0S1521. URL http:///w//pooottaococcoo/oo///a
Robert.
Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Stalislical Sociely. Series B (Methodological), 58(1):267-228 288, 1996. ISSN 00359246. URL http:://ww..jtttee org/stable/
2346178
Robert Tibshirani and Jerome Friodman. A Pliable Lasso. Journal of Computational and Grophical Statistics, 29(1):215 225, 2020. ISSN 15372715. doi: 10.1080/10618600.0011 10611600.2019.1648222
Silviu Marian Udrescu and Max Tegnark. AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16), 2020. ISSN 23752548. doi: 20.1126/sciadv.aay//
Harsha Vaddireddy, Adil Rasheed, Anne E. Staples, and Omer San. Feature engineering and symbolic regression methods for detccting hidden physics from sparse sensor observation data. Physics of Fluids, (1):015113, jan 2020. ISSN 1070-6631. doi: 11166/1.5113351. URL 10. 10.1063/1.51366
W Von Der Linden, R Preuss, and V Dose. The prior-predicivve value: A paradigm of nasty caati-dimenssonal integrals. In Marimum Entropy and Bayesian Methods Garching, Germany 1998, pp. 319 326. Springer, 1999. URL http:://innnssppingge. http://///.//cppeere///00/0////
Matthias Werner, Andrej Junginger, Philipp Hennig, and Georg Martius. Informed Equation Learning. pp. 1 1-24, may 2021. URL http:/ httt:://ant //trx/vaavvv.///////aa http//aai...rrabs/////aaa
15

---

Published in Transactions on Machine Learning Research (11/2023)
Liu Yang, Xuhui Meng, and George Em Karniadakis. B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data. Journal of Compulational Physics, 425: 109:10 jan 2021. ISSN 00219991. doi: 10. 10.1016/1.cp.22200109999 URL http:://iikkighuu.eeess.
com/
Sheng Zhang and Guang Lin. Robust data-driven discovery of governing physical laws with error bars. Proceedings of the Royal Sociely A: Malhemalioal, Physical und Engineering Sciences, 474(2217):201 80305, doi:
sep 2018. ISSN 1364-5021.
URL httpe://roggaaadietyppulddiiiig hrggiiiiingggg///////g//
101/io..098/rspa.. 2018.03010
Sheng Zhang and Guang Lin. SubTSBR to tackle high noise and outliers for data-driven discovery of differential equations. Journal of Compulational Physics, 428:1 32, 2021. ISSN 10902716. doi: 10.1016/j.
jcp.2020.009962.
Peng Zheng, Travis Askham, Steven L. Brunton, J. Nathan Kutz, and Aleksandr Y. Aravkin. A Unified Framework for Sparse Relaxed Regularized Regression: SR3. IEEE Access, 7:1404-1423, 2019. ISSN 2169- 3536. doi:
2886628. URL
Junxian Zhu, Canhong Wen, Jin Zhu, Heping Zhang, and Xueqin Wang. A polynomial algorithm for best-subsct
selection problem. Proceedings of the National Academy of Sciences, 117:33117 33123, 12 2020. ISSN 0027- 8424. doi: 10.1073 /1100///ms...// pppp:/aaaa2220.0... 204421117 URL httpp://ppaa..orr////////////010
16

---

Published in Transactions on Machine Learning Research (11/2023)
A Linear regression
In this section, we show how linear regression can be applied to equation learning, which sets the basis of this work. Details to regression can be found in Montgomery et al. 2001
Starting point are N observations (iy, where the index denotes data points, and the index j denotes features. We assume that the response (dependent) variable Y is given as a function of explanatory (independent) variable x,
== [(x) oZ,
(11)
where f(x) defines the model, Z is a standard normal random variable, and o² is the variance of the noise term. The explanatory observations are often organized in terms of a design matrix X, with features columns and datapoints in rows, Xij = Zij-
We assume that f(x can be given in terms of p basis functions kn (x)),
f(x) Σ (x),
(12)
with weights Wh. known as basis function expansion. Similar to the design matrix, we can define a basis function design matrix K with elements Kin = kn (x:).
The ordinary least squares (OI.S) estimates for wn are known to be
i arg min RP
Kw 2112 (K" K) -²
(13)
with the g-norm
[w/lg
Σ |wm]19]
1/9
(14)
Predictions based on the OLS estimates are then given by
i = K'
(15)
From the normality of linear regression, it is known that the estimates v follow normal distribution with the mean given by the true values for w and the variance given by 6² (y 2) (1-0) properties later to empirically define the prior in the Bayesian description.
We will use these N
In view of Eq. 13), once choice of k (a:) is made, the actual learning of the model is straight forward. The difficult part is the choice of (x):: the one hand we require sufficient expressivity of the model to minimize bias, on the other hand we want to avoid overfitting to minimize variance of predictions. This bias-variance trade-off essentially dictates the number of kn (x),, i.e. the effective dimension of feature space or complexity. Apart from the appropriate model size, we also seek the "correct" " ki (x)), in the sense that the true ((x) is recovered from data generated by Eq. 11
A common approach to avoid overfitting is regularization,
à arg min E RP
Kw w + X Al|wolgg]
(16)
where || ISn, w.. 19/ and X is the Lagrange parameter that sets the strength of the E, penalty. Common choices for are = (Ridge regression), q q= 1 (the standard, sparsity promoting choice known as the LASSO), or combinations such as clastic net. A special case of regularization is q = 0, for which ||l]|| Σ, 8w..u the number of non-zero weight estimates the regression procedure with this penalty is often called best subset selection and requires specialized optimization algorithms.
17

---

Published in Transactions on Machine Learning Research (11/2023)
A standard measure for goodness of fit is the coefficient of determination, R², which relates the variance explained by the prediction & to the variance of the observations y,
y 9) y 9)
R²
(17)
y yy
where - Di is the empirical mean. Assuming standardized y, R² can be simplified to
R² 1- yy y y T Kw . 'T y + 2 T K K
(18)
(y 70 -w)
1 y" y y Kr
(19)
(y 7) (y 4)
1 y" w K y
(20)
(y () (y- 7)
1 y J y w K y
(21)
y"y y
y² K (K - K)- -1 K y
(22)
y y
where plugged in Eq. 15] in the 1st line, in the 2nd line we used that Kyy =T in the 3rd line that K = 55 3j K 0x = Ei Dr K UR T in the 4th line we assumed ÿ = 0 due to standardization, and in the last line we used to (KTK)-IK' y- For sparse weight estimates u R² is extremely efficient to compute.
A value of R² close to signifies good predictions ÿ. However, it is well known that R can always be brought closer to 1 by increasing the number of features in the model, thus essentially lacking any overfitting penalty.
For the purpose of adequate model selection, it is helpful to formulate regression in a Bayesian setting. The main step to this end is defining the likelihood distribution for Y. In the simple case of Eq. 11 Y is normally distributed,
pulylw, MM = NN,,,,,
(23)
where the mean vector is given by / = fx,). Here, we have included the model M as a condition for the likelihood. In general, the model M can be specified in any way that determines the likelihood apart from its parameters, e.g. the underlying distribution assumption, but here M is equivalent to the choice n of terms being part of the model equation f(x) in Eq. (1).. More specifically, is boolean vector of length p matching an ordered list of the terms of the basis function dictionary, indicating with that the corresponding (erm is part of the model. The model size is therefore given by m ンー= f1j
Maximization of the corresponding log- likeiiood function reproduces the OLS result 13 Encoding existing
information on the weights Wn as a prior distribution ppr(w, σ), Bayes' formula implics for the posterior distribution
Pu(yl w. M) Pp)( P(ylM)
Ppo(w: (ly)
(24)
where the normalization factor P) is known as the evidence or prior-prodictive value. At this point, we have substituted the N-sized data y for the larget variable Y under the assumption factorizing likelihood distribution.
The marginal likelihood M) may be interpreted as the likelihood of the observation y given the model M regardless the choice of parameters. Hence, taking p(y)M) as the likelihood function the level of models, we may use Bayes' formula again to obtain
(M P(y) Ppr (M) p(y)
(25)
18

---

Published in Transactions on Machine Learning Research (11/2023)
where (P)( is prior for set of models. Without further information the model set, we may assume constant Ppr (M), and find that the model evidence P(y]M) is proportional to the probability of a model, Therefore, if we could maximize p(y(M) over all M, we would in fact identify the model M that most likely explains the observations y.
Poo given the observations
From Eq. 24 it follows that the evidence is given by
p(ylM)= = f do f dw pulylu, o, M) Ppr (w,(
(26)
which in general is not solvable analytically, and also poses a particularly tough numerical challenge (Von
Der Linden el al. 1999: Knuth et al.) 2015). Fortunately, by choosing Ppr(w, a) conjugate to pn (ylw, a, M), the integral becomes solvable analytically. The conjugate prior, however, is not necessarily the sensible choice from the inference point of view. In fact, making a good choice for the prior is a much debated problem Fortuin 2021). Here, we demonstrate that for the purpose of linear equation learning, the conjugate prior is suitable choice, if hyper-parameters are distilled from data. The question whether other choices for the prior would perform siguificantly better is left for future research.
The conjugate prior for the likelihood
[3 is the gamma gamma-normal distribution O'Hagan & Kendall 1994
det E T|, E, k, (2 (
7D/2+k-1 e- i (Iw
(27)
with mean vector l and precision matrix Σ for the weights w, and shape k and scale for the precision == 1/02. Plugging Eq. (27) and Eq. 23 into Eq. (26) and performing the integration, we obtain for the
log-evidence per data-point the closed expression
1 Inp(yl/) = 2N In
det Σ
1
k
E 2
)
In 2T
+
In
N
det
A
2
2
N
1
N
In 1 +
In
r
+k
In I(k)
(28)
N
N
2
N
with
A = K"K + +.
(29)
b = KT y Ep,,
(30)
{ y y + Ep 6T A-
(31)
We detail the calculations in App. B
Owing to the normality of linear regression, and from standardizing the data, it is reasonable to assume the following parameters for the prior: For the mean vector, we choose μ = w and the precision matrix Σ is taken to be diagonal with elements diag(E) "..mm"K.j) resulting in normal distributions broadened by a factor N to make the prior more uninformative. The gamma distribution entering Eq. (27) has the mode (k- 1)), which we set to due to standardized y. The scale is set to & 1/2 which appears to be broad enough for an uninformative prior.
I-p/N
For completeness, we mention a few more selection criteria. The adjusted R² Montgomery et al. 2012
N - N m
R2
(1 R²)
(32)
equips the usual R² with an overfitting penalty. The Akaike information criterion (AIC) measures the loss of information by using the inferred model instcad of the (unknown) true model (Murphy 2012), and similarly but derived from the model evidence (26) in the big data limit, follows the Bayesian (Schwarz) information criterion,
AIC 2pm 10 σ)- 2m, BIC Pm (y) δ) 2m In N
(33)
19

---

Published in Transactions on Machine Learning Research (11/2023)
Similarly, but derived from the model evidence 34 in the big data limit, is the Bayesian (Schwarz) information criterion,
BIC = Pu ō) 2m In N.
(34)
12) a prominent application of lincar equation learning is the sparse identification of dynamical systems,
Apart from the equations that can directly be written in the form of Eq.
a(1) ((x((),
(35)
where i(() denotes the time derivative of x(1). To map this problem to Eq. 12 the response variable can be computed from finite differences y/ = شسسس for a fixed time step At. In different approach, the weak form of the differential equations can be utilized; we refer the reader to Messenger & Bortz (2021a) for details.
A restricttion for the regression models to stay linear in ils parameters is that the parameters of basis functions only enter as weights w. Basis functions like
e0 In(a + x), cos ax, x", (a + x)
(36)
with internal parameter not suitable.
This restriction might seem quite limiting, On the other hand, the function s(x((t, w) defining a dynamical system typically is linear in its parameters w. The reason for that is that functions shown in Eq. (36) often reproduce when differentiated, which can be used to eliminate these functions, retrieving the standard form shown Eq. (12 and Eq. [11.. being linear paramcters. Some special functions like Bessel, Hankel, Struve and Meijer functions are even defined as solutions of differential equations.
In general, by considering the differentiated response variable,
dy dr
y(x++1) v(22)
(37)
Fi41
22
if necessary to higher order, can learn surprisingly broad class of equations relating y and æ, even relations that do not exist in closed form. Here, we restrict ourselves to dynamical systems and leave the full exploration of learning in this broad model class for future work.
While many equations with non-linear parameters can rewritten linear form by differentiation as explained above, some functions like x", (a+=)* can only be reduced fractions or require many differentiations which can give rise to numerical issues. Therefore, if we were able to include fractions the basis function expansion, we could expand the model class even further. The problem is that basis functions like 1$ are divergent for certain values of and odd n.
It is, however, possible to modify the regression model to also incorporate fractions. In its simplest form, may consider
En Wn hin (x) + oz Σm= kim (x))
(38)
with different (sparse) weights Un but same set of basis functions for the denominator. For the next step, we assume that (x) is part of the numerator, that is wi + O, and we can rewrite
(x) Um W1
Wn (x).
ykm (x)
(39)
WT
n=
In this form, kj takes the role of the response variable, and we have a second set basis functions given by ykm- For given model of this form, the weights also follow deterministicaly from Eq. (13), only wi neods to be determined from 1-dimensional numerical root-finding algorithm. We leave this ansatz for future research.
20

---

Published in Transactions on Machine Learning Research (11/2023)
A similar idea has been proposed in Kaheman et al. (000)) where tu regularized objective function needs to be minimized for each possible basis function taking the role of lą above. Our comprehensive search strategy naturally includes this procedure as straight forward possibility, which is planned to be investigated in future work.
B
Exact Bayesian model evidence
The model is given by
V4 Σ wn WnKKin + 21/ VT
(40)
where Kin = (x:)) is the basis function design malrix, wn are the weights, r= 1//22 is the precision, and 2/ ~ ~V 1). For the whole vector of N responses, we can use the multivariate normal for the likelihood,
N/22
T (y ply K, w) = (2T) exp (- (y-Kw))
Kw)):
(41)
where the precision matrix is diagonal with identical the diagonal. Since the weight parameters wn, enter quadratically, we can rewrite this expression in normal form for w,
N/2
Kl w.
N/2 exp (--5)
x exp
T
(w - -0)² KT K (w (w)
(42)
2
with the residual sum of squares
S= (y Kw)) (o-KKi)(Yy-K (z)
(43)
= 'リ 11 T KTy y
(44)
TTy y y "KK (TkKK)
(45)
The mixed lerms cancel after plugging in KTy = K T K from the known OLS solution v (K*
The above is of the form of a gamma distribution for r multiplied with normal distribution for ap conditioned on r. If we use prior of the same form, we keep the form for the posterior, and thus have found the conjugate prior.
As a prior for the weights w, we choose
TP/2 /dot E
p(w μ. E)
(2x)p/2 exp (-7
T (w 11) T (w 11)).
(w-H)).
(46)
2
where E the precision matrix with 7 split off, and μ is the mean vector the multivariate normal prior. Splitting off r technical means that specifying is relative the unknown т, but does not need to be known for that, as we integrate over all possible values.
For the posterior,
interested in the quadratic form involving w,
(w- T0)TT KTK K (w-
: + (w ((m-1)) 11) Σ (w -u--10 p
(47)
= T A 2w b
(48)
21

---

Published in Transactions on Machine Learning Research (11/2023)
with
A = K + E,
(49)
b = KJ Ki + Ep
K"y En,
(50)
w K K + 4 En
= W" y Ez Σp
(51)
-y"K(K" k 1 E/..
(52)
Put into this form, we can use the Gaussian integral /d e 30² Aoto"o= (2xr) det A
2200A-1 to marginalize,
ply r) p(y K, ne, E)
(53)
= f dw plyl K, w. () p(w n. E)
(54)
VTN+D det E e #(++)) / dw e Aw-""
(55)
V2
VTN+P det E
V2T)
((5+))
ebb" A-bb
(56)
((2x)
P A
7 det E
E 1-16)
(57)
V (2T) det
For the r-integration, we choose the gamma distribution
p) I(k)o xk- exp (- (TT/D)
(58)
as the (conjugate) prior for T, and define
E = - y 11 E1 A
(59)
The remaining r-integral follows then from z0º dr -r - = - as
p puy K, μc E, k. 1)
(60)
% dr p 0) p(y 7)
(61)
det Σ
00
T( (2T) 7
det A
/
(62)
r +k)
del Σ
{
1 9
2 -K
(63)
r(k)o* (2T)
det A
2
and for the log-evidence per data-poiis taint we obtain
1
det Σ
N
Inp(y)
=
2NV
In
det A
2
In 2m
k
C
1
k
In
In
2
N
2
N
N
+
r
+k
In [
(64)
N
2
N
22

---

Published in Transactions on Machine Learning Research (11/2023)
C Details and illustration of comprehensive search equation learning
To our knowledge; all equation learning approaches that consider the full candidate model space include optimization step in various representation spaces of equations. Here, we propose a stralegy that does without any mumerical optimization algorithms, and instead considers candidate models individually in an almost exhaustive manner. Since already small dietionaries of basis functions can lead to tremendous numbers of candidate models, a combination of cheap selection criteria and successive reduction of model space with a suitable slopping criterion is required. We demonstrate how the simple criterion R2 can be used for such a comprehensive search.
In first step, dictionary of basis functions is generated using Alg. 1 These basis functions consist of all possible products of available features zy. In these products, the factors are raised to all possible combinations of powers (line 9), where we restrict individual powers to M1 (line 4) and the combined power of term to M2 (line 7). For example, for M1 =3 and M2=5, 5, the term at43 would not be allowed since 4 >> M1, and the term xf2223 would not be allowed because 2+1+3>M;
The comprehensive search (CS) strategy we propose is based on this basis function expansion and described by Alg. 3 It begins by considering all regression models with mn non-zero weights wy (line 9), c.f. Eq. 12 To this end, the auxiliary Alg. 2 produces list of models with top R values by looping through all candidate models of size m. These models are returned as an index matrix M, indicating selected terms with 1 and deselected terms with a O, where each column stands for a candidate model (lines 5,11). The models are sorted in descending order with respect to R² (line 13).
Back to Alg. 3 we successively increase m starting from m=1 (lines 7-8). We found that for a fixed m value, R² performs particularly well in identifying the best model out of the millions of models (see for instance the left plot in figure Aj. To infer value for m with just R², we create feature rating matrix F defined as the weighted counts of terms being selected across s top models, where the weight is given by R². Based on F, we chock for terms selected by R² for two successive model sizes m and m-1 (lines 16-21). If for both m and m-1 the same terms are selected consistently, we choose these two terms as part of the inferred model M by CS-R² and conclude the search.
As for larger values of m the number of candidate models can casily reach hundreds of millions, we implement another strategy to divide out the list of candidate models. If terms have not been selected for two successive model sizes and m- we remove these terms from the design matrix K (lines 14-15). In this way, we reduce the model space as we go along.
The rationale for this selection and elimination strategy being solely based on R is the following: If the true model has m terms, and we are testing all models with m = - terms, then the models with largest R² will consistently be composed of the m 2 terms that contribute the most to explaining the variance of y. The other true terms will be selected sporadically but at least once, terms that have not been selected at all can hence be removed from the candidate models. Testing in the next stage all models with m m terms, one more term will be consistently selected. The same holds for testing models with m² terms, but when testing models with +1+ terms, no new term can contribute consistently to explaining more of the variance of y. While the R²² measure will increase for models with m*+1 terms, compared to models with m terms, no extra term will consistently be selected. Therefore, once no new term is selected consistently when incrementing the number m of terms, we may conclude that all contributing terms have been found. An illustration of this strategy can be found in Fig. where a typical case with m= = is shown.
In a final step, the list of top models from the R² evaluation can be combined and each tested with the Bayesian model evidence P(y]M) (c.f. App. B). In Fig. 5 the selective power of p(y)lM is demonstrated. Fig. 0 schematically summarizes the CS approach.
The perpparrametess of our procedure are the number 8 of models used to count consistent selection of basis functions k, (z), together with the threshold Cmin that (normalized) count of a feature must exceed to be selected, the number t of the best models in terms of R² that are combined a new list of top models for re-evaluation with p(y]MM. and the maximum number of iterations, Mmax. We found that the universally best values are p/2, where p is the number of basis functions, Gmin = 0.75, = 25, and mmax = 8.
23

---

Published in Transactions on Machine Learning Research (11/2023)
9.2
0.08
10
0.07
feature index 8
feature Index in
20
feature index
R*
0.06
30
0.05
8.6
0.04
8.4
0.03
70
8.2
70
model ranking (R", terms)
model ranking (R², terms)
model ranking (R*, terms)
An example of how candidate models with m =3=3, 3, m= and m= terms with largest R² in each category tend to consistently choose the true terms the model. The indices of all terms in the dictionary (the basis functions k, (a)) are shown on the left. vertical axis. Each candidate model along the horizontal axis is represented by squares indicating which terms make up the respective model. The ground truth model in this example has 3 lerms, indicated by two lighter squares the left and the horizontal lines. The case where the candidate model is the true model is indicated by vertical line in the middle plot. The R² value in logarithmic scale is shown as line with closed circles, the values are given by the right vertical axis. It. can be seen that the true model is chosen by I22 among all models with 3 terms.
Figure 4:
0.570
10
0
10
feature index n
20
0.565
feature index n
M)
20
30
30
0.560
40
40
In
50
50
0.555
60
60
70
70
0.550
model ranking (top R²)
model ranking (CS-p(M)
Figure 5: In the plot on the left, the f best models in terms of R2 from each of the m-sized candidate models in figure have been combined and sorting according to their n² value. Clearly, the models with more terms are favored over models with fewer terms, and the true model is not selected. The plot on the rights shows the best models in terms of the Bayesian model evidence P(ylM), where the true model is selected illustrating the adequate overfitting penalty held by M)y)
24

---

Published in Transactions on Machine Learning Research (11/2023)
CS CP-P(M) M
M&M
M
top top
top
top
..e
count kula))
rating
selected
Ka(a))
removed
1f Dew selected
if selected selected
new selocted
top
top
top
top
RR
R2
R
R2
2
m==3
candidate models
Figure 6: Schematic representation of algorithm CS-R² and CS-pM. Starting point is the set of candidate models built from basis functions (x). For fixed model size m, R² is computed for all (*") candidate models. Incrementing m, top models in terms R² are collected, from which basis functions k, are rated based on counts of (x) these top models. Typically, incrementing m, new k, (s) with high rates are found and selected for the inferred model M by algorithm CS-IR" as long as m is smaller than the true model size. The iteration in therefore terminates if no new (s) selected and M is returned. Basis functions (n that are hardly selected removed and not considered for building candidate models in following iterations. The top models selected by R² are scored again by p(M) of which the best model (!)) is the put algorithm CS-p(M).
25

---

Published in Transactions on Machine Learning Research (11/2023)
The final aspect to be discussed here is the complexity of your CS approach. As an empirical method with unpredictable stopping criterion and pruning, rigorous analysis appears intractable. However, the complexity of a direct computation of (2( ~ y"K K (K"K) K"y y with N m matrix K and N x 1 vector c.f. Eq. 6 can be determined as follows: The matrix product KJK has complexity Omm²N), matrix-vector products K and KTy y have complexity O(mN), the malrix-inverse of KT has complexity O(m²), such that in total we have a complexity of (m² + N)). Luckily, as discussed in RUscr4512 2018)) efficient. numerical implementations can bring it down O(m (m....Nu77) i.e. near-linear scaling, which confirmed with an own numerical analysis shown in Fig. 7
26

---

Published in Transactions on Machine Learning Research (11/2023)
(a)
=
N=100
0.14
average runtime per model [ms]
N 400
0.12
1000
2000
0.10
0.08
0.02
4
7
candidate model size
N
= 100
10a
N - 400
number candidate models
1000
(P)
107
105
10a
candidate model size m
(c)
N-100
average runtime per iteration [min]
102
- 200
N 1000
- 2000
10-1
-0
10-1
candidate model size
Figure 7: Runtime analysis of the R2 elimination the comprehensive search approach based on 100 random polynomials and p 72 basis functions on standard laptop using about cores. Chart (a) shows that the computation of n22 per model stays well below 1 ms with slight increase with model size Despite the near-exponential increase of model space dimension shown in chart (b), the efficient computation of R² and the effect of pruning keeps the runtime per iteration below an exponential increase and for practically relevant model sizes of m < well below a minute.
27


Reviews:
[{"review_id": "eJqVVryhcn", "reviewer": "Reviewer_F8VX", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strengths include\n\n(i) clarity and conciseness of exposition \\\n(ii) an interesting new angle of incorporating the Bayesian evidence. \n\nThe weaknesses include \n\n(i) repetition of text between the manuscript and appendix \\\n(ii) incomplete or naive arguments (e.g. the authors claim to reduce the complexity of the sparse regression problem by only considering all models of sparsity $m<p$, so considering ${p \\choose m}$ total models, instead of all $2^p$ models, when the former is still NP-hard) \\\n(iii) While it is useful to identify $R^2$ as a criterion, this is by no means a new approach, and the authors heavily rely on it for their method, reducing the novelty of the manuscript.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present a Bayesian version of the system identification problem employing classical notions from statistics. In particular, the authors modify the SINDy method to include a forward-backward stepwise regression algorithm to find the coefficient weights, which emplyed the $R^2$ statistic and approximation of the model evidence under assumptions of normality.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "S denotes merely stengthens paper, C denotes critical for acceptance.\n\n- S: (Sec 2) Nardini et al are incorrectly cited they use NN to denoise the data and STLSQ for sparse regression. i.e. NNs are not employed for model selection.\n- S: (Sec 2) the \"improved noise robustness\" category should include reference to weak formulations of SINDy (e.g. WSINDy)\n- C: (Sec 2) Review of Bayesian-based sparse regression is lacking, e.g. \"Robust data-driven discovery of governing physical laws with error bars\", Zhang, Lin\n- S: Eq 1: for identifying dynamical systems one rarely has an iid normal error structure. The authors should comment on this\n- S: (after eq 6, and throughout) it should be made more clear that $k_n(\\pmb{x})$ refers to a set of functions\n- C: beginning of sec. 3.4, \"a non-greedy algorithm would consider...\", this is slightly incorrect, convex relaxations of sparse regression problems are not greedy and do not exhaustively consider all possible sparsity patterns. The authors should clarify what is meant by \"greedy\" vs \"non-greedy\", as typically sparse regression algorithms that build the solution term by term are referred to as greedy (e.g. OMP), then there are thresholding methods (STLSQ, IHT), and convex relaxations (LASSO). Thresholding is not typically considered a greedy approach as it is here.\n- C: (sec 3.4) Considering all ${p \\choose m}$ options is still a combinatorially hard problem. For instance, $p=20$, $m=5$ gives ${p \\choose m}=15504$. In fact, many sparse regression algorithms assume that $m$ is a user input (see OMP, COSAMP). Since this part of the algorithm will not scale well, the authors should comment on how to adapt it to more realistic settings. Is the \"additional pruning step\" meant to circumvent this? If so more explanation is needed. How many iterations of not being selected are allowed before the term is omitted completely? Can it ever be regained?\n- C: Compared to the algorithms in the appendix (which should be listed in the main body) section 3.4 is lacking key details. In contrast, Appendix C is missing any mention of this key \"pruning step\", so there appears to be a mismatch between the algorithm as stated in the body and in the appendix.\n- S: much of appendix A is a repeat of the manuscript\n- S: the switch from __y__ to $y$ in eq. 18-22 is unmotivated\n- S: \"is essentially is\" after eq 22\n- C: after eq 25, \"we see that the evidence p(y) is in fact the model-likelihood p(y|M)\" this is not clear to me, the authors should explain why p(y) equals p(y|M). In fact, if this is true then the posterior equals the prior always.\n- C: after eq 25, in general the authors should define $\\mathcal{M}$. Is $\\mathcal{M}$ a set of terms, so an indicator vector in {0,1}$^p$ (equal to the sparsity pattern of $\\pmb{w}$)? Or is it the entire set of possible functions that are used in the regression? Is it a function space? A \"constant prior\" $p_{pr}(\\mathcal{M})$ may not make sense in some of these interpretations. \n- S: eq 29-31 are repeated in the line above, authors should pick one.\n- S: after eq 31, \"...is $M$ is...\"\n- S: much of the appendix should be given in the rest of the paper, for instance the discussion in eq 32-35 etc.\n- C: eq 35 and mention of finite differences throughout should include mention of methods to circumvent the errors of this approach (namely using a weak formulation)\n- S: \"The problem is that basis functions like $1/(1+x^n)$ are divergent for certain values of x, and are also quite limiting in their form\", perhaps the authors should say \"certain values of $n$\" since $n$ even leads to bounded functions. While this notion is interesting, I would caution away from claiming these functions are \"limiting\" (they appear in many biochemical models)\n- S: \"needs to be determined from a 1-dimensional numerical root-finding algorithm\" this is interesting, could the authors provide the details? Conditioning would seem to be a concern here.\n- S: $y^T$ should be bold in equation 45", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bVqdL9c6rJ", "reviewer": "Reviewer_ihmg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper has many weaknesses.\n\n1 The proposed algorithm only works in the model class containing the true model. This is hardly the case.\n\n2 The proposed LG algorithm has the same order complexity as a full search. Because the binomial (p p/2) also grows exponentially in p. This algorithm would not work for larger p or a solution that contained many non-zero entries. \n\n3 Lasso is as good as the proposed method for the Lorenz attractor. \n\n4 For the toy data in Section 4.1, I am surprised that LASSO performs so poorly. I have the impression that lambda is set too high. But the authors do not say how they adjust lambda. \n\n5 The R2 is used for model selection. But the authors do not say whether test or training data is used for computing R2. If it is the former, how much data does the test set contain? If it is the latter, how do they avoid overfitting? \n\n6 The authors talk about feature selection. But there are many relevant missing literature studies. For example https://www.jmlr.org/papers/special/feature03.html and https://bayes.wustl.edu/MacKay/pred.pdf.\n\n7 The paper is mostly written in a frequentist manner, but it uses some concepts from Bayesian statistics that do not make sense. For example, the solution for the regression problem uses regularized least squares, but the author proposes a Bayesian model selection procedure. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors propose a not-so-greedy approach for model selection. The algorithm only works when the true model is a linear combination of basis functions and the number of non-zero weights is small. \n\nThe authors propose two algorithms: SR and LG. SR first follows a forward procedure to include potential basis functions and then uses a backward step. The LG algorithm, which is the strongest algorithm, only works for polynomials. It first takes all the polynomials of degree one and keeps the best model. To the obtained model, it adds all the polynomials of degree two and keeps the best model.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "I do not think this paper is interesting to the ML community and no change in it would make me reconsider this decision. For me, this paper should have been an editorial rejection. ", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9fGoPRZxIu", "reviewer": "Reviewer_mC8n", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n* The idea of the paper is straightforward: a heuristic algorithm is proposed based on comparing existing methods, and experiments are conducted to compare the new method with previous ones.\n\nWeakness:\n* The writing can be improved. \n* If you are comparing with LASSO and SR, I don't understand why there is a need to discuss those basis functions instead of just focus on linear regression. \n* The R2 related discussion seems quite obvious and relatively redundant. \n* \"Greedy\" v.s. \"Less greedy\" doesn't seem like an intriguing conceptual comparison to me. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a new method for sparse linear regression based on Bayesian model selection, inspired by stepwise regression. Essentially, the algorithm does some less greedy operation when selecting features since it focuses on a larger search space. The algorithm is purely heuristic based and the authors provided some experiments that shows they are better than all existing ones. ", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "I would suggest the following:\n* Write the algorithm more formally (even it is heuristic, it should be written in a semi-pseudo algorithm format).\n* Analyze the computation complexity further. \n* Explain more on the comments related to existing literature: e.g., \"Greedy algorithms have proven to be more successful than LASSO in equation learning\" -- where is the proof? and is it suggesting that SR is always better than LASSO for linear regression? ", "broader_impact_concerns": "No specific broader impact concerns. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers were a bit split: one was tentatively positive, another found merit in the subfield of equation learning but had many comments that they felt were most suited for a major revision, and the third found the assumptions of equation learning quite restrictive (i.e., that the true model is in the class, and related to this, that the residual was iid Gaussian).

I read the paper myself, keeping the TMLR [evaluation criteria](https://jmlr.org/tmlr/editorial-policies.html) in mind:
- *Are the claims made in the submission supported by accurate, convincing and clear evidence?*
- *Would at least some individuals in TMLR's audience be interested in knowing the findings of this paper?*

*Papers should be accepted if they meet the criteria, even if the contribution or significance of the work is modest*.  

Valid grounds for rejection are lack of evidence, incorrectly claimed novelty, and re-implementing an idea that has been reproduced before.  I didn't see any of these issues.

Under these basic criteria, I find the paper does pass the threshold.  This is a new idea, and there is numerical evidence to substantiate some of the claims.

However, the reviewers raised good points, and in particular reviewer F8VX requested significant changes.  Hence I am recommending an accept-with-revisions.
- Please address reviewer F8VX's comments
- There are some formatting issues, for example `\citet` and `\citep` are mixed up for many of the references.
Decision: Accept (poster)