GSR: A Generalized Symbolic Regression Approach

As highlighted in Section 2, the goal of SR is to search the function space to find the model that best fits the mapping between the independent variables and the target variable (i.e. the mapping between $\mathbf{x}_{i}$ and $y_{i}$, for all $i$). Since the main objective of SR is to recognize correlations and find non-trivial interpretable models (rather than making direct predictions), we modify the goal of SR; we instead search the function space to find the model that best describes the mapping between the independent variables and a transformation of the target variable (i.e. the mapping between $\mathbf{x}_{i}$ and some transformation or function of $y_{i}$, for all $i$). Formally, we propose modifying the goal of SR to search for appropriate (model) functions from a space of all possible mathematical expressions $\mathcal{S}$ defined by a set of given mathematical functions (e.g., $\exp$, $\ln$, $\sin$, $\cos$) and arithmetic operations (e.g., $+$, $-$, $\times$, $\div$), which can be described by the following optimization problem:

where $f^{*}$ and ${\color[rgb]{0,0,1}g^{*}}$ are the optimal analytical functions. In other words, instead of searching for mathematical expressions of the form $y=f(\mathbf{x})$ as is usually done in the SR literature, the proposed GSR approach attempts to find expressions of the form $g(y)=f(\mathbf{x})$. We illustrate this concept in Table 1.

Although the main goal of GSR is to find expressions of the form $g(y)=f(\mathbf{x})$, we may encounter situations where it is best to simply learn expressions of the form $y=f(\mathbf{x})$ (i.e. $g(y)=y$). For instance, consider the ground truth expression $y=\sin(x_{1})+2x_{2}$. In this case, we expect to learn the expression exactly as is (i.e. $g(y)=y$ and $f(\mathbf{x})=\sin(x_{1})+2x_{2}$) as long as the right basis functions (i.e. $\sin(x)$ and $x$ in this case) are within the search space, as we will see in the next sections.

Making predictions. Given a new input feature vector $\mathbf{x}_{*}$, predicting $y_{*}$ with GSR is simply a matter of solving the equation $g(y)=f(\mathbf{x}_{*})$ for $y$, or equivalently, $g(y)-f(\mathbf{x}_{*})=0$. Note that $f(\mathbf{x}_{*})$ is a known quantity and $y$ is the only unknown. If $g(\cdot)$ is an invertible function, then $y_{*}$ can be easily found using $y_{*}=g^{-1}\big{(}f(\mathbf{x}_{*})\big{)}$. If $g(\cdot)$ is not invertible, then $y_{*}$ will be the root of the function $h(y)=g(y)-f(\mathbf{x}_{*})$. Root-finding algorithms include Newton’s method. Whether the function $g(\cdot)$ is invertible or not, we might end up with many solutions for $y_{*}$ (an invertible function, which is not one-to-one, can lead to more than one solution). In this case, we choose $y_{*}$ to be the solution that belongs to the range of $y$ which can be determined from the training dataset.

Now that we have presented the goal of our proposed GSR approach (summarized by Equation 2), we need to constrain the search space of functions $\mathcal{S}$ to reduce the computational challenges and accelerate the convergence of our algorithm. Inspired by mcconaghy2011ffx; chen2017elite as well as classical system identification methods (brunton2016discovering), we confine $\mathcal{S}$ to generalized linear models, i.e. to functions that can be expressed as a linear combination (or as a weighted sum) of basis functions (which can be linear or nonlinear). In mathematical terms, for a given input feature vector $\mathbf{x}_{i}$ and a corresponding target variable $y_{i}$, the search space $\mathcal{S}$ is constrained to model functions of the form:

where $\phi_{j}(\cdot)$ and $\psi_{j}(\cdot)$ are the basis functions applied to the feature vector $\mathbf{x}_{i}$ and the target variable $y_{i}$, respectively, $M_{\phi}$ and $M_{\psi}$ denote the corresponding number of basis functions involved, respectively. In matrix form, the minimization problem described in Equation 2 is equivalent to finding the vectors of coefficients $\bm{\alpha}=[\alpha_{1}\,\cdots\,\alpha_{M_{\phi}}]^{T}$ and $\bm{\beta}=[\beta_{1}\,\cdots\,\beta_{M_{\psi}}]^{T}$ such that:

where

Note that if we examine the minimization problem as expressed in Equation 4, we can indeed minimize $||\mathbf{X}\bm{\alpha}-\mathbf{Y}\bm{\beta}||^{2}$ by simply setting $\bm{\alpha}^{*}=0$ and $\bm{\beta}^{*}=0$ which will not lead to a meaningful solution to our GSR problem. In addition, to avoid reaching overly complex mathematical expressions for $f(\cdot)$ and $g(\cdot)$, we are interested in finding sparse solutions for the weight vectors $\bm{\alpha}^{*}$ and $\bm{\beta}^{*}$ consisting mainly of zeros which results in simple analytical functions containing only the surviving basis functions (i.e. whose corresponding weights are nonzero). This is closely related to sparse identification of nonlinear dynamics (SINDy) methods (brunton2016discovering). To this end, we apply $L_{1}$ regularization, also known as Lasso regression (tibshirani1996regression), by adding a penalty on the $L_{1}$ norm of the weights vector (i.e. the sum of its absolute values) which leads to sparse solutions with few nonzero coefficients. In terms of our GSR method, Lasso regression automatically performs basis functions selection from the set of basis functions that are under consideration.

Putting the pieces together, we reformulate the minimization problem in Equation 4 as a constrained Lasso regression optimization problem defined as

where $\lambda>0$ is the regularization parameter, and

To solve the GSR problem, we first present our approach for solving the constrained Lasso problem in Equation 6, assuming some particular sets of basis functions are given. We then outline our genetic programming (GP) procedure for finding the appropriate (or optimal) sets of these basis functions, before discussing our matrix-based encoding scheme (to represent the basis functions) that we will use in our GP algorithm.