Interaction-transformation symbolic regression with extreme learning machine

Abstract

Symbolic Regression searches for a mathematical expression that fits the input data set by minimizing the approximation error. The search space explored by this technique is composed of any mathematical function representable as an expression tree. This provides more flexibility for fitting the data but it also makes the task more challenging. The search space induced by this representation becomes filled with redundancy and ruggedness, sometimes requiring a higher computational budget in order to achieve good results. Recently, a new representation for Symbolic Regression was proposed, called Interaction-Transformation, which can represent function forms as a composition of interactions between predictors and the application of a single transformation function. In this work, we show how this representation can be modeled as a multi-layer neural network with the weights adjusted following the Extreme Learning Machine procedure. The results show that this approach is capable of finding equally good or better results than the current state-of-the-art with a smaller computational cost.

Introduction

Regression analysis [1] is a process used to describe the distribution of a dependent variable, also called target or outcome, given the values of a set of independent variables, called predictors. This relationship can then be used for prediction, the understanding of the system being studied, or even the inference of causal relationships. When describing the regression models ¹, $\mathbf{X}$ is the $n\times d$ matrix of n samples, each one composed of d predictors and ${\mathbf{X}}_i$ is the i-th sample. The corresponding target values of each sample is described by a vector of n elements called $\mathbf{y}$. A regression model is often described as:

$$y_i=f({\mathbf{X}}_i)+{∊}_i,$$

where $f({\mathbf{X}}_i)$ is a function of the predictors that outputs a corresponding target value and ${∊}_i$ is a residual or normally distributed random error centered on 0 and with a variance ${\sigma }^2$. A simple example of a regression model is the linear model:

$$y_i=f_{\mathit{linear}}({\mathbf{X}}_i)+{∊}_i={\beta }_0+{\beta }_1·x_1\dots {\beta }_d·x_d+{∊}_i,$$

with $x_j$ representing the j-th variable for the sample ${\mathbf{X}}_i$. This model is often called interpretable since the impact of each predictor to the target value is clearly described by $\beta$. With this model, a regression algorithm has the task of finding the optimal $\beta$, also called free parameters, that minimizes the prediction error.

As with the linear regression model, many regression algorithms [2] try to adjust the free parameters of a regression model described by a closed function form. Examples of such models are Polynomial Regression [3] and Multilayer Perceptron [4], [5]. Depending on the complexity of the model, we may lose the interpretation capabilities. In such cases, the model is termed a black-box model.

A different approach to regression analysis is given by Symbolic Regression [6], [7], [8], [9], as described in Section 2, where both the function form and the free parameters are adjusted concomitantly. This model is represented by an expression tree and this task is often performed with the use of Genetic Programming algorithms [6], [10], an evolutionary meta-heuristic that evolves a population of expressions. The main motivation behind this approach is that the returned model is expected to be close or equal to the generator function of the studied process, thus improving the understanding of the system.

The main challenge of applying this algorithm is that the explored search space is huge, rugged, noisy, and with many local optima. This makes the performed search inefficient and often leading to bloat [11], [12], [13], [14], i.e. when the generated models are comparable to a black-box model w.r.t. interpretability.

Recently, a new representation for Symbolic Regression was proposed in [15], named Interaction-Transformation (IT) expression, which will be detailed in Section 2. This representation restricts the search space by removing function forms that can lead to bloating. An Interaction-Transformation function is a function of the form:

$$f({\mathbf{X}}_i)={\beta }_0+{\displaystyle \sum _{j=1}^n}{\beta }_j·g_j(p_j({\mathbf{X}}_i)),$$

with the j-th interaction function described as:

$$p_j({\mathbf{X}}_i)={\displaystyle \prod _{k=1}^d}{x}_k^{s_{\mathit{jk}}},$$

where $g_j:{\mathbb{R}}^d\to \mathbb{R}$ is called a transformation function, $p_j({\mathbf{X}}_i)$ is an interaction function, $s_{\mathit{jk}}$ is the strength of the interaction of the k-th predictor for the j-th interaction, and n is the number of additive terms of the expression.

In this same paper, the author introduced the algorithm called SymTree as a simple search heuristic that incrementally expands the current best expression by adding new interactions or experimenting with different transformation functions. Some experiments performed on this paper with synthetic benchmark functions pointed out that SymTree was capable of finding regression models closer to the original generating function than those found by traditional nonlinear and symbolic regression algorithms. This study was further extended in [16] to 20 different equations from physics and engineering resulting in the discovery of the correct generating function in 14 cases and a similar function in 5 of them. One problem with the SymTree algorithm is that, in the worst case scenario, the complexity for t steps of the algorithm is $O(d^{2t})$, for a d-dimensional data.

Following Eq. 1, we can model a layered computational model, with the first layer adjusting the strength of interactions ($s_{\mathit{jk}}$), the second layer just serving as a pass-through to the final layer, called transformation layer, in which the free parameters ${\beta }_j$ are adjusted. In this context, we want to find $\mathbf{S},\beta$ such that:

$$y_i=f({\mathbf{X}}_i,\mathbf{S},\beta )+{∊}_i,$$

minimizes the approximation error.

With this perceptron model, we can fix the size of the expression and adjust the free parameters. But, to use a gradient descent algorithm to find the optimal values, we depend on the gradient of the error function that, for the model described in Eq. 1, is undefined if we have any $x_j⩽0$.

As such, in this paper, we propose the use of an Extreme Learning Machine [17], [18] approach with the values of the parameters of the first layer being generated at random and the values of the last layer adjusted by a gradient descent algorithm with $l0$ and $l1$ regularization. The details of this new approach will be explained in Section 3.

A set of experiments is described in Section 4 to verify the performance of different variations of the proposed algorithm and compare them with the current state-of-the-art in Symbolic Regression. The results reported in Section 5 show that the proposed approach outperforms the state-of-the-art while being up to 30 times faster than its closest competitor. Following the encouraging results, in Section 6 we describe some ideas for further improvements together with some final discussions about the proposed algorithm.