Freek Stulp
email: freek.stulp@ensta-paristech.fr

Olivier Sigaud
email: olivier.sigaud@isir.upmc.fr

Many Regression Algorithms, One Unified Model -A Review

Keywords: regression, locally weighted regression, Gaussian mixture regression, radial basis function networks, Gaussian process regression

. The aims of this paper are to provide an overview of many regression algorithms, and to demonstrate how the function representation whose parameters they regress fall into two classes: a weighted sum of basis functions, or a mixture of linear models. Furthermore, we show that the former is a special case of the latter. Our ambition is thus to provide a deep understanding of the relationship between these algorithms, that, despite being derived from very different principles, use a function representation that can be captured within one unified model. Finally, step-by-step derivations of the algorithms from first principles and visualizations of their inner workings allow this article to be used as a tutorial for those new to regression.

Introduction

Regression is the process of learning relationships between inputs and continuous outputs from example data, which enables predictions for novel inputs. This relationship is represented as a function f : X → Y , which predicts, for instance, a person's height from their age. Here, the input space X (age), is known as the independent variable, and the output space Y (height) as the dependent variable [START_REF] Fisher | Statistical methods for research workers[END_REF]. In the example, the training data consists of concrete age and height measurements for a set of people.

Regression is a form of supervised learning where the output space is continuous, i.e. Y ⊆ R M . In parametric regression, one assumes that the function f is well represented by a specific parameterized model, for instance a linear model f (x) = a x. With a linear model, the model parameters are the slopes a. The aim of parametric regression is to find the parameters of the model that minimize some error on the training examples.

Another example of a parameterized model is a Radial Basis Function Network [START_REF] Park | Approximation and radialbasis-function networks[END_REF], where the function is modeled as a weighted sum of basis functions f (x) = E e=1 w e φ e (x). If we assume that the basis functions φ e=1...E have pre-specified centers and widths, the model parameters that are to be determined through parametric regression are the weights w e=1...E , see Figure 1.

The generic scheme for parametric regression is depicted in Figure 2. The input to the regression algorithm is the training data and a set of algorithmic meta-parameters, including for instance learning rates. Each regression algorithm assumes a certain type of model, e.g. linear least squares assumes a linear model. The output of the algo- rithm is a vector of model parameters, which are determined by minimizing an error measure on the training data. Evaluating the model to make predictions for novel inputs requires both the model (e.g. f (x) = a x) and its model parameters (e.g. a = [2 1] ). A detailed discussion of the differences between model parameters and meta-parameters is given in Section 2.6. In this article, we take a model-centric view on regression, which means that we classify and analyze algorithms based on the model they assume, rather than the algorithmic procedure that is used to optimize the parameters of this model. Our first contribution is to show that the models used in a wide variety of regression algorithms (listed in Table 1) fall into two main classes: a mixture of linear models or a weighted sum of basis functions.

Our second contribution is to demonstrate that the latter class of models (weighted sum of basis functions) is a special case of the former one (mixture of linear models). As a consequence, and rather strikingly, all the algorithms in Table 1 despite having being derived from very different principles -use parameterized functions that can be described by one unified model. This has been visualized in Figure 3. Thus, these regression algorithms should not be thought of as using their own distinct model customized to the algorithmic procedure, but rather as using models that are spe-cial cases of the unified model. Such a perspective provides a deeper understanding of the relationship between these algorithms, and is a necessary step towards model-based machine learning, as proposed by [START_REF] Bishop | Model-based machine learning[END_REF], i.e. the idea of the automated selection of the adequate machine learning algorithm given the formal description of a specific learning problem.

training data Figure 3: The algorithms described in this paper have a model that is either a mixture of linear models (upper block) or a weighted sum of basis functions (lower block). Since the former is a special case of the latter, the output of all algorithms can be described by one unified model. Full algorithm names are given in Table 1.

(
Despite our model-centric view, we do describe and explain the algorithmic procedures used in different regression algorithms, including (regularized) least squares, expectation-maximization, backpropagation, decision tree learning, and Gaussian process regression. This is necessary to understand why an algorithm assumes a certain type of model, and how that model relates to the unified model we propose. These explanations however, should not distract from the fact that our main interest is in the underlying model that the algorithms assume (linear model, rbfn, model tree, Gaussian mixture model, Gaussian process), and that all these models are special cases of the unified model.

Explaining the algorithms also allows this article to be used as a tutorial on regression; we provide an overview of many algorithms, show their derivations from first principles, visualize their inner workings so that novices may acquire an intuitive understanding, and provide network representations for readers with a background in artificial neural networks. Using one unified, easy to Model: mixture of linear models lwr Locally Weighted Regression (Atkeson and Schaal, 1995) rfwr Receptive Field Weighted Regression [START_REF] Schaal | Receptive field weighted regression[END_REF]) lwpr

Locally Weighted Projection Regression [START_REF] Vijayakumar | Locally weighted projection regression: An O(n) algorithm for incremental real time learning in high dimensional space[END_REF] xcsf . [START_REF] Butz | Context-dependent predictions and cognitive arm control with XCSF[END_REF]) gmr Gaussian Mixture Regression [START_REF] Hersch | Dynamical system modulation for robot learning via kinesthetic demonstrations[END_REF][START_REF] Calinon | Robot programming by demonstration[END_REF]

M5

Model Trees [START_REF] Quinlan | Learning with continuous classes[END_REF] Model: weighted sum of basis functions rbfns Radial Basis Function Networks [START_REF] Park | Approximation and radialbasis-function networks[END_REF]) krr Kernel Ridge Regression [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF]) gpr

Gaussian Process Regression [START_REF] Williams | Gaussian processes for machine learning[END_REF]) svr Support Vector Regression (Vapnik, 1995) irfrls Incr. Random Features Regularized Least Squares [START_REF] Gijsberts | Incremental learning of robot dynamics using random features[END_REF]) i-ssgpr

Incr. Sparse Spectrum Gaussian Process Regr. [START_REF] Gijsberts | Real-time model learning using incremental sparse spectrum gaussian process regression[END_REF]) CART Regression Trees [START_REF] Breiman | Classification and Regression Trees[END_REF]) elm Extreme Learning Machine (Huang et al., 2006b) BProp Backpropagation [START_REF] Werbos | Beyond regression: New tools for prediction and analysis in the behavioral sciences[END_REF] Table 1: List of algorithms presented in this article.

understand model highlights relationships between algorithms; a key to acquiring more global understanding of regression methods. It is not our aim to be exhaustive, in terms of presenting all regression algorithms and their variants. This would distract from our actual aim, which is to highlight the similarities and differences between those algorithms whose underlying model is a special case of the unified model. For further reading, we provide references to other tutorials and books, including the work of [START_REF] Smola | A tutorial on support vector regression[END_REF], [START_REF] Bishop | Pattern recognition and machine learning[END_REF], [START_REF] Williams | Gaussian processes for machine learning[END_REF], and [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF].

The rest of this article is structured as follows. In the next section, we present (regularized) linear least squares, the classical algorithm for linear regression. In Section 3, we describe how linear least squares has been extended in two different directions. These two directions have lead to two classes of algorithms, which generate either a mixture of linear models or a weighted sum of basis functions. In Section 3.3, we show that these two different classes of models can in fact be unified, the latter being a special case of the former. Algorithms that use these two classes of models are presented in Section 4 and 5 respectively. In Section 6, we discuss the impact of the unifying perspective, before concluding with Section 7.

Linear Least Squares

In this section, we present linear least squares (lls) regression. We describe some variants of the basic algorithm, including different regularization methods, multivariate lls, and recursive least squares. The aim of the following section -Section 3 -is then to show how linear least squares is used in the context of non-linear regression.

First, we formalize some of the concepts related to parametric regression, as presented in the introduction. The aim of supervised learning is to learn a function f :

X → Y from N training examples 1 {(x n , y n )} N n=1
, where ∀n, x n ∈ X ∧ y n ∈ Y . In regression, the output space is continuous (Y ⊆ R M ). Multivariable (resp. multivariate) regression refers to input spaces X (resp. output spaces Y ) with a dimensionality of more than one2 .

In the case of linear regression, the family of model functions is linear, which means that f (x) is represented as a line (1D), plane (2D) or hyperplane (≥3D). Formally, if the data is assumed centered so that the hyperplane goes through the origin, the underlying model is

f (x) = a x, ( 1
)
where a is a D × 1 vector that contains the slopes of the linear function, as visualized in Figure 4. A popular method for determining a is the linear least squares (lls) algorithm. It uses the vertical distance between observed values y n and the predictions f (x n ), which are known as the residuals:

r n = y n -f (x n ) (2) = y n -a x n .
(for linear models)

The residuals are visualized as vertical lines in Figure 5. In the lls case, the sum of the squared residuals

S(a) = N n=1 (y n -f (x n )) 2 (4) = N n=1 (y n -a x n ) 2 (5)
is minimized, which in matrix form is

S(a) = (y -Xa) (y -Xa), (6)
with X being the N × D design matrix that combines all the input vectors, and the N × 1 vector y combining all scalar target values:

X =      x 1,1 x 1,2 • • • x 1,D x 2,1 x 2,2 • • • x 2,D . . . . . . . . . . . . x N,1 x N,2 • • • x N,D      , y =      y 1 y 2 . . . y N      , (7)
where each row corresponds to one input/target example. The slopes of the linear model that minimize (6) are

a * = arg min a S(a) (8) = arg min a (y -Xa) (y -Xa), (9)
which is a continuously differentiable unconstrained optimization problem. The solution of this problem is to determine when the derivative

S (a) = 2(a(X X) -X y) (10)
is 0, which is the case when a * = (X X) -1 X y.

(11)

Thus, given the input data in X and target data in y, the linear function that best fits the data -in terms of minimizing the sum of squared residualsis f (x) = a * x.

Before continuing, we want to draw explicit attention to the difference between the model used in this section -the linear model in (1) -and the algorithm used to determine the model parameters of the model -the computations in (11). Note that to make a prediction for a novel input x new , we need only to know the model (linear) and its parameters (a), but not the algorithm used to compute a.

Least Squares and Least Deviations

The L p -norm z p of a D-dimensional vector z is defined as

z p = D d=1 |z d | p 1 p . ( 12
)
The L 2 -norm is commonly known as the Euclidean distance, and the L 1 -norm as the Manhattan distance.

As the name implies, linear least squares uses the sum of squares of the residuals as the error measure to minimize. This corresponds to the square of the L 2 -norm of the residuals, and we may rewrite (9) as

a * = arg min a ( y -X a 2 ) 2 (13) = arg min a    N n=1 |y n -a x n | 2 1 2    2 (14) = arg min a N n=1 |y n -a x n | 2 . ( 15
)
Alternatively, one may minimize the L 1 -norm of the residuals

a * = arg min a y -X a 1 (16) = arg min a N n=1 |y n -a x n | . ( 17
)
which is known as "least absolute deviations" [START_REF] Bloomfield | Least absolute deviations curve-fitting[END_REF]. Regression with the L 1 -norm is more robust, but a disadvantage is that ( 16) cannot be solved analytically, so an iterative optimization approach is required.

Regularized Least Squares

Potential singularities in (X X) may make it difficult to invert this matrix, which can result in very large values in the vector a. A solution to this issue is to explicitly penalize large weights, which results in Regularized Linear Least Squares (rgls), where (9) becomes

a * = arg min a λ 2 a 2 + 1 2 y -X a 2 , ( 18
)
where λ is a parameter that determines the tradeoff between small values in a and small residuals.

Using the L 2 -norm for a

a 2 = D d=1 |a d | 2 1 2 (19)
is called Thikonov regularization. rgls with Thikonov regularization is known as Ridge Regression (rr). Using the same derivation as for ( 11), the analytical solution to this least squares minimization problem is

a * = (λI + X X) -1 X y. (20)
It is also possible to use the L 1 -norm

a 1 = D d=1 |a d |, (21)
which has the advantage that values in a will tend to go to zero. In the context of basis function networks (Section 5), the L 1 norm can thus be used as a feature selector. Here again, using the L 1 -norm does not have an analytical solution as in (20). Several algorithms for finding solutions are described in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Schmidt | Least squares optimization with l1-norm regularization[END_REF]. In the context of rgls, the L 1 -norm is known as the "least absolute shrinkage and selection operator" (LASSO).

Another norm is used in Support Vector Regression (svr) methods, which we discuss in Section 5.5. These three norms are visualized in Figure 6.

Adding offsets to the linear model

In the above, we have assumed that the data is centered, in which case a linear model f (x) = a x suffices. If the data is not centered, an offset should be added so that the linear model becomes f (x) = ]. In this case, the design matrix becomes

X =      x 1,1 x 1,2 • • • x 1,D 1 x 2,1 x 2,2 • • • x 2,D 1 . . . . . . . . . . . . 1 x N,1 x N,2 • • • x N,D 1      (22)
The result of fitting a linear model with offsets to our example data -which are sampled from the function y = 0.5x + sin(x) -is depicted in Figure 5, where the estimated function is f (x) = 0.11x+1.18.

Multivariate regression

In multivariate regression, the output space is multi-dimensional (Y ⊆ R M ), where the output is thus a vector, a becomes a matrix A, and the linear model is

f (x) = A x. ( 23
)
The design matrix of inputs is the same matrix X, but the targets are now stored in a matrix Y, and the least squares solution becomes

A * = (X X) -1 X Y. ( 24
)
This solution is decoupled between the different output variables, and corresponds to M separate univariate least squares solutions:

a * m = (X X) -1 X y m . (25)
In the rest of this paper, we will mostly use a 1-D output, for ease of presentation and visualization. All examples are easily extended to the multivariate case by applying (24) or (25).

Batch vs. Incremental Learning

Equations ( 11) and ( 20) provide a batch method for finding the optimal weights w * , by considering

the training examples all at once. In contrast, incremental methods are able to update the weights (and other parameters of the model) incrementally, one training example at a time.

Throughout this paper, we present several models for which both batch and incremental algorithms exist. For instance, the weights of a single-layer feedforward network can be trained through Extreme Learning Machines (batch) or backpropagation (incremental). The weights of a mixture of linear models can be updated through locally weighted regression (batch) or receptive field weighted regression (incremental).

For least squares, the batch method in ( 11) and (20) requires the inversion of a matrix, which has complexity O(N 3 ). In fact, incremental methods also exist for linear least squares itself. The most well-known is the Recursive Least Squares (rls, [START_REF] Plackett | Some theorems in least squares[END_REF], which updates the linear model from each new example without recomputing (11) at each step.

Alternatively, one can use the Sherman-Morrison formula, which enables the inverse to be updated incrementally, reducing the inversion complexity to O(N 2 ). However, this method is sensitive to rounding errors; a numerically more stable option consists in updating the Cholesky factor of the matrix using the qr algorithm [START_REF] Gijsberts | Real-time model learning using incremental sparse spectrum gaussian process regression[END_REF].

Another incremental approach to solving these equations consists in optimizing the weights of a matrix that stands for (X X) -1 using gradient methods, as in the backpropagation algorithm for instance (see Section 5.10).

Note that the distinction between batch and incremental learning is related to, but not the same as the distinction between offline and online learning. A batch algorithm uses all the training data at once. An incremental algorithm updates an existing model incrementally, using one training example at a time. In offline learning, all examples are stored and accessible, whereas in online learning a training example is forgotten once it has been processed. Batch learning is always off-line, because it needs access to all the training examples at once, and they must thus be stored in memory. Online training always requires an incremental algorithm, because only one data example can be considered at a time, and the algorithm must thus be able to update the model incrementally. However, incremental learning can also be performed in an offline setting, by providing the stored examples to the algorithm one by one, but without deleting them from memory. The term 'online' is often (inaccurately) used to denote 'incremental' [START_REF] Wilson | The general inefficiency of batch training for gradient descent learning[END_REF].

Model parameters and Meta-parameters

Figure 2 illustrates that regression algorithms assume a certain type of model. Each algorithm is designed to determine the optimal values of the parameters of this model, given an optimization criterion. Meta-parameters are algorithmic parameters that the user has to provide as an input to the algorithm. An understanding of linear least squares allows us to now give specific examples for an algorithm's model parameter and meta-parameters. This distinction is important for later sections.

The model used in linear least squares is a linear model, i.e. f (x) = a x. The model parameters are thus a. If we use an offset f (x) = a x + b, the model parameters are a and b. Model parameters are all the parameters that are required to make a prediction for a novel output.

Linear least squares has no meta-parameters; there is no parameter the user has to tune to make the algorithm work. But in linear least squares with Thikonov regularization, the user has to tune the parameter λ, which determines the trade-off between having small values in a and having small residuals. The parameter λ is thus a metaparameter of regularized linear least squares.

The important insight is that regularized lls and 'standard' lls use the exact same model, and thus have the same aim: optimize the values in the model parameter vector a. 'Standard' lls will always optimize a in the same way, and the values in a will always be the same for the same training data. With regularized lls however, the resulting values in a depend not only on the training data, but also on the meta-parameter λ.

The main contribution of this article is to show that all the algorithms in Table 1 use the same unified model, and they thus all have the same aim: optimize the values in the 'unified model parameter vector'. They algorithms only differ in the procedure used to determine the values of the parameters of the unified model, and in the meta-parameters that the user must tune.

From Linear Least Squares to Non-linear Regression

Linear least squares (lls) regression can be extended in two distinct ways to non-linear regression, yielding two classes of algorithms:

• Algorithms that perform multiple weighted lls regressions, using different input-dependent weighting functions. The resulting model is a mixture of linear models. This model is used in lwr, rfwr, lwpr, xcsf, gmr, and M5 (see Table 1 for full names).

• Algorithms that project the input space into a feature space using a set of non-linear basis functions, and performing one lls regression in this projected feature space. The resulting model is a weighted sum of basis functions. This model is used in rbfns, krr, gpr, svr, irfrls, i-ssgpr, CART, and elm.

In this section, we first present two algorithms -Locally Weighted Regression (lwr) and regression with Radial Basis Function Networks (rbfn) -which are representative for each class. We have chosen these two algorithms as representatives of their class for pedagogical reasons, because they are the most straight-forward extensions to linear least squares, and thus require only little explanation beyond the previous section. Based on the models used in these illustratory algorithms, we then introduce the unified model. In Section 4 and 5, we go through the algorithms in Table 1 one by one, we highlight their relationships and demonstrate how they all use the same unified model to represent the parameterized function.

Locally Weighted Regression (lwr)

When performing the fitting in lls, we may prefer to give certain data points more weight than others. Thus, for each training example (x n , y n ), we define a corresponding weight w n . With these weights, the sum of residuals is now computed as

S(a) = N n=1 w n (y n -a x n ) 2 . ( 26
)
When minimizing S(a) with this weighting, it is more important to minimize residuals that have large weights. For instance, examples (x n , y n ) that have a weight of w n = 0 do not contribute at all. Defining W to be an N × N diagonal weight matrix with W nn = w n , and rewriting (26) in matrix form

S(a) = (y -Xa) W(y -Xa), (27)
the weighted linear least squares solution becomes

a * = (X WX) -1 X Wy. (28)
The weights for each sample are typically defined as a function of the input space through a function φ, parameterized with θ

w n = φ(x n , θ). (29)
A commonly used weighting function is the multivariate Gaussian:

φ(x n , θ) = g(x n , c, Σ) with θ = (c, Σ) (30) g(x, c, Σ) = exp -1 2 (x -c) Σ -1 (x -c) , ( 31
)
where θ is a tuple containing the parameters of a multivariate Gaussian. An example of a Gaussian weighting function with c = π and Σ = 1.0 is plotted in graph B2 of Figure 7. In B1 , the weights w n are visualized by the thickness of the vertical lines representing the residuals. When comparing A1 and B1 , it becomes clear that weighting the examples can lead to a very different parameterization of the linear model. Again, we draw explicit attention to the fact that the model in linear least squares and weighted linear least squares is exactly the same (1); only the algorithm used to determine the model parameters a is different, i.e. compare (11) and (28).

Locally Weighted Regression (lwr) is an extension of weighted linear least squares, in which E independent weighted regressions are performed on the same data (in the design matrix X), but with E different weight matrices W e .

for e = 1 . . . E

a e = (X W e X) -1 X W e y. (32)
In practice, this usually means using the same type of weighting function, but with different parameter vectors θ e=1...E , for instance a set of E Gaussians with different centers c e=1...E , as depicted in graph Each data point in graph D3 is acquired by computing the value of both φ(x, c 1 ) and φ(x, c 2 ) (which can be read from D2), and plotting y against these two values. Because the 2D feature space φ(x, c 1 ), φ(x, c 2 ) is of a higher dimension that the 1D input space x, we now have a 3D plot.

The resulting model is a weighted sum of these linear models, where the weights are determined by the functions φ(x, θ e ):

f (x) = E e=1 φ(x, θ e ) • (a e x). (33)
If we include the offsets b e in the linear model and use the extended design matrix as explained in Section 2.3, we acquire the model

f (x) = E e=1 φ(x, θ e ) • (a e x + b e ), (34)
which is visualized in Figure 8.

x 1 x D . . . 1 φ 1 φ 2 φ E . . . y b 1 b2 bE a1, 1 a 1, 2 a 1 , E a D , 1 aD ,2 a D ,E w 1 w 2 w E
Figure 8: Function model in Locally Weighted Regression, represented as a feedforward neural network. The functions φe(x) generate the weights we from the hidden nodes -which contain linear sub-models (a e x + be) -to the output node.

Here, φe is an abbreviation of φ(x, θe)

Because each of the E local least squares regressions is performed independently of the others, the summation in ( 33) and ( 34) requires the weighting functions to be normalized, i.e. E e=1 φ(x, θ e ) = 1. For instance, the normalized Gaussian weighting function is:

φ(x, θ e ) = g(x, c e , Σ e ) E e =1 g(x, c e , Σ e )
with θ e = (c e , Σ e ).

(

) 35
In summary, the lwr algorithm performs several independent weighted linear least squares regressions on the same data, but each with a different weighting function, usually localized around a different part of the input space. The result is a weighted sum of linear models, where the models are linear, and where the weights are determined by the (normalized) weighting functions.

Regression with Radial Basis Function Networks (rbfns)

Instead of differently weighting each input sample -as in weighted linear least squares and lwran alternative extension to lls is to project each input sample into a higher-dimensional feature space using a set of basis functions. Thus, instead of a linear function model f

(x) = a x = D d=1 a d • x d we now have f (x) = E e=1 a e • φ(x, θ e ). ( 36
)
Here, a is no longer of size D×1 (with D the dimensionality of the input examples) but rather E × 1 (with E the number of basis functions). In the literature, it is more customary to use w ('weights') instead of a; a custom which we adopt in this section.

f (x) = E e=1 w e • φ(x, θ e ) ( 37
)
This model is a weighted sum of basis functions.

The term basis function is used because the subfunctions φ(x, θ e ) are the 'building blocks', or 'basis', from which f (x) is constructed.

When the basis functions are radial, it is called a Radial Basis Function Network (rbfn). This means that the output of the basis function depends only on the distance to a center, i.e. φ(||x -c||). A frequently used radial basis function is a Gaussian function, as in (31). Two Gaussian basis functions with centers c 1 and c 2 are plotted in D2 in Figure 7.

In a linear model f (x) = a x, the model is linear in the parameters a and linear in the input x. In contrast, ( 37) is non-linear in the input x, but it is still linear in the parameters w. Thus, we can readily apply lls to acquire the parameters w, even if the resulting function f (x) is not linear.

To apply lls to optimize w, we define the projected version of the design matrix -the feature matrix -which is N × E (one column for each basis function feature) instead of N × D as in ( 7)

Θ =      φ(x1, θ1) φ(x1, θ2) • • • φ(x1, θE) φ(x2, θ1) φ(x2, θ2) • • • φ(x2, θE) . . . . . . . . . . . . φ(xN , θ1) φ(xN , θ2) • • • φ(xN , θE)      (38)
and we get the least squares solution

w * = (Θ Θ) -1 Θ y. ( 39
)
This transforms non-linear regression in the input space into linear regression in a feature space, where the basis functions are responsible for projecting the input space into the feature space. The feature space is generally of higher dimensionality than the input space (i.e. E > D). In graph D3 of Figure 7, we visualize the projection of the data from the 1D input space into a 2D feature space using the two Gaussian basis functions. The plane which minimizes the residuals is also shown. The slopes w of this plane correspond to the weights of the basis functions, and w 1 φ(x, c 1 ) and w 2 φ(x, c 2 ) are plotted in D1 .

Note that the function can be approximated relatively well with only these two basis functions; this is because their position has been manually tuned for this illustratory example. Several methods for automatically tuning the number of basis functions and their parameters are presented in Section 5.2.

In summary, non-linear regression can be achieved by using a set of basis functions to project the input space into a feature space, and to perform lls in the feature space. The result is a weighted sum of linear models, where the models are the basis functions, and where the weights w are determined by a single linear regression. An example of such a model is a radial basis function network.

A Unified Model

When listing both models -a mixture of linear models (34) and a weighted sum of basis functions (37) -together

f (x) = E e=1 φ(x, θ e )•(b e + a e x) (40)
f (x) = E e=1 φ(x, θ e )• w e , (41)
it becomes clear that ( 41) is a special case of ( 40)

with a e = 0 and b e ≡ w e . This is also illustrated in Figure 9.

Going from a mixture of linear models (top lwr network in Figure 9) to a sum of weighted basis functions (bottom rbfn network in Figure 9) leads to an interesting shift of interpretation. In an rbfn, the sub-models in the hidden nodes are basis functions, whereas in lwr, these sub-models are linear models. In an rbfn, the weights are the (degenerate) lines (b e=1...E ), whereas in lwr, these weights are determined by the basis functions. Thus, the roles of the function implementing the sub-models/weights are switched between the representations. Due to this required shift in interpretation, the relationship between the two models may not be immediately obvious.

A closer inspection of the differences between the models also reveals why multiple lls regressions are necessary in lwr, whereas only one lls regression is required in rbfn. In (40), the model is not linear in the parameters b e and a e , and lls cannot be applied. In summary, the model used in lwr is a mixture of linear models, and in rbfn it is is a weighted sum of basis functions, the latter being a special case of the former with a = 0. In the next two sections, we list several algorithms that yield a mixture of linear models (Section 4) and a weighted sum of basis functions (Section 5). Throughout, it is important to keep in mind that the model of all these algorithms are special cases of the unified model in (40).

b2 bE a1,1 a 1,2 a 1 ,E a D ,1 aD, 2 a D, E w 1 w2 w E f (x) = E e=1 φ(x, θe)(a e x + be) | a = 0 ↓ x1 xD . . . 1 φ1 φ2 φE . . . y b 1 b2 bE w 1 w2 w E f (x) = E e=1 φ(x,

Model: Mixture of Linear Models

In Section 3.1 we summarized that lwr yields a function whose underlying model is a mixture of linear models, where the sub-models are linear, and where the weights are determined by the (normalized) weighting functions. In this section, we describe several algorithms which yield the same model, including Receptive Field Weighted Regression (rfwr), Locally Weighted Projection Regression (lwpr), xcsf, Gaussian Mixture Regression (gmr), and M5. Before turning to other algorithms, we briefly revisit lwr to discuss its algorithmic metaparameters. In the lwr algorithm [START_REF] Atkeson | Memory-based neural networks for robot learning[END_REF], the number of linear models, and the position and center of the weighting functions are fixed and must be specified by the user. Assuming diagonal covariance matrices for Gaussian weighting functions, this leads to 2 • E • D parameters to be tuned, i.e. a D-dimensional center and covariance matrix diagonal for each of the E weighting functions.

When the dimensionality of the input space is not too large, for instance when using lwr in the context of dynamical movement primitives [START_REF] Ijspeert | Dynamical movement primitives: learning attractor models for motor behaviors[END_REF], it may suffice to specify the number of weighting functions per dimension (i.e. D integers), and space the centers equidistantly within the bounds of the input data. The diagonals of the covariance matrices can be determined automatically by specifying the height at which the (unnormalized) weighting functions should intersect, see Figure 10 for an illustration. This reduces the number of meta-parameters for lwr to D+1. Alternatively, cross-validation over the training set may be used to determine these parameters. Although proposed almost two decades ago, lwr is still a competitive, widely-used algorithm, due to its simplicity. Fitting results may, however, become brittle if the input space dimensionality is high and not enough data is available.

Algorithm: Receptive Field Weighted Regression

The rfwr algorithm is the incremental variant of lwr, and it automates the choice of several model parameters [START_REF] Schaal | Receptive field weighted regression[END_REF]. The main differences to lwr are:

• rfwr is incremental, rather than batch, and uses Recursive Least Squares (rls) instead of a batch ls algorithm.

• New linear models and weighting functionscalled receptive fields -are added automatically by the algorithm when needed, i.e. when some input point is not covered by any receptive field.

• The centers and covariance matrices of the receptive fields are also adapted.

This flexibility comes at the cost of having to tune several meta-parameters, such as thresholds for removing or adding receptive fields as well as first and/or second order learning rates. These metaparameters control the dynamics of incremental learning; these dynamics are often difficult to predict, let alone the influence the meta-parameters have on the dynamics. Therefore, setting these meta-parameters is not always easy nor intuitive.

Algorithm: Locally Weighted Projection Regression

The lwpr algorithm [START_REF] Vijayakumar | Locally weighted projection regression: An O(n) algorithm for incremental real time learning in high dimensional space[END_REF] is an extension of rfwr, in which a lowdimensional affine projection is realized before fitting the linear models. For 1-dimensional input spaces, no such low-dimensional projection is possible, and lwpr becomes equivalent to rfwr (Vijayakumar and [START_REF] Vijayakumar | Locally weighted projection regression: An O(n) algorithm for incremental real time learning in high dimensional space[END_REF]. In lwpr, linear models are fitted using an iterative version of Partial Least Squares called nipals [START_REF] Geladi | Partial least squares regression: A tutorial[END_REF].

As outlined by [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF], lwpr is particularly interesting to perform regression when the data lies in a limited domain -because it adds receptive fields only in this domain -in a space with a high dimensionality -because it uses nipals to infer reduced linear models. Therefore, it has often been used to learn mechanical models of robots along trajectories, see [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF] for a survey. lwpr has many of the same meta-parameters as rfwr, and thus the same difficulties in tuning these meta-parameters apply.

Algorithm: XCSF

The xcsf algorithm [START_REF] Butz | Learning sensorimotor control structures with XCSF: redundancy exploitation and dynamic control[END_REF]) also uses a mixture of linear models. The algorithm manages a population of classifiers that contain a condition part and a prediction part. The condition part is usually defined as a Gaussian function which determines the area of influence of the prediction part. The prediction part is usually a linear model. The population of classifiers is optimized using a genetic algorithm, and the corresponding local linear models are learned using rls. Once the population of classifiers stabilizes, a condensation operator removes all the unnecessary classifiers.

A key feature of xcsf comes from the insight that the local linear models and corresponding Gaussian functions do not need to be defined over the full input space, an insight which can give rise to much more compact models, see [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF] for details.

Additionally, lwpr uses nipals instead of rls whereas xcsf can build local linear models in a predetermined subdomain of the input space. As lwpr, xcsf suffers from the necessity to empirically tune several meta-parameters, though some of them can be set more intuitively in xcsf than in lwpr. See [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF] for further discussion and Droniou et al. (2012a) for an empirical comparison.

Algorithm: Gaussian Mixture Regression

The underlying assumption in Gaussian Mixture Regression (gmr) is that the data in the joint input×output space can be well represented by a set of Gaussians, which is known as a Gaussian Mixture Model (gmm). For an illustration of a gmm consisting of three Gaussians, see Figure 11. Left: Gaussian Mixture Model (gmm) with 3 Gaussians on the 20 example data points. This model is acquired through unsupervised learning on the joint input×output space. Right: Gaussian Mixture Regression, which yield the mean and variance in the output space Y when conditioning on the input space X.

Training: unsupervised learning. A notable feature of gmr is that the training phase consists of unsupervised learning. It is performed by fitting a Gaussian Mixture Model (gmm) to the data with the Expectation-Maximization (em) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] algorithm. k-means clustering is commonly used to provide a first initialization of the centers of the Gaussians. Because em is an unsupervised learning algorithm, there is no distinction between an input x n and a target y n example, and they are concatenated into one vector z n = [x n y n ] . The gmm represents a model of the density of the vectors z n as a weighted sum of E Gaussian functions:

p(z n ) = E e=1 π e N (z n ; µ e , Σ e ), with E e=1 π e = 1. (42)
The em algorithm adjusts the priors π e and the parameters of the Gaussian functions (the means µ e and covariance matrices Σ e ) that define this model, see [START_REF] Ghahramani | Supervised learning from incomplete data via an EM approach[END_REF] for details. The only meta-parameter of the training phase is E, the number of Gaussians.

Regression: conditioning on the input. As noted by [START_REF] Ghahramani | Supervised learning from incomplete data via an EM approach[END_REF], the learned density can be exploited in several ways, since we can estimate y = f (x), x = f -1 (y), or any other relation between two subsets of the elements of the concatenated vector [x n y n ] .

In regression, we are interested in predicting y = E(y|x), the expectation of y given x. To do so, µ e and Σ e can be separated in input and output components as follows

µ e = [µ e,X , µ e,Y ] and Σ e = Σ e,X Σ e,XY Σ e,Y X Σ e,Y , (43)
For instance, if the dimensionality of the input space X is 2, and that of the output space Y is 1, then Σ e,X is a 2 × 2 sub-matrix, and Σ e,Y would be a scalar. The matrix Σ e would then be of size 3 × 3, being the covariance matrix of the full input × output space. Given the decomposition in (43), the expected output y given an input x is then computed as

y = E e=1 h e (x)(µ e,Y + Σ e,Y X Σ -1 e,X (x -µ e,X )), (44)
with:

h e (x) = π e N (x; µ e,X , Σ e,X ) E l=1 π l N (x; µ l,X , Σ l,X ) . ( 45
)
As a Bayesian method, an advantage of gmr over other lwr methods is that the variance can also be computed in the estimate of y with, see [START_REF] Calinon | On learning, representing, and generalizing a task in a humanoid robot[END_REF]:

var(y) = E e=1 h e (x) 2 Σ e,Y -Σ e,Y X Σ -1 e,X Σ e,Y X . ( 46
)
Relationship to the unified model. When we rewrite (44) for the univariate case and set a e = Σ e,Y X Σ -1 e,X and b e = µ e,Y -Σ e,Y X Σ -1 e,X µ e,X we get

y = E e=1 φ e (x)(a e x + b e ), ( 47
)
where Note that when representing gmr with the unified model, we can no longer compute the variance var(y) or the resulting confidence interval. This is because the matrices required to do so (Σ e,X , Σ e,X and Σ e,Y X ) have been compiled into a e and b e . Thus, for performing the regression from x to y the unified model is sufficient to capture the expected output of the gmm, but not to compute var(y).

Incremental gmr. [START_REF] Cederborg | Incremental local online Gaussian mixture regression for imitation learning of multiple tasks[END_REF] present a fast and incremental version of gmr named ilogmr. They add to gmr a dedicated structure for storing all the data points so as to quickly access the points that are close to a query point. The idea consists in performing local regression at the query point on demand, using only a small set of close points to build a small gmr model with only 2 or 3 Gaussian features. The resulting algorithm, ilo-gmr, is as accurate as gmr yet performs faster. e,Y and offsets be = µ e,Y -Σ e,Y X Σ -1 e,X µ e,X .

Algorithm: M5 (Model Tree Regression)

Model trees are decision trees which do not have class labels in their leaves, but rather linear models [START_REF] Quinlan | Learning with continuous classes[END_REF][START_REF] Wang | Induction of model trees for predicting continuous classes[END_REF]. The classical algorithm for growing model trees is M5 [START_REF] Quinlan | Learning with continuous classes[END_REF]. In each input dimension, it defines a decision boundary between each pair of adjacent input examples in that dimension. If there are N data points and the input dimensionality is D, there are D(N -1) decision boundaries. Each decision boundary splits the original set S into two subsets T 1 and T 2 . From these boundaries, the boundary that minimizes the standard deviation reduction is selected with

SDR = std(S) - i |T i | |S| × std(T i ), (49)
where std denotes the standard deviation and |...| denotes the cardinal of a set. Equation ( 49) is similar to the information gain criterion used to grow decision trees. This splitting procedure continues recursively until a given halting condition. After growing the tree, it may optionally be pruned to reduce its size. Each of E leaves in the tree is associated with a different disjoint subset T e=1...E of the original training set S, i.e. S = E e=1 T e and ∀eT e ⊂ S.

Within our unified model, the weighting functions of a model tree are not Gaussian, but rather box-car functions, which are 1 if the condition in the leaf holds, and 0 otherwise. A visualization of the model tree function representation for the example data is shown in Figure 13. As Figure 13 illustrates, there may be discontinuities when going from one linear model to the next. This is a consequence of using discontinuous weighting functions, i.e. the box-car functions. To avoid discontinuities, a smoothing procedure is usually applied when querying a model tree [START_REF] Quinlan | Learning with continuous classes[END_REF][START_REF] Wang | Induction of model trees for predicting continuous classes[END_REF].

One of the main advantages of a model tree is that its function representation is a tree, which can be easily visualized and interpreted. But this advantage holds mainly with trees of limited size, and with a large input space with many irrelevant features [START_REF] Stulp | Refining the execution of abstract actions with learned action models[END_REF].

The main meta-parameters for growing model trees with M5 are the minimum number of training examples per leaf (n), a splitting threshold on the ratio of the standard deviation of the entire data set and the standard deviation of the leaf under consideration (σ), and whether the tree should be pruned or not (P ). Finally, a smoothing parameter k determines how much the output is smoothed when querying the model tree.

Summary

A mixture of linear models is used in lwr, rfwr, lwpr, xcsf, gmr, and M5. Thus these algorithms only differ in the way they tune the weighting functions and the linear models, and in the meta-parameters the user has to provide as an input to the algorithm.

In this family, lwr is at one extreme, as all the parameters of the weighting functions are predetermined by the user and the local linear models are learned with a one batch least squares method, whereas in lwpr and xcsf the number of linear models and their weighting functions are adaptive but constrained through algorithmic metaparameters.

Model: Basis Function Network

In Section 3, we described two directions in which standard least squares has been extended. One extension leads to a model that consists of a mixture of linear models; algorithms that yield this model, lwr amongst others, have been discussed in the previous section.

In this section we turn to algorithms that yield a weighted sum of basis functions. This model is a special case of a mixture of linear models with a = 0, as demonstrated in Section 3.3 and repeated below in (50).

f (x) = E e=1 φ(x, θ e )• w e , (50)
In Radial Basis Function Networks (rbfns), presented in Section 3.2, the basis function φ(x, θ e ) is constrained to be radial, i.e. the output of φ depends only on the distance to a center: φ(||x -c||). We refer to the most generic case -when there are no constraints on φ -as a Basis Function Network (bfn). The different models presented in this section are all bfns, but with different basis functions: rbfn uses radial basis functions, neural networks use mostly sigmoid or hyperbolic functions, irfrls uses cosine functions, krls and gpr use kernel functions, and regression trees use box-car functions.

Because the basis function weights in all these regression methods can be found through least squares, in Section 5.1 we first describe three types of design matrices used in the least squares solution.

Treating them in a separate section avoids repeating these formulae for each of the algorithms; they are essentially the same for all.

Least Squares for Basis Function Networks

For all algorithms in this section, least squares regression can be used to learn the weight associated with each basis function. The standard and regularized solutions to the least square problem from Section 2 are repeated in ( 51) and ( 52):

w * = (Z Z) -1 Z y (51)
w * = (λI + Z Z) -1 Z y. (regularized) (52)

We now describe three forms that Z may take: the design matrix, the feature matrix, or the Gram matrix.

Design Matrix

In linear least squares, the matrix Z is of size N × D, i.e. each row contains one input example, with each column representing a different dimension of the input space. Informally, Z may be considered as containing the 'raw' input data.

Z = X =      x 1,1 x 1,2 • • • x 1,D x 2,1 x 2,2 • • • x 2,D . . . . . . . . . . . . x N,1 x N,2 • • • x N,D .      (53)

Feature Matrix

As discussed in Section 3.2, projecting the input space into a higher dimensional feature space with E basis functions, yields an N ×E feature matrix (one column for each basis function feature) to which Equations ( 51) or ( 52) are readily applicable.

Z = Θ(X) =      φ(x 1 , θ 1 ) φ(x 1 , θ 2 ) • • • φ(x 1 , θ E ) φ(x 2 , θ 1 ) φ(x 2 , θ 2 ) • • • φ(x 2 , θ E ) . . . . . . . . . . . . φ(x N , θ 1 ) φ(x N , θ 2 ) • • • φ(x N , θ E )      (54)

Gram Matrix (for Kernel Functions)

A kernel is a special type of basis function which relates two elements of the input space to each other, see [START_REF] Williams | Gaussian processes for machine learning[END_REF] for their mathematical properties. In this case, Z ≡ K(X, X) is of size N × N , because the kernel is evaluated on all combinations of input examples. This matrix is known as the Gram matrix.

Z = K(X, X) =      k(x 1 , x 1 ) k(x 1 , x 2 ) • • • k(x 1 , x N ) k(x 2 , x 1 ) k(x 2 , x 2 ) • • • k(x 2 , x N ) . . . . . . . . . . . . k(x N , x 1 ) k(x N , x 2 ) • • • k(x N , x N )      . ( 55
)
Note that the same kernel function k is used for all entries of the Gram matrix. For convenience, we therefore drop the kernel parameters θ.

Because the Gram matrix is symmetrical, (51) can be simplified as

w * = (K K) -1 K y (56) = K -1 y, (57)
where we use K as an abbreviation for K(X, X).

The same is possible for the regularized case, if we replace the general Z matrix with the kernel matrix K in (52), after some algebra we get

w * = (λI + K) -1 y. ( 58
)
For more details about this computation, and in particular for a discussion about using the kernel trick [START_REF] Saunders | Ridge regression learning algorithm in dual variables[END_REF] in this context, see [START_REF] Williams | Gaussian processes for machine learning[END_REF].

Preview

Table 2 provides a preview of the algorithms presented in this section, listing which type of basis functions they use, the type of matrix used for linear least squares, and the least squares variant they use for batch regression. This table does not list incremental methods for learning the model parameters.

Algorithm

Features?

Size of Z? Regularized?

rbfn RBFs FM (E × N ) Yes krr kernels Gram (N × N ) Yes gpr kernels Gram (N × N ) No irfrls cosine FM (E × N ) Yes i-ssgpr cosine FM (E × N ) Yes elm sigmoid FM (E × N ) No
Table 2: Design of all weighted basis function algorithms (FM: feature matrix, RBFs: radial basis functions).

Meta-parameters and model complexity. The table above already gives some insight into the types of meta-parameters these algorithms have. For instance, in krr and gpr each data point is the center of one kernel. Thus, these centers depend on the training data alone, and must not be set as metaparameters. When applying linear least squares to rbfns or in irfrls, i-ssgpr or elm on the other hand, the placement of basis functions is arbitrary, and must thus be chosen in advance.

There are several strategies for choosing the centers of the basis functions: 1) setting them manually; 2) equidistantly spacing in the range of the input instances along each dimension; 3) randomly sampling from the input range. 4) randomly sampling from the input instances; 5) clustering the input examples and using the cluster centers as basis function centers; In the first strategy, the centers c e=1...E are part of the meta-parameters, and must be set by the user. In the other cases, only E is a meta-parameter, and the centers are determined automatically. Similar strategies exist for cosine basis functions (irfrls and i-ssgpr choose the phases randomly) and sigmoid basis functions (elm chooses the bias randomly).

The choice of the number of basis function E affects the complexity of the model, and thus bias/variance trade-off. krr and gpr take an extreme position on model complexity, by using all the training data as part of the model, thus letting the "data speak for itself". For this reason, such algorithms are known as non-parametric regression methods [START_REF] Williams | Gaussian processes for machine learning[END_REF]. Although this article focuses on parametric regression, we have nevertheless included krr and gpr because of their similarity to rbfns and other basis function networks presented in this section. The disadvantage of this approach in gpr and krr is that letting each data point speak requires patience when there is a lot of data; it involves the inversion of a N × N matrix, which has complexity O(N 3 ).

Algorithm: Regression with Radial Basis

Function Networks (revisited)

Before turning to other algorithms, we briefly revisit regression with rbfns to discuss its algorithmic meta-parameters. In rbfns, the basis functions φ are chosen to be symmetric around some center c e , which is therefore always a parameter of φ, i.e. φ(x, c e , . . . ). Commonly used radial basis functions (rbfs) include the Gaussian (31). Different strategies for choosing these centers were presented in the previous section. Each of these versions yields a slightly different algorithm, and they do not all have the same meta-parameters. And although these algorithmic variations all use the exact same model, the algorithms may tune the values of the model parameters in different ways, as discussed in Section 2.6.

For the relationship between between rbfns, regularization networks, generalized splines, and di-mensionality reduction, we refer to the work by [START_REF] Poggio | Networks for approximation and learning[END_REF].

Algorithm: Kernel Ridge Regression

In Kernel Ridge Regression (krr), also called Kernel Regularized Least Squares, the basis functions φ are generated from a kernel function k(x, x ), which takes two vectors from the input space as input. Kernel functions are such that their output is maximal when x = x and decreases as the distance ||x -x || increases, which provides a locality property. An example of such a function is the Gaussian kernel.

Using these kernel functions, regression consists in finding the weights w of the function

f (x) = N n=1 w n • k(x, x n ). ( 59
)
To perform the regression with this approach, we first compute the Gram matrix in (55), and then perform regularized ridge regression with (58). How does (59) relate to the unified model in (41)? krr turns out to be a special case where the number of basis function is identical to the number of data points, and where each data point is the center of a basis function. The corresponding network model is depicted in Figure 14. Thus, while rbfn allows arbitrarily placed basis function centers, krr centers one basis function around each data point.

Algorithm: Gaussian Process Regression

In Gaussian Process Regression (gpr), the key idea is that the output data y = {y 1 , . . . , y n }, or any subset of the output data, can be thought of as one sample from a multivariate (n-variate) Gaussian function [START_REF] Williams | Gaussian processes for machine learning[END_REF][START_REF] Ebden | Gaussian processes for regression: A quick introduction[END_REF]. Often, it is assumed that the mean of this distribution is 0, i.e. the Gaussian Process is represented as y ∼ N (0, Σ).

(60)

The covariance matrix Σ is determined by the covariance function or kernel. A typical choice is the Gaussian function

k(x, x ) = σ 2 f exp -1 2 (x -x ) W -1 (x -x ) (61) = σ 2 f g(x, x , W), (62)
as for smooth functions, we expect in general higher covariance between points that are closer. Here, σ f represents signal variance, and W the covariance matrix of the Gaussian kernel. The Gaussian function g, was previously defined in (31). For a given covariance function k and N training points {(x n , y n )} N n=1 the corresponding Gaussian Process is

y ∼ N (0, K(X, X)), ( 63
)
where K is again the Gram matrix, as in (55).

Predicting y q for a novel input x q is done by assuming that the novel output y q is also sampled from a multivariate Gaussian with

[ y yq ] ∼ N (0, K(X,X) k(xq,X) k(xq,X) k(xq,xq)
), and (64) k(x q , X) = [k(x q , x 1 ), . . . , k(x q , x n )].

(65)

Conditioning this joint distribution to predict y q given X, y, and x q yields another multivariate Gaussian [START_REF] Ebden | Gaussian processes for regression: A quick introduction[END_REF], which predicts the probability of observing a certain y q :

yq|X, y, xq ∼ N (k(xq, X)K(X, X) -1 y, k(xq, xq) -k(xq, X)K(X, X) -1 k(xq, X) ). ( 66
)
Thus, the best estimate for y q is the mean, and the uncertainty in y q is captured by the variance as

y q = k(xq, X)K(X, X) -1 y (67) var(yq) = k(xq, xq) -k(xq, X)K(X, X) -1 k(xq, X) . ( 68
)
Figure 15 visualizes the mean y q and variance var(y q ) of the Gaussian process over the input domain, given 10 and 20 data points of our example training data. We see that for areas where no training data is available, the mean tends to go towards zero (i.e. to the mean of the Gaussian Process) and The open parameters of a Gaussian Process, e.g. the parameters σ 2 f and W in ( 62), are known as the hyperparameters. In our terminology, these correspond to the meta-parameters of the gpr algorithm for regression. Rather than tuning the hyperparameters by hand, it is customary to tune them automatically by minimizing their log likelihood on a training dataset, i.e. argmin θ log p(y|X, θ), where θ contains all the hyperparameters. For a zero-mean Gaussian process with Gaussian kernel as in (61) for instance, θ contains σ 2 f and W. In this case, the only meta-parameters of the gpr algorithm are those related to the optimization of the hyperparameters of the Gaussian process.

Relation to the Unified Model. When computing the mean y q , K(X, X) and y depend only on the training data, not the novel input x q . Therefore, K(X, X) -1 y can be compacted in one weight vector, which does not depend on the query x q . This substitution of K(X, X) -1 y with w * corresponds exactly to the linear least squares solution in (57). The substitution yields

y q = k(x q , X)K(X, X) -1 y (69) = k(x q , X) • w * apply (57) (70) = [k(x q , x 1 ) . . . k(x q , x N )] • w * (71) = N n=1 w n • k(x q , x n ). ( 72
)
Thus, we see that the mean of the gpr is the same weighted sum of basis functions as in krr (see ( 59)), and the same network representation as in Figure 14 applies. As krr, gpr may thus be considered as a special case of rbfn, where each data point is the center of one basis function, and each basis function has the same width, determined by the covariance function. Such links have been previously made by [START_REF] Williams | Gaussian processes for machine learning[END_REF]; [START_REF] Williams | Computation with infinite neural networks[END_REF]; [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF]. Therefore, gpr is also a special case of the unified model in (41). From this perspective, the only difference between krr and gpr is that the former uses regularized least squares, and the latter standard least squares (with λ = 0). Figure 16 is a repetition of Figure 15, except that we now visualize the basis functions (bottom row) and weighted basis functions (top row), just as was done for rbfn in Figure 7 D1 and D2 . As gpr is a Bayesian method, it provides an estimate of the variance var(y q ). However, representing the gpr model with the unified model means that we can no longer compute the variance var(y). This is because the matrices required to do so K(X, X) and y have been reduced to the vector w * . Thus, for performing the regression from x to y the unified model is sufficient to capture the gpr model, but to compute var(y) or the confidence interval, it is not. The same holds for gmr, see Section 4.5.

Non-zero mean function. Prior knowledge about a function may be introduced in the mean of the Gaussian process. For instance, if we know the average outside temperature at a specific location is 20 • C, we may rather set the mean of the Gaussian process to µ = 20. The mean itself may depend on the input, so we write µ(x), which leads to a Gaussian Process:

[ y yq ] ∼ N ( µ(X) µ(xq) , K(X,X) k(xq,X) k(xq,X) k(xq,xq) ), (73)
which, when conditioned on X, y, and x q yields (the variance estimate stays the same)

y q = µ(x q ) + k(x q , X)K(X, X) -1 (y -µ(X)). ( 74
)
Note the similarity for computing the most likely estimate of in Gaussian mixture regression in ( 44).

In the unified model, adding a non-zero mean to the Gaussian process thus means the weights w * become K(X, X) -1 (y -µ(X)), and a bias µ(x q ) needs to be added.

Sparse and online versions of gpr. Given the cubic computation cost in the number of examples n due to the inversion of K(X, X), gpr is often not applicable online. There are two approaches to avoiding the full inversion of K(X, X). The sparse approach consists in limiting n by adequately choosing the points to remember and forgetting the rest. A good overview of this approach can be found in Quiñonero [START_REF] Quiñonero Candela | A unifying view of sparse approximate Gaussian process regression[END_REF]. For instance, it has been recently used in Sparse Online Gaussian Processes (sogp), in the context of learning control policies from demonstrations [START_REF] Grollman | Sparse incremental learning for interactive robot control policy estimation[END_REF].

The mixture of experts approach rather consists in splitting the global domain of the function into smaller regions, as done in the lwr family of methods. Such combination of ideas from both families are taken in Local Gaussian Process Regression (Nguyen-Tuong et al., 2009), see [START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: a survey[END_REF] for details, and more recently in [START_REF] Meier | Local gaussian regression[END_REF].

Algorithm: Support Vector Regression (svr)

The L 1 and L 2 loss functions penalize any residual that is larger than 0. In Support Vector Regression (svr) a more tolerant approach is taken. Any residual r n = y n -f (x n ) that is smaller than a given threshold |r n | < is not penalized. The region in which the residuals are zero, visualized in Figure 17, is called the -insensitive tube. This leads to the -insensitive loss function

L (r) = 0 if r ≤ |r| - if r > ( 75
)
Figure 17: Support Vector Regression. Data points within the -insensitive tube have zero penalty, and are gray.

Those outside are the support vectors, which are black. Adapted from [START_REF] Smola | A tutorial on support vector regression[END_REF].

This loss function cannot be minimized analytically, and the problem is rather stated as a constrained optimization problem [START_REF] Smola | A tutorial on support vector regression[END_REF])

minimize 1 2 a 2 + C N n=1 (ξ + n + ξ - n ), (76)
subject to

     r n ≤ + ξ + n -r n ≤ + ξ - n ξ + n , ξ - n ≥ 0 (77)
which uses the slack variables3 ξ + n and ξ - n , which represent the distance above or below theinsensitive tube, cf. Figure 17. Informally, this problem can be formulated as "Try to find small slopes a, such that the distances above ξ + n or below ξ - n the tube are minimal; ideally points lie within the tube, because this leads to ξ + n = 0 and ξ - n = 0, and thus zero penalty for that point." C plays the same role as λ in Thikonov regularization; it determines the trade-off between minimizing the residuals and minimizing the norm of a.

This constrained optimization problem is solved by defining a Lagrange function that incorporates both the term to be minimized and the constraints 3 . The final result for linear svr, in which the under-lying model is f (x) = a x + b, is:

a = N n=1 (α + n -α - n )x n , (78)
where α + n and α - n are the Lagrangian multipliers in the Lagrange function for the first two constraints in (77). With this solution, the linear function becomes:

f (x) = a x + b (79) = N n=1 (α + n -α - n )x n x + b (80) = N n=1 (α + n -α - n ) x n , x + b. (81)
The slopes a are thus represented4 as a linear combination of the training examples x n=1...N .

Note that the data points inside the -insensitive tube do not contribute to the computation of (80), because α + n and α - n are 0 for these points, and (0 -0)x n will thus always be 0, independent of the values in the vector x n . Since only the examples outside the tube contribute to (81), i.e. only the example input vectors support it, they are known as support vectors.

Using the L 2 norm in krr and gpr leads to a dense representation, in which all example data points are required to make a prediction, and must be stored in the model parameters. Using the L norm leads to a sparse representation, in which only a limited number of data points -the support vectors -must be stored. Thus, only the support vectors are part of the model parameters.

Linear least squares is readily extended to nonlinear regression by projecting the input space into a feature space, as explained in Section 3.2. When applying this approach to svr by using a basis function φ k=1...K (x) to project in the feature space (and using w rather than a to denote weights instead of slopes), the solution becomes

w k = N n=1 (α + n -α - n )φ k (x n ), (82)
f (x) = K k=1 N n=1 (α + n -α - n )φ k (x n ) φ k (x) + b (83) = N n=1 (α + n -α - n ) φ(x n ), φ(x) + b (84) = N n=1 (α + n -α - n )k(x n , x) + b. (85)
This last step is the "kernel trick", which is based on the following insights: 1) svr only requires dot products in the feature space φ(x n ), φ(x) .

2) certain types of kernels k(x, x ) -see Smola and Schölkopf (2004, Section 2.3) for their properties -correspond to a dot product in a feature space. 3) φ(x n ), φ(x) may be replaced with k(x, x ), which avoids computing the potentially high-dimensional feature vector φ(x) explicitly. The representation in ( 85) is compatible with the unified model in ( 37) with E = N (every data point x n is the center of a basis function φ n≡e ) and w e = (α + e -α - e ). The offset b can be pulled into the weight vector by defining w 0 = b and k 0 (.) = 1.

The meta-parameters of svr are the tolerance parameter and the trade-off parameter C. An additional parameter ν can be introduced to deal with the trade-off between the value of and model complexity, giving rise to the ν-svr algorithm [START_REF] Schölkopf | New support vector algorithms[END_REF].

Algorithm: irfrls

To circumvent the kernel expansion problem present in krr and gpr, [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] propose to approximate the Gram matrix K with a set of random features, giving rise to Random Features Regularized Least Squares (rfrls). [START_REF] Gijsberts | Incremental learning of robot dynamics using random features[END_REF] present an incremental version of rfrls, named irfrls hereafter and i-rfrr in [START_REF] Gijsberts | Real-time model learning using incremental sparse spectrum gaussian process regression[END_REF]. Though the authors present their algorithm starting from krr, from the perspective of this paper the relationship to rbfns is stronger. Indeed, with respect to rbfns, the key idea behind irfrls is that any function can be approximated arbitrarily well as a weighted sum of cosine functions, as outlined by the theory underlying the Fourier Transform (see Droniou et al. (2012b) for a presentation of irfrls based on its relationship to Fourier Transform). Thus, instead of using Gaussian kernels as in rbfns, f (x) is as a set of E cosine functions

f (x) = E e=1 w e • φ(x, ω e , ψ e ) (86) = E e=1 w e • cos(ω e x + ψ e ). ( 87
)
The basis functions are defined by randomly drawing their multidimensional frequency ω and phase ψ using ω ∼ N (0, 2γI) and ψ ∼ U(0, 2π).

As in all other rbfn-like methods, the main meta-parameter in irfrls is the number of features E, that is easy to tune considering a trade-off between the model accuracy and the time complexity. The other two meta-parameters of irfrls are λ, the regularization parameter and γ, the variance over the period of the random cosine features. Experiments on several databases confirm that rfrls and irfrls can be set close to krr with a small enough processing time. Indeed, irfrls shows interesting practical performances (Droniou et al., 2012a) and its computation time is independent of the learning data.

To summarize, irfrls is a basis function network approach that uses randomly sampled features, as elm and certain variants of rbfn do. But it uses cosine features instead of radial or sigmoid features.

Algorithm: I-SSGPR

In order to approximate a function, krr uses (59) whereas gpr uses (67). A quick way to describe issgpr is to say that it applies to gpr the same idea that irfrls applies to krr: instead of using one kernel per data point, it draws a collection of random cosine feature functions and then it solves the resulting -non-regularized -linear least square problem into the projected domain. Actually, issgpr also inherits from the gpr framework the capability to optimize most of meta-parameters using the marginal likelihood [START_REF] Williams | Gaussian processes for machine learning[END_REF], so the only meta-parameter that remains to be tuned is E.

Algorithm: CART (Regression Trees)

Regression trees [START_REF] Breiman | Classification and Regression Trees[END_REF] are a special case of model trees for which the slopes of the linear model are zero (a = 0), i.e. the leaf contains only one continuous value (b). Thus model trees and regression trees follow our classification between a mixture of linear models (model trees), and a weighted sum of basis functions (regression trees), the latter being a special case of the former with degenerate linear models with a = 0. model trees and regression trees are grown using very similar algorithms (M5 and "CART learning", respectively), and have essentially the same advantages and disadvantages, we refer to Section 4.6 for the full discussion. [START_REF] Orr | Combining regression trees and radial basis function networks[END_REF] use regression trees to initialize a radial basis function network, thus highlighting the similarity between these two representations.

Algorithm: Extreme Learning Machine (elm)

The Extreme Learning Machine (elm) is a batch regression algorithm for training the weights of a Single-Hidden Layer Feedforward Network (slfn) (Huang et al., 2006b[START_REF] Huang | Extreme learning machines: a survey[END_REF]). An slfn is an Artificial Neural Network (ann) that has one input layer, one hidden layer, and one output layer, as depicted in Figure 19. The representation 5 of an slfn is

f (x) = E e=1 w e ψ D d=1 v ed x d , (88)
= E e=1 w e ψ (z e ) , (89)
where z e = D d=1 v ed x d is the input to the activation function ψ, which is often chosen to be sigmoid (ψ(z e ) = (1+e -ze ) -1 ) or hyperbolic (ψ(z e ) = tanh(z e )). The model is easily extended to have multiple output nodes [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. The elm algorithm for determining the weights of the slfn is essentially the same as when using rbfns for regression, i.e. they both perform batch linear least squares on a projected feature space. In the elm approach, the basis function parameters v e=1...E are first set randomly, and are then used to generate the feature matrix Θ(X) (Huang et al., 2006b). Least squares (39) is then used to determine the parameters w.

y v11 v 21 v E 1 v 1 D v2D v E D w 1 w2 w E
Some extensions to elms include growing the number of hidden nodes incrementally (Huang et al., 2006a) -whereby hidden nodes may be preselected from a random set (Huang and Chen, 2008) -or extending them to the complex domain (Huang et al., 2008).

The slfn model ( 88) used in elms is readily made compatible with the unified model by pulling the weights from the input layer to the output layer into the activation function, i.e.

f (x) = E e=1 w e φ(x, v e ). ( 90
)
Thus, from the model-centric perspective taken in this paper, a slfn is a special case of a basis func-tion network6 (and thus the unified model), where the basis function must have the form

φ(x, v e ) = ψ D d=1 v ed x d = ψ ( v e , x ) , (91)
where v e , x denotes the inner product between v e and x.

Algorithm: Backpropagation

An incremental approach to optimizing the connection weights in an artificial neural network is the backward propagation of errors [START_REF] Werbos | Beyond regression: New tools for prediction and analysis in the behavioral sciences[END_REF], or simply backpropagation. This algorithm first propagates an input x n forwards through the network (from the input to the output layer), and computes the error 1 2 (y n -f (x n )) 2 at the output node. It then propagates this error backwards through the network, and updates the weights in each layer through gradient descent. Although backpropagation is usually associated with artificial neural networks that use sigmoid or hyperbolic activation functions, the algorithm is also applicable to radial basis function networks [START_REF] Poggio | Networks for approximation and learning[END_REF]Schwenker et al., 2001a).

In Figure 20, we depict one neuron of an ann, agnostic of whether it is in a hidden or an output layer. The local rule for updating one weight of one such neuron is:

w ji ← w ji + αδ j u i (92)
where δ j is the error between the actual output of the neuron and the target output, u i is the input from neuron i to neuron j, and α is the learning rate. This rule is acquired7 by computing the partial derivative of the error of the neuron w.r.t. the weight ∂E ∂wji . If the node is an output node, the error δ j depends on the residual r n = y n -f (x), i.e.

δ j = ψ j (z j )(y n -f (x n )). ( 93
)
If the node is a hidden node, δ j is defined in terms of errors δ k of neurons in the subsequent layer, i.e. This recursive rule for computing δ is what propagates errors from the back (where the errors correspond to the residuals at the output neuron) to the front of the network (the input layer).

δ j = ψ j (z j ) K k=1 δ k w kj . ( 94
) u 1 u 2 u D . . . z j = D i=1 w ji u i ψj(zj) o j w j 1 w j2 w j D
So far, we have presented backpropagation as an incremental algorithm, i.e. it is able to update the weights incrementally with (92), one training example at a time. However, a batch version exists also, which updates the weights at once, given all the available data:

w ji ← w ji + α N n=1 δ n,j u n,i (batch version). (95)
In mini-batch learning, the weights are updated with (95), but using sub-batches of size K < N [START_REF] Wilson | The general inefficiency of batch training for gradient descent learning[END_REF]. Incremental and standard batch training can thus be seen as special cases with K = 1 and K = N respectively. [START_REF] Wilson | The general inefficiency of batch training for gradient descent learning[END_REF] show that using batch or mini-batch for backpropagation on a speech recognition task provides no advantage over incremental learning in terms of learning accuracy or speed.

Backpropagation for slfns

In the slfn depicted in Figure 19, there is only one hidden layer, and the output node does not have an activation function (it is identity, and thus ψ (z e ) = 1). When applying the local update rule (92) to update the weight w e , we have δ = y n -f (x n ) from ( 93 (96)

If we keep the weights v e=1...E fixed, and only the weights to the output are thus to be trained, ( 96) is an incremental gradient-based alternative to the batch least squares used in elm.

For the weights from the input layer to the hidden layer, δ d = ψ (z d )δ e w e and u d = x d , which yields the update rule

v ed ← v ed + αψ (z d )δ e w e x d .
(97)

Comparison of elm and Backpropagation. [START_REF] Huang | Extreme learning machines: a survey[END_REF] show that training elms with random features and least squares is much faster than with backpropagation. It has lower computational cost than using gradient back-propagation, despite the cubic cost of computing the inverse (X X) -1 . [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF] analyze the advantages of using random features instead of backpropagation. They show that, although being much easier to obtain, random features can generate close to optimal solutions if numerous enough. This analysis is also at the heart of the design of irfrls and i-ssgpr (see Sections 5.6 and 5.7). A complementary analysis comes from [START_REF] Widrow | The no-prop algorithm: A new learning algorithm for multilayer neural networks[END_REF], who examine under which conditions training a network benefits from backpropagation or not. Schwenker et al. (2001a) use a combination of both least squares batch learning and incremental backpropagation.

References to empirical comparisons between elm, backpropagation, and other algorithms are given in Section 6.2.

Backpropagation for Multi-Layer Feedfor-

ward Networks (mlffs) The network visualized in Figure 19 is easily extended to incorporate more than one hidden layer. Following the notation by [START_REF] Poggio | Networks for approximation and learning[END_REF], the function representation of such multi-layer feedforward networks (mlffs) is

f (x) = E e=1 weψ D d=1 v ed ψ . . . ψ C c=1 u dc xc . . . . ( 98
)
An advantage of the backpropagation algorithm is that it optimizes all of the weights in the network from output to input with the generic update rule (92), and thus readily extends to multi-layer networks. In the multivariate, high-dimensional case, [START_REF] Friedman | Projection pursuit regression[END_REF] propose an advanced gradient descent algorithm for networks with several layers. In contrast, the least squares approach used in elms can only be used to tune the weights from the final hidden layer to the output layer, because these are the only weights which are linear with respect to the output.

Deep Neural Networks (DNN) [START_REF] Bengio | Learning deep architectures for AI[END_REF] are networks, typically feedforward, that have many hidden layers. This allows more abstract features to be represented in the hidden layers. Gradient backpropagation is also used for tuning the weights of the network, though a first stage of unsupervised learning is commonly used to preshape those features [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF][START_REF] Bengio | Representation learning: A review and new perspectives[END_REF], to avoid the vanishing gradient problem [START_REF] Hochreiter | Gradient flow in recurrent nets: the difficulty of learning long-term dependencies[END_REF]. DNNs are more commonly used for classification, where they outperform their shallow counterpart. Whether this superiority extends to regression is still an open question. Even for classification, understanding the reasons for the superiority of DNNs is an active topic of investigation. Researchers who address these questions tend to use elementary basis functions that are piecewise linear rather than nonlinear, because it provides a better grip on their composition [START_REF] Pascanu | On the number of inference regions of deep feed forward networks with piece-wise linear activations[END_REF]. These investigations are too recent to be properly covered in this review, so the reader is referred to [START_REF] Pascanu | On the number of inference regions of deep feed forward networks with piece-wise linear activations[END_REF]; [START_REF] Alain | What Regularized Auto-Encoders Learn from the Data Generating Distribution[END_REF]; [START_REF] Dauphin | Identifying and attacking the saddle point problem in high-dimensional non-convex optimization[END_REF] for more references to this quickly evolving domain.

Our unified model can in principle be used to represent a DNN, if we consider all the layers up to the final layer of a DNN to be the implementation of a rather complex basis function (consisting of weighted compositions of the basis functions in preceding layers). However, this shallow perspective does not capture the dynamics of learning the weights in preceding layers through backpropagation, nor does it provide insight in the relations to other types of regression algorithms.

Backpropagation for Recurrent Neural Networks (rnns)

In feedforward networks, only connections to nodes in subsequent layers are allowed. Recurrent neural networks do not have this constraint, and connections to nodes in previous layers or nodes within a layer are possible. This leads rnns to have an internal state, which changes over time even if the input does not change. Backpropagation through time is a generalization of backpropagation, enabling this gradient-based method to be applied to rnns [START_REF] Werbos | Generalization of backpropagation with application to a recurrent gas market model[END_REF]. Due to the cycles in rnns, and the internal state that arises from them, rnns cannot be captured by the unified model.

Comparison of rbfns and mlffs

Bishop (1995, Section 5.8) summarizes some of the differences between rbfns and mlffs. rbfns use radial basis functions, which take the distance to a center as an input. Each basis function is thus a localized representation, and only a few basis functions are active for a given input. In contrast, an mlff is a distributed representation, because the sigmoid function -which splits the input space into an active and inactive spaces -activates many of the hidden units. Furthermore, rbfns typically have only one layer of hidden units, whereas mlffs often have many hidden layers [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. Finally, rbfns weights are typically trained in two phases -randomly setting or preshaping features and then performing batch least squares -although backpropagation has been applied also (Schwenker et al., 2001a;[START_REF] Poggio | Networks for approximation and learning[END_REF]. On the other hand, mlffs are typically trained with backpropagation, because least squares can only be applied to the last layer of mlffs. [START_REF] Dorffner | A unified framework for MLPs and RBFNs: Introducing conic section function networks[END_REF] discusses these differences in greater details, before introducing a more general class of conic basis functions that encompasses radial basis functions and sigmoid activation functions. The paper also discusses the interest of combining both families of functions in the same approximator, as in [START_REF] Cohen | A hybrid projection-based and radial basis function architecture: initial values and global optimisation[END_REF].

Apart from standard activation functions such as Gaussian, hyperbolic or sigmoidal functions, many other families of basis functions can be used. One can even evolve the activation function of neurons with evolutionary techniques, sometimes leading to higher performance than sigmoidal activation functions [START_REF] Augusteijn | Evolving transfer functions for artificial neural networks[END_REF]. However, as noted by [START_REF] Dorffner | A unified framework for MLPs and RBFNs: Introducing conic section function networks[END_REF], finding the right class of activation functions might be very problem dependent or give rise to the curse of dimensionality. One can also use neural units whose activation function is of higher order than a sigmoidal function applied to a weighted sum of inputs. The representational power and learning complexity of several instances of this family of networks, such as product unit networks, is studied in [START_REF] Schmitt | On the complexity of computing and learning with multiplicative neural networks[END_REF].

Summary

All the algorithms listed in this section use the same function model

f (x) = E e=1 φ(x, θ e )• w e , (99)
which is a special case of the unified model in ( 41) with a = 0. However, different algorithms use different instantiations of (99), as listed in Table 2.

In krr and gpr, the number of basis functions is taken to be N whereas in all other algorithms this is an arbitrary parameter E. Additionally, in gpr and krr, the basis functions are centered on the data points whereas in rbfns, they can be either regularly placed to pave the input space or drawn randomly. In all models but anns, irfrls and i-ssgpr, it is standard to use Gaussian basis functions (that combine the kernel and radial basis function properties), whereas irfrls and i-ssgpr use cosine basis functions and anns more often use sigmoids or hyperbolic functions.

Some types of anns are special cases of the unified model (slfn) and some are not (rnns). The backpropagation algorithms is applicable to all anns, whereas least squares is only applicable to slfns. The latter case, known as the elm, is very similar to training rbfns. For further connections between these methods, the reader is referred to [START_REF] Smola | A tutorial on support vector regression[END_REF], who examine the relationship between svr and anns.

The Unified Model Revisited

All the algorithms discussed in this paper use the same generic function model in ( 41), but different algorithms use special cases of this generic model. In this section, we reconsider the relationships between the two main classes of models. Furthermore, we list, for the different algorithms, which model parameters must be pre-specified by the user, and which are determined by the algorithm itself. This model-based perspective on the algorithm leads to several observations, some straightforward, and some rather intriguing.

• Least squares is special case of weighted least squares, with w = 1

• Weighted least squares is a special case of Locally Weighted Regression with E = 1

• Gaussian process regression, which uses standard least squares, is the Bayesian interpretation of krr, which rather uses regularized least squares.

• The model used in Gaussian process regression is a radial basis function network with a potentially infinite number of hidden basis functions [START_REF] Williams | Computation with infinite neural networks[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks[END_REF].

• irfrls is a basis function network that uses (non-radial) randomly generated cosine basis function. Its excellent performance in practice raises interesting questions, because irfrls differs from all other algorithms in this paper but artificial neural networks in two important ways: it uses global (cosine) rather than local (radial) basis functions, and it generates their parameters randomly, rather than from the data.

• regression trees are to rbfns (both using a weighted sum of basis functions) what model trees are to lwr (both using a mixture of linear models), except that the trees use box-car basis functions. Thus, even though regression and model trees are learned using the same algorithm, the underlying model they learn belongs to a different class.

Meta-Parameters and Empirical Comparisons

Regression algorithms differ in the number of meta-parameters that need to be set. For instance, for lwr, the number of basis functions, as well as their centers and widths, must be specified by the user. For gmr, gpr and i-ssgpr, on the other hand, only the number of basis functions must be set. For completeness, Table 3 summarizes the meta-parameters for the different algorithms. The first column lists the regression algorithms, and the second column the underlying algorithm used for the linear model estimation(s). Column 3-5 lists the model-parameters related to the number, position, and size of the basis functions. Model-parameters that must be pre-specified by the user -and are thus also meta-parameters of the algorithm -are highlighted in gray. Model-parameters that are automatically determined by the algorithm are left white. For example, in elm the number of basis functions is a meta-parameter that must be specified by the user (gray cell), whereas the centers of the basis functions are determined automatically through random sampling (white cell). The final column lists any further meta-parameters which are used for training, but not stored in the model parameters.

The settings of these meta-parameters can greatly influence the result of fitting, and they often depend on properties of the data. Exhaustive empirical comparisons are beyond the scope of this paper; our aim is to show that a wide variety of algorithms use different instances of the same unified model. In Table 4, we provide references to papers in which comparisons between the algorithms presented in this article are made.

Conclusion

In this article, we have described two main classes of models that regression algorithms usea weighted sum of basis functions and a mixture of linear models -and shown that the former is a special case of the latter. Therefore, rather than viewing every algorithm as using a different model, we see them as variations on a common theme, whereby these variations include the use of different basis functions, or using degenerate linear models with a = 0. This view thus forms a strong bridge between seemingly disparate regression algorithms in the field of machine learning, among which many ann models.

Going from a weighted sum of basis functions and a mixture of linear models requires a shift of interpretation. In the former, the sub-models are basis functions, and the weights are (degenerate) lines. But in a mixture of linear models, the sub-models are linear models, and the weights are determined by the basis functions. Due to this shift in interpretation, it may not be immediately obvious that the two models are the same.

A further insight is that, in a mixture of linear models, both the basis functions φ(x n , θ e ) and the linear model a e x + b e depend on x. Since we cannot regress against both at the same time. This is why lwr must perform multiple weighted linear regressions. By contrast, with a sum of weighted basis functions, the weights are linear with respect to the output of the basis functions, and can thus be learned in one least squares regression.

What are the implications of having a unified model for all these algorithms? First of all, it allows us to clearly separate questions about models from those about algorithms. As to the models, is it better to use global (e.g. sigmoid, cosine) or local (e.g. radial) features? When is it better to assume local linearity -and thus use a model based on locally weighted lines -or use weighted basis functions? Which types of models are easier to optimize in online learning, e.g. through backpropagation or reinforcement learning? What is the effect of mixing different types of basis functions [START_REF] Augusteijn | Evolving transfer functions for artificial neural networks[END_REF]? Orthogonal question related only to algorithms are: How should features be generated [START_REF] Rahimi | Weighted Sums of Random Kitchen Sinks: Replacing minimization with randomization in learning[END_REF]) -should they be randomly generated, preshaped through unsupervised learning, or optimized through constrained optimization or backpropagation? In which context should we use which norm? This work provides a conceptual framework in which to address these questions more clearly.

From our perspective, research topics whose full potential has not yet been explored include • using weighted linear models with global features for weighting,

• deriving a locally weighted regression variant of Support Vector Regression based on its specific Reference Algorithms compared [START_REF] Atkeson | Memory-based neural networks for robot learning[END_REF] lwr, mlffs [START_REF] Schaal | Receptive field weighted regression[END_REF] rfwr, mlffs, ls [START_REF] Vijayakumar | Locally weighted projection regression: An O(n) algorithm for incremental real time learning in high dimensional space[END_REF] lwpr, rfwr (Schwenker et al., 2001a) elm, BProp, svr [START_REF] Williams | Gaussian processes for machine learning[END_REF] lwpr, gpr, ls [START_REF] Grollman | Sparse incremental learning for interactive robot control policy estimation[END_REF] sogp, lwpr [START_REF] Nguyen-Tuong | Model learning with local gaussian process regression[END_REF] lgp, lwpr, svr, gpr [START_REF] Cederborg | Incremental local online Gaussian mixture regression for imitation learning of multiple tasks[END_REF] lgp, lwpr, gmr, svr, gpr [START_REF] Lammert | Locallyweighted regression for estimating the forward kinematics of a geometric vocal tract model[END_REF] lwr, mlffs [START_REF] Huang | Extreme learning machines: a survey[END_REF] elm, BProp, svr [START_REF] Gijsberts | Real-time model learning using incremental sparse spectrum gaussian process regression[END_REF] gpr, lwpr, i-ssgpr (Droniou et al., 2012a) xcsf, lwpr, irfrls [START_REF] Munzer | Non-linear regression algorithms for motor skill acquisition: a comparison[END_REF] lwr, gmr, rbfn, irfrls Table 4: List of articles, and the algorithms between which they make empirical comparisons.

norm, which does not seem to exist,

• mixing initial batch learning with subsequent incremental learning, as done by (Schwenker et al., 2001b),

• and investigating which types of regression algorithms yield models instances that are easily optimized in a reinforcement learning context [START_REF] Munzer | Non-linear regression algorithms for motor skill acquisition: a comparison[END_REF].

This latter topic is highly relevant to robotics, where skill learning is often split in two phases: 1) imitation: initializing a model of human motion from observations of human movement [START_REF] Calinon | Robot programming by demonstration[END_REF][START_REF] Stulp | Reinforcement learning with sequences of motion primitives for robust manipulation[END_REF] through offline batch regression 2) reinforcement learning: optimize the model (incrementally or with mini-batch) with respect to a cost function through trial and error exploration [START_REF] Stulp | Reinforcement learning with sequences of motion primitives for robust manipulation[END_REF][START_REF] Marin | Learning cost-efficient control policies with xcsf: generalization capabilities and further improvement[END_REF].

A very pragmatic advantage of a unified model arose whilst implementing several of the algorithms presented here. The classes representing the models of the different algorithms became so similar that it was an obvious step to implement a routine for casting them all to a unified model, which greatly simplified the code, and made it easier to exchange algorithms. This is not merely a low-level implementation advantage though. On a more abstract level, such modularity is an important pre-condition for realizing the vision of model-based machine learning [START_REF] Bishop | Model-based machine learning[END_REF]. research and innovation programme under grant agreement No 640891. We thank Vincent Padois, Nicolas Perrin, Ryan Lober and Mathieu Lefort for proofreading earlier versions of this article. Thank you Ajaya Kumar Pani for pointing out an error in the first paragraph (we switched "dependent" and "independent"), which we have corrected in this version.

Figure 1 :

 1 Figure 1: Radial Basis Function Network (rbfn) with D inputs and E radial basis functions. Unlabeled connections have weight 1.

Figure 2 :

 2 Figure 2: Generic flowchart for parametric regression.

Figure 4 :

 4 Figure 4: Network representation of a linear model.

Figure 5 :

 5 Figure 5: Illustration of least squares. Black dots represent the training examples, and the thick line is the learned function f (x). Vertical lines represent residuals.

Figure 6 :

 6 Figure 6: Visualization of the L 1 , L 2 , and svr norm.

C2 of Figure 7 .

 7 Performing weighted least squares E times for the different weighting functions leads to E local linear models a e x+b e , which are visualized as dashed lines in C1 . By choosing E Gaussian weighting functions φ e (x) = g(x, c e , Σ e ) with different centers c e in the input space, each linear model is thus responsible for fitting a different part of the data, as depicted in C1 in Figure7. This 'local responsibility' is the reason for the name Locally Weighted Regression.

Figure 7 :

 7 Figure 7: Illustration of four regression algorithms. Black dots represent the 20 training examples, and the thick (red) line is the learned function f (x). Vertical lines represent residuals, where the thickness of these lines indicates their weight. Dashed lines represent local models. The lower row visualizes the basis functions, where for lwr the thin line represents the unnormalized basis functions (normalization is not necessary for the other algorithms). The right graph plots the feature space for rbfn.Each data point in graph D3 is acquired by computing the value of both φ(x, c 1 ) and φ(x, c 2 ) (which can be read from D2), and plotting y against these two values. Because the 2D feature space φ(x, c 1 ), φ(x, c 2 ) is of a higher dimension that the 1D input space x, we now have a 3D plot.

Figure 9 :

 9 Figure 9: The rbfn model is a special case of the lwr model, with model parameters a = 0. Here, φe is an abbreviation of φ(x, θe)

  4.

Figure 10 :

 10 Figure 10: Meta-parameters for lwr: the number of weighting functions (E = 3 above) whose centers are spaced equidistantly in the range of the input data, and the intersection height of the unnormalized weighting functions (0.35 above).

  Figure 11:Left: Gaussian Mixture Model (gmm) with 3 Gaussians on the 20 example data points. This model is acquired through unsupervised learning on the joint input×output space. Right: Gaussian Mixture Regression, which yield the mean and variance in the output space Y when conditioning on the input space X.

  h e (x) ≡ φ(x, θ e ) are normalized Gaussian functions φ(x, θ e ) = π e g(x, c e , Σ e ) E e =1 π e g(x, c e , Σ e ) (48) as in (35), with the centers, widths and priors stored in the parameter vector θ e = (c e , Σ e , π e ). The model in (47) is equivalent to the unified model in (40). The relationship between the gmm and the unified model is illustrated in Figure 12.

Figure 12 :

 12 Figure 12: Relationship between gmm and the unified model. The basis functions are the (normalized) projections of the Gaussians in the gmm on the input space. The linear models have slopes ae = Σ e,Y X Σ -1e,Y and offsets be = µ e,Y -Σ e,Y X Σ -1 e,X µ e,X .

Figure 13 :

 13 Figure 13: Model tree representation, with linear models (top), and non-overlapping box-car weighting functions (bottom), one for each leaf in the tree.

Figure 14 :

 14 Figure 14: The function model used in krr as a network.

Figure 15 :

 15 Figure15: Mean and standard deviation (y q ± var(yq) ) of the Gaussian process given 10 (left) or 20 (right) example data points.

Figure 16 :

 16 Figure 16: Repetition of Figure 15, this time emphasizing that the mean in gpr is also a weighted sum of basis functions.

Figure 18 :

 18 Figure 18: Regression tree model, with degenerate linear models (top), and box-car basis functions (bottom).

Figure 19 :

 19 Figure19: Single-Hidden Layer Feedforward Network with D inputs and E hidden nodes. A bias may be introduced by adding an input neuron with input 1, but it has been left out for clarity.

Figure 20 :

 20 Figure20: Detailed view of one node in a artificial neural network. If the neuron is in a hidden layer, the output of the neuron o j will be the input for another neuron.

  ), and u e = D d=1 v ed x d , which results in w e ← w e + α(y n -f (x)) D d=1 v ed x d .

6. 1 .

 1 Classification based on Special Cases of the Unified Model Figure 21 presents a classification of the regression algorithms presented in this paper, based on what type of special case of the unified model they use. Note that this perspective considers only the model used to represent the function, and treats the algorithmic process as a black-box. For instance, least squares and weighted least squares output the exact same representation (one linear model), even if they use different methods to compute it, and may thus yield different parameterizations of this model.

  Figure 21: Classification of regression algorithms, based only on the model used to represent the underlying function.

  inputs/targets)

	uses mixture of linear	models model (Section 4)	M5 meta gmr meta xcsf meta lwpr meta rfwr meta lwr meta	M5 gmr xcsf lwpr rfwr lwr	unified model + model params	Evaluate input (novel)
	uses weighted sum of	basis functions (Section 5)	elm meta i-ssgpr meta irfrls meta svr meta gpr meta krr meta rbfn meta	elm i-ssgpr irfrls svr gpr krr rbfn		output (prediction)
			CART meta	CART

Table 3 :

 3 Gray: meta-parameters of the algorithms (partial model specifications in columns 3-5 and further algorithmic parameters in column 6) which the user must set. White: model parameters that the algorithm determines itself. Abbreviations: adap.=adapted, indiv.=from individual training examples, rand.=randomly sampled, CO=constrained optimization.

	Algorithm	Linear model estim.	Number of BFs (E)	Position of BFs (ce)	Size of BFs (Σe)	Algorithmic parameters
	Model: Mixture of Linear Models
	lwr	LS	fixed	fixed	fixed
	rfwr	RLS	adap.	indiv. adap.	many
	lwpr	NIPALS adap.	indiv. adap. many
	xcsf	RLS	adap.	adap.	adap. many
	gmr	EM	fixed	adap.	adap.
	M5	LS	adap.	adap.	adap.	n, σ, P , k
	Model: Weighted Sum of Basis Functions
	rbfn	LS	fixed	fixed	fixed
	krr	(RG)LS	indiv. indiv. fixed	λ
	gpr	LS	indiv. indiv.	fixed	k(.)
	svr	CO	adap.	indiv. fixed
	irfrls	RLS	fixed	rand.	rand.	γ, λ
	i-ssgpr	RLS	fixed	rand.	rand.
	CART	LS	adap.	adap.	adap.	n, σ, P , k
	elm	LS	fixed	rand.	rand.
	BProp	N/A	fixed	adap.	adap. α

Throughout the paper, we denote scalars as lowercase symbols (x), vectors as bold lowercase symbols (x) and matrices as bold uppercase symbols (X).

2 To facilitate presentation and visualizations, we mostly use Y ⊆ R 1 in this paper, i.e. the output is 1-dimensional, and the target values yn are scalars rather than vectors. Multivariate regression is presented in Section 2.4.

The full understanding of the solution requires knowledge of "slack variables" and "Lagrangian multipliers". Explaining these concepts is beyond the scope of this article, and we refer to the tutorial by[START_REF] Smola | A tutorial on support vector regression[END_REF] for the full derivation.

Methods for computing b are presented in[START_REF] Smola | A tutorial on support vector regression[END_REF].

We follow the convention by[START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF], where w ji denotes the weight from layer i to layer j. For consistency with the rest of the paper, we use d for the input layer, e for the hidden layer, and use only one output node for ease of presentation. Furthermore, we use different symbols for the weights in the different layers (v for input→hidden, w for hidden→output) to facilitate the explanation of applying least squares to learning w.

The relationship between radial basis function networks and multi-layer feedforward networks is further discussed in Section 5.10.4.

For a full derivation of the backpropagation rules, we refer to[START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. Here, we provide only a brief summary of the resulting rules.

Acknowledgments

This work was partially supported by the European Commission, within the CoDyCo project (FP7- ICT-2011-9, No.600716), and has received funding from the European Unions Horizon 2020