LS-SVR as a Bayesian RBF network

TL;DR

This paper establishes a formal connection between LS-SVR with RBF kernels and Bayesian RBF networks, enabling potential improvements by integrating Bayesian methods into LS-SVR.

Contribution

It explicitly states the theoretical similarities between LS-SVR and Bayesian RBF networks, providing a formal basis for future methodological enhancements.

Findings

Confirmed the theoretical correspondence through experiments

Demonstrated potential for improving LS-SVR with Bayesian techniques

Validated the approach on standard regression benchmarks

Abstract

We show theoretical similarities between the Least Squares Support Vector Regression (LS-SVR) model with a Radial Basis Functions (RBF) kernel and maximum a posteriori (MAP) inference on Bayesian RBF networks with a specific Gaussian prior on the regression weights. Although previous works have pointed out similar expressions between those learning approaches, we explicit and formally state the existing correspondences. We empirically demonstrate our result by performing computational experiments with standard regression benchmarks. Our findings open a range of possibilities to improve LS-SVR by borrowing strength from well-established developments in Bayesian methodology.

Figures

Equations

$$f({\mathbf{x}}_i)=⟨\mathbf{w},\phi ({\mathbf{x}}_i)⟩+b,\text{mbox}with\mathbf{w}\in {\mathbb{R}}^Q,b\in \mathbb{R},$$

$$J(\mathbf{w},e)=\frac{1}{2}\Vert \mathbf{w}{\Vert }_2^2+\gamma \frac{1}{2}\sum _{i=1}^N{e}_i^2,$$

$$y_i=⟨\mathbf{w},\phi ({\mathbf{x}}_i)⟩+b+e_i,i=1,2,\dots ,N,$$

$$L(\mathbf{w},b,\mathbf{e},\mathbf{\alpha })=\frac{1}{2}\Vert \mathbf{w}{\Vert }_2^2+\gamma \frac{1}{2}\sum _{i=1}^N{e}_i^2-\sum _{i=1}^N{\alpha }_i[⟨\mathbf{w},\phi ({\mathbf{x}}_i)⟩+b+e_i-y_i]$$

$$\left\{\begin{array}[]{ll}\frac{\partial L}{\partial\bm{w}}&=0~{}~{}\rightarrow~{}~{}\bm{w}=\sum_{i=1}^{N}\alpha_{i}\varphi(\bm{x}_{i}),\\ \frac{\partial L}{\partial e_{i}}&=0~{}~{}\rightarrow~{}~{}\sum_{i=1}^{N}\alpha_{i}=0,\\ \frac{\partial L}{\partial b}&=0~{}~{}\rightarrow~{}~{}\alpha_{i}=\gamma e_{i},\\ \frac{\partial L}{\partial\alpha_{i}}&=0~{}~{}\rightarrow~{}~{}\langle\bm{w},\varphi(\bm{x}_{i})\rangle+b+e_{i}-y_{i}=0.\end{array}\right.$$

$$\Psi \left[\begin{array}{c}b\\ \mathbf{\alpha }\end{array}\right]=\left[\begin{array}{c}0\\ \mathbf{y}\end{array}\right],\text{where }\Psi =\left[\begin{array}{cc}0& 1^{\top }\\ 1& \mathbf{\Omega }+{\gamma }^{-1}I\end{array}\right].$$

$${\mathbf{\theta }}_{\text{LS}}={\Psi }^{-1}\left[\begin{array}{c}0\\ \mathbf{y}\end{array}\right].$$

$$k({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp \{-c^2\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }^2\},$$

$$f(\mathbf{x})=\sum _{i=1}^N{\alpha }_{i}k(\mathbf{x},{\mathbf{x}}_i)+b.$$

$$h_i(\mathbf{x})=k(\Vert \mathbf{x}-{\mathbf{x}}_i\Vert ).$$

$$f(\mathbf{x})=\sum _{i=1}^N{\alpha }_i{h}_i(\mathbf{x})+b,$$

$$\psi \left[\begin{array}{c}b\\ \mathbf{\alpha }\end{array}\right]=\mathbf{y},\text{where }\psi =[1\mathbf{\Omega }],$$

$${\mathbf{\theta }}_{\text{RBF}}=({\psi }^{\top }\psi {)}^{-1}{\psi }^{\top }\mathbf{y}.$$

$$p(y_i|{\mathbf{x}}_i,\mathbf{\theta })=\mathcal{N}(y_i|f({\mathbf{x}}_{\mathbf{i}}),{\sigma }^2),$$

$$p(\mathbf{\theta }|\mathcal{D})=\frac{p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta })}{\int p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta })\mathrm{d}\mathbf{\theta }}\propto p(\mathcal{D}|\mathbf{\theta })p(\mathbf{\theta })$$

$$p(y|\mathbf{x},\mathcal{D})=\int p(y|(x),\mathbf{\theta })p(\mathbf{\theta }|\mathcal{D})\mathrm{d}\mathbf{\theta },$$

$$\mathbb{E}\{y\}=\int yp(y|\mathbf{x},\mathcal{D})\mathrm{d}y.$$

$$\displaystyle y_{n}\sim\mathcal{N}\Bigg{(}\sum_{i=1}^{N}\alpha_{i}k(\bm{x}_{n},\bm{x}_{i})+b,\,\sigma^{2}=1\Bigg{)},$$

$$[b,\mathbf{\alpha }{]}^{\top }\sim \mathcal{N}(\mu ,\Sigma ),$$

$${\Sigma }^{-1}=\left[\begin{array}{cc}ϵ& {\gamma }^{-1}{1}^{\top }\\ {\gamma }^{-1}1& 11^{\top }+2{\gamma }^{-1}\Omega +{\gamma }^{-2}I\end{array}\right],$$

$$\mu ={\gamma }^{-1}\Sigma \left[\begin{array}{c}0\\ y\end{array}\right].$$

$$\frac{ϵ}{1+ϵs_1({\Psi }^{\top }\Psi {)}^{-1}{s}_1^{\top }}({\Psi }^{\top }\Psi {)}^{-1}{s}_1{s}_1^{\top }{\Psi }^{-1}\left[\begin{array}{c}0\\ \mathbf{y}\end{array}\right],$$

$${\mathbf{\varphi }}_i=[1,k({\mathbf{x}}_i,{\mathbf{x}}_1),\dots ,k({\mathbf{x}}_i,{\mathbf{x}}_N){]}^{\top }.$$

$$\displaystyle p\big{(}b,\bm{\alpha}\,|\,\mathcal{D},\sigma^{2}=1/2\big{)}$$

$$=\mathcal{N}(\mathbf{\theta }|{\mu }^{\star },{\Sigma }^{\star }),$$

$$\Sigma^{\star}=\Bigg{(}\Sigma^{-1}+\Phi^{\top}\Phi\Bigg{)}^{-1}\text{ and }\mu^{\star}=\Sigma^{\star}\Bigg{(}\Phi^{\top}\bm{y}+\gamma^{-1}\begin{bmatrix}0\\ \bm{y}\end{bmatrix}\Bigg{)}.$$

$${\mu }^{\star }={\Sigma }^{\star }{\left[\begin{array}{cc}0& 1^{\top }\\ 1& \mathbf{\Omega }+{\gamma }^{-1}I\end{array}\right]}^{\top }\left[\begin{array}{c}0\\ \mathbf{y}\end{array}\right]={\Sigma }^{\star }{\Psi }^{\top }\left[\begin{array}{c}0\\ \mathbf{y}\end{array}\right].$$

$$\begin{bmatrix}N&\mathbf{1}^{\top}\Big{(}\Omega+\gamma^{-1}I\Big{)}\\ \Big{(}\Omega+\gamma^{-1}I\Big{)}\mathbf{1}&\mathbf{1}\mathbf{1}^{\top}+\Big{(}\Omega+\gamma^{-1}I\Big{)}^{\top}\Big{(}\Omega+\gamma^{-1}I\Big{)}\end{bmatrix}.$$

$$\Phi =\left[\begin{array}{cc}1& \Omega \end{array}\right],$$

$${\Phi }^{\top }\Phi =\left[\begin{array}{c}{1}^{\top }\\ {\Omega }^{\top }\end{array}\right]\left[\begin{array}{cc}1& \Omega \end{array}\right]=\left[\begin{array}{cc}N& 1^{\top }\Omega \\ {\Omega }^{\top }1& {\Omega }^{\top }\Omega \end{array}\right].$$

$${\Sigma }^{-1}={\Psi }^{\top }\Psi +\left[\begin{array}{c}\sqrt{ϵ}\\ 0\\ ⋮\\ 0\end{array}\right]\left[\begin{array}{c}\sqrt{ϵ},0,\dots ,0\end{array}\right]-{\Phi }^{\top }\Phi ,$$

Full text

1 Introduction

Statistical learning theory has been studied for general function estimation from data since the late 1960’s [25]. However, it was only widely adopted in practice after the introduction of the learning algorithms known as Support Vector Machines (SVMs) [26]. Using the so-called kernel trick, which replaces dot products between features and model parameters by evaluations of a kernel function, SVMs can learn nonlinear relations from training patterns by solving a convex optimization problem [18].

An important variant of the SVM is the Least Squares Support Vector Machine (LS-SVM) [23], which is obtained by making all data points support-vectors. LS-SVM avoids the constrained quadratic optimization step of standard SVMs by replacing the training procedure with one that reduces to solving a system of linear equations, which can be performed via ordinary least squares.

The first SVM formulation was derived for classification tasks, but it has been readily adapted to tackle regression problems, being usually named Support Vector Regression (SVR) [7]. Similarly, the regression counterpart of LS-SVM is the LS-SVR [22].

It is well known that the LS-SVR model is closely related to the kernel ridge regression (KRR) formulation [17, 4]. As detailed in the forthcoming sections, those models only differ by the inclusion of a bias term. Interestingly, one can find additional relations by recalling that KRR itself can be seen as the maximum a posteriori (MAP) solution of a Gaussian process (GP) model [16]. Standard SVMs [19, 20] and LS-SVMs [12] also present a close relation to the MAP solution of a GP classifier. Similar relations have also been explored in regression settings with the SVR model [9, 3]. On the other hand, it can be shown that GPs are a limit for Bayesian neural networks (BNNs) with infinitely many hidden neurons with Gaussian priors in their weights [15]. The kernel functions related to such limit was derived in [27] for standard multilayer perceptron (MLP) with sigmoid activation functions and Gaussian radial basis functions (RBF) neural networks.

Many works have explored theoretical relations between learning models. For instance, Bayesian interpretations for SVMs are proposed in [28] as an effort to push forward the interface between large-margin learning and Bayesian nonparametrics. Recently, authors in [8] have presented a method for mapping a GP prior to the weights’ prior of a BNN. The authors emphasize how their proposal is able to bring properties of a variety of kernels without explicitly considering the weights parameter space. In [10], it is presented how a residual convolutional neural network with an appropriate prior over the weights and biases behave similarly to a GP. They use such correspondence to introduce a convolutional GP kernel with few hyperparameters. The authors in [14] have also tackled the relation between deep neural networks and GPs, deriving an equivalent kernel function for the latter. Those works motivate us to explore further the missing links between learning approaches.

Given the previous observations, our work aims to explore the relations between LS-SVR and Bayesian RBF neural networks. Such theoretical direct connection is not readily available in the literature, and it is required to provide ways to improve least squares based kernel methods with recent advances in the domains of Bayesian methods and neural networks.

Some works, such as [24, 21], have already considered the Bayesian formalism in the context of LS-SVM models, but only for the sake of hyperparameter tuning. In [5] the authors aim to provide confidence intervals for the model predictions. However, frequentist machinery is used, which makes is dependent on the parameter estimation procedure. The probabilistic interpretation of the LS-SVM classifier presented in [12] is the closest work in our direction. However, in that work the authors propose a Bayesian inference framework for the model bias and weights in the LS-SVM primal space. An inference procedure for the model hyperparameters is also presented in [12]. As we will detail later, our main result considers the task of performing inference in the variables of the LS-SVR dual space, in the regression setting.

Thus, in the present work, instead of only considering some evidence-based inference step, we tackle the regression setting with LS-SVR using a RBF kernel and directly relate its model formulation with the RBF neural network within the Bayesian view. Fig. 1 illustrates the aforementioned learning approaches connections, as well as the one we cover in the remaining sections.

As follows we summarize both models and derive relations between their formulations. Afterwards, we perform computational experiments to empirically assess our theoretical results. Finally, we draw concluding remarks and discuss the implications of our findings.

2 Preliminaries

In the subsequent subsections, we review the basics of LS-SVR and RBF networks while setting up relevant notation.

2.1 LS-SVR in a nutshell

Let a dataset be composed by $N$ inputs $\bm{x}_{i}\in\mathbb{R}^{D},\forall i\in\{1,\cdots,N\}$, and their correspondent outputs $y_{i}\in\mathbb{R}$. We can map the inputs into the outputs using some unknown nonlinear function $f$:

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the dot-product, $\bm{w}$ is a vector of weights, $b$ is a bias and $\varphi(\cdot):\mathbb{R}^{D}\rightarrow\mathbb{R}^{Q}$ is a nonlinear map into some $Q$-dimensional feature space. In the LS-SVR framework, we consider the minimization of the following functional [17, 22]:

[TABLE]

subject to the equality constraints

[TABLE]

where $e_{i}=y_{i}-f(\bm{x}_{i})$ is the error due to the $i$-th input pattern and $\gamma>0$ is a regularization parameter.

The related Lagrangian function of the optimization problem is

[TABLE]

where the $\alpha_{i}$ are the Lagrange multipliers. The Karush-Kuhn-Tucker optimality conditions are given by

[TABLE]

The optimal solution $\bm{\theta}_{\text{LS}}$ depends only on $\bm{\theta}=\begin{bmatrix}b\\ \bm{\alpha}\end{bmatrix}$ and is obtained by solving the following linear system:

[TABLE]

The solution if given by:

[TABLE]

In Eq. 6 we have defined the vectors $\bm{y}=[y_{1},y_{2},\ldots,y_{n}]^{\top}$, $\bm{1}=[1,1,\ldots,1]^{\top}$, $\bm{\alpha}=[\alpha_{1},\alpha_{2},\ldots,\alpha_{n}]^{\top}$. We have also defined the kernel matrix $\bm{\Omega}\in\mathbb{R}^{N\times N}$, whose entries are $\Omega_{i,j}=k(\bm{x}_{i},\bm{x}_{j})=\langle\varphi(\bm{x}_{i}),\varphi(\bm{x}_{j})\rangle$. The so-called kernel function $k(\cdot,\cdot)$ is usually chosen to be a RBF (or Gaussian), expressed by

[TABLE]

where $\|\cdot\|$ represents the Euclidean norm and $c$ is a kernel hyperparameter.

Predictions for an input $\bm{x}$ with the LS-SVR model can be performed using the optimized $\bm{\alpha}$ and $b$ values:

[TABLE]

2.2 RBF networks in a nutshell

Radial basis function networks are neural models with a simple structure [13]: a given input pattern $\bm{x}\in\mathbb{R}^{D}$ is nonlinear transformed by a hidden layer where each unit $i$ is given by

[TABLE]

The model output is obtained by a linear transformation:

[TABLE]

where $\alpha_{i}$ is the $i$-th weight of the output layer and $b$ is the bias term. Although some classic RBF formulations do not present a bias term [13], many authors do include it in the output layer [2]. We note that the latter case is more general, since we could make $b=0$.

If we choose $h_{i}(\bm{x})=k(\|\bm{x}-\bm{x}_{i}\|)=\exp\{-c^{2}\|\bm{x}-\bm{x}_{i}\|^{2}\}$, the output of the RBF network would be the same as in Eq. 8. Note that one could opt to build a hidden layer with $M<N$ hidden units, due to inherent data redundancies. In that case, a subset of the training data or the centroids found by a clustering method, e.g., the K-means algorithm, could be used as the centers of the hidden layer [13].

From Eq. 9 we can see that the optimal solution for $\begin{bmatrix}b\\ \bm{\alpha}\end{bmatrix}$ is obtained from the following linear system:

[TABLE]

where we have defined the matrix $\bm{\Omega}\in\mathbb{R}^{N\times N}$, whose entries are computed by $\Omega_{i,j}=k(\|\bm{x}_{i}-\bm{x}_{j}\|)$. The solution $\bm{\theta}_{RBF}$ is given by:

[TABLE]

In the present paper we are interested in the Bayesian treatment of the RBF network. In that formulation, we begin by defining the observation model

[TABLE]

where $f(\bm{x},\bm{\theta}|\mathcal{D})$ is the same as in the right side of Eq. 9. We complete the model by placing a prior distribution $p(\bm{\theta})$ on the parameters $\bm{\theta}=(b,\bm{\alpha})^{\top}$ and proceed to compute the posterior given the training data:

[TABLE]

Once in possession of the posterior $p(\bm{\theta}|\mathcal{D})$, we can compute the distribution of the network output $y$ for an input $\bm{x}$ [15]:

[TABLE]

for which the mean, i.e., the expected network output, is given by:

[TABLE]

We note that one could also place priors on kernel parameters, on the variance $\sigma^{2}$ and on the centers locations in the hidden layer [1], but for the following sections choosing only $p(\bm{\theta})$ is enough.

3 Main result

We now present a result that shows an equivalence in limit between the estimation procedure of the LS-SVR and the MAP estimation of RBF networks. For this purpose, we introduce the concept of a Bayesian $\epsilon$-LS-SVR.

Definition** (Bayesian $\epsilon$-LS-SVR).**

Given $\epsilon>0$, a Bayesian $\epsilon$-LS-SVR is a Bayesian RBF network model defined using all training points as centroids in the hidden layer and is described as

[TABLE]

where the covariance matrix $\Sigma$ is such that

[TABLE]

and the mean vector $\mu$ is given by:

[TABLE]

And show that, there is an analytic form for the exact difference between MAP estimates in an Bayesian $\epsilon$-LS-SVR and its traditional counterpart.

Lemma.

Given $\epsilon>0$, the parameters $\bm{\theta}_{\text{LS}}$ of an $LS-SVR$ and the MAP estimates $\mu^{\star}$ of a Bayesian $\epsilon$-LS-SVR defined under the same aforementioned conditions, the gap $\bm{\theta}_{\text{LS}}-\mu^{\star}$ is exactly:

[TABLE]

where $s_{1}$ is a $(N+1)$-dimensional column vector of zeros except for a $1$ in its first element.

Proof.

Let $\bm{\theta}=[b,\bm{\alpha}]^{\top}$ and $\Phi\in\mathbb{R}^{N\times(N+1)}$ with its $i$-th line being denoted as $\bm{\phi}_{i}$ such that

[TABLE]

The posterior of the aforementioned model is given by:

[TABLE]

with parameters $\mu^{\star}$ and $\Sigma^{\star}$:

[TABLE]

As the mean and mode of a multivariate normal random variable coincide, $\mu^{\star}$ is also the MAP estimate for $\bm{\theta}=[b,\alpha]^{\top}$. Note that:

[TABLE]

We can also write $\Psi^{\top}\Psi$ as

[TABLE]

Observe that $\Phi$ is essentially the kernel matrix $\Omega$ padded with ones on the left, i.e.:

[TABLE]

so that:

[TABLE]

Notice that the matrix $\Sigma^{-1}$ is given by

[TABLE]

and thus, the matrix $(\Sigma^{\star})^{-1}$ becomes

[TABLE]

which can we rewritten as:

[TABLE]

where $s_{1}$ is a $(N+1)$-dimensional column vector of zeros except for a $1$ in its first element.

Applying Sherman-Morrison’s matrix inversion identity on $(\Sigma^{\star})^{-1}$ and rearranging the terms we get the expression below:

[TABLE]

Thus, we can rewrite the vector $\mu^{\star}$ as

[TABLE]

Finally, we are able to compute the gap $\bm{\theta}_{\text{LS}}-\mu^{*}$:

[TABLE]

Since $\Psi$ is square and its inverse $\Psi^{-1}$ exists, which is true given its definition in Eq. (6), we get $\Psi^{-1}=(\Psi^{\top}\Psi)^{-1}\Psi^{\top}$ and we obtain the desired result. ∎

Analyzing this gap in the limit as $\epsilon$ shrinks, we can derive the following equivalency result, which incurs directly in the subsequent corollary.

Theorem.

As $\epsilon$ decreases towards zero, the parameters $\bm{\theta}_{\text{LS}}$ of LS-SVR coincide with the MAP estimate $\mu^{\star}$ of the Bayesian $\epsilon$-LS-SVR defined over the same dataset with same hyper-parameters.

Proof.

Using Lemma 1, as $\epsilon\xrightarrow[]{+}0$, the gap $\bm{\theta}_{\text{LS}}-\mu^{\star}$ tends to

[TABLE]

which equals zero. Thus, as $\epsilon$ decreases towards zero, the gap $\bm{\theta}_{\text{LS}}-\mu^{\star}$ vanishes and $\mu^{\star}$ coincides with the LS-SVR estimates. ∎

Corollary.

For any $\bm{x}\in\mathbb{R}^{D}$, as $\epsilon$ decreases towards zero, the LS-SVR prediction using Equation 8 coincides with the ones obtained using the same rule but using the MAP estimates of the respective Bayesian $\epsilon$-LS-SVR.

4 Discussion and Limitations

The results detailed in Section 3 lay down the theoretical similarities of the LS-SVR model with the RBF kernel and the Bayesian RBF network. We were able to show how the outputs of both models can be matched, which itself depends on the same parameter inference procedure. Our probabilistic presentation, for instance, incorporates the $\epsilon$ hyperparameter in the introduced $\epsilon$-LS-SVR formulation, which corresponds to the precision of the prior given to the model bias. In the limit $\epsilon\rightarrow 0$, the bias has an improper uniform prior.

As mentioned in the Introduction, differently from [12] for the LS-SVM, we consider a Bayesian approach directly in the dual formulation of the LS-SVR. We emphasize that this choice is comprehensive enough and sufficient for the proposed analysis, given that kernel-based methods in general apply only dual variables to perform their learning and predictions. Furthermore, although not explored in the current work, the hyperparameter inference level pursued in [12] could be readily considered by choosing a (hyper)prior for $\epsilon$ and the regularization hyperparameter $\gamma$. However, the effective number of parameters analysis performed in [12] is not straightforwardly ported to our approach, since we bypass the primal formulation. The latter also prevents us to make more detailed geometric interpretations in the primal space.

One could question if the similarities between the expressions we have derived actually translate to model equivalence. Full model equivalence is a concept difficult to define. However, we have demonstrated that both formulations have the same parameter estimation problem, the same (Bayesian) model selection approach and the same predictive rules. Although we have strayed from the support vectors analogy usually defined in the primal weight formulation, which causes the above mentioned interpretation limitations, the lack of such analogy does not disqualify our claims.

Finally, while our derivations implied the use of the RBF kernel, the standard LS-SVR framework allows for any valid kernel function to be used. Even though the choice of a non-RBF kernel causes the RBF network correspondence to be broken, it is important to note that the resulting variant would present similar inference expressions and be comparable to a Bayesian generalized linear model [11]. A more thorough analysis of the latter implication is left for future investigations.

5 Illustrative experiments

In this section we report some empirical results to illustrate the presentation detailed so far. First, we show how the MAP estimate for the $\epsilon$-LS-SVR can get arbitrarily close to the standard LS-SVR weight vector solution $\bm{\theta}_{\text{LS}}$ as we choose small $\epsilon$ values. Then, we perform a simple experiment to show how to make Bayesian predictions with the formulated $\epsilon$-LS-SVR model.

5.1 Comparison of the $\epsilon$-LS-SVR and LS-SVR models

To empirically illustrate our main theoretical result, we performed a set of experiments using $4$ datasets a available at UCI machine learning repository [6]. We progressively decrease the $\epsilon$-LS-SVR model $\epsilon$ value and perform MAP estimation for its parameters. After that, we compute the Average Root Square Error (ARSE) between the estimates of both methods, given by $\text{ARSE}=\frac{1}{N}\sum_{i=1}^{N}\left|\mu_{i}^{\epsilon\text{-LS-SVR}}-y_{i}^{\text{LS-SVR}}\right|$, where $\mu_{i}^{\epsilon\text{-LS-SVR}}$ is the predicted mean of the $\epsilon$-LS-SVR and $y_{i}^{\text{LS-SVR}}$ is the prediction of the standard LS-SVR. We set $\gamma=0.5$ and $c^{2}=1$.

The obtained results are presented in Fig. 2. As can be noticed, the estimates of both methods converge to the same values as $\epsilon$ decreases. A similar result is obtained when comparing the parameters estimated by the two learning strategies. Fig. 3 shows the Euclidean distance between the weight vectors obtained by the $\epsilon$-LS-SVR and LS-SVR models. These results coincide with the behaviour outlined in the theorem that we have introduced in Section 3.

5.2 Bayesian predictions with the LS-SVR model

Once in possession of the $\epsilon$-LS-SVR posterior, one can compute statistics of an arbitrary test function. For example, we can compute the expected value and variance of the LS-SVR prediction rule, in Eq. (8), as

[TABLE]

We illustrate these capabilities by modeling the univariate (normalized) sinc function:

[TABLE]

For this purpose, we created a dataset by taking $1200$ input points uniformly drawn in the range $x\in[-2\pi,2\pi]$ and evaluating the sinc function in each of them. We used $c=\gamma=1$ and $\epsilon=10^{-4}$. After the $\epsilon$-LS-SVR posteriori has been computed, we used Eqs. (15)–(16) to predict the values of the sinc function for new values in the range $x\in[-3\pi,3\pi]$. The results are shown in Fig. 3. Note that, as expected, variances are lower in the regions contemplated in the observed data than outside it. This experiment has been repeated ten, times yielding similar results.

6 Conclusion

We have presented a theorem showing a relationship between LS-SVR and Bayesian RBF networks. In a nutshell, the result provides a way to build a Bayesian RBF network which, when taking MAP estimates, presents theoretical similarities to a LS-SVR model. Such relationship opens a path to improve LS-SVR models with well-established advances in Bayesian methodology.

As a byproduct, we have introduced $\epsilon$-LS-SVR, a very particular Bayesian RBF network, which outlines new foundations for building Bayesian LS-SVRs. The posterior from our model can be used, for instance, to inform decision making processes (e.g. computing expected utilities) or to provide distributions for a new output $y^{\prime}$ given a new input vector $\bm{x}^{\prime}$.

A natural improvement would be to automate the selection of the regularizer $\gamma^{-1}$, the bias prior precision $\epsilon$ and the observation noise $\sigma^{2}$. This can be easily achieved by placing priors on these parameters. Besides freeing us from laborious parameter tuning, priors can also contain the subjective knowledge of a domain expert. While the introduction of such priors might break conjugacy, we can still resort to full Bayesian inference (e.g., using Markov Chain Monte Carlo), which would allow us to still recognize our uncertainty about the estimated parameters, in the form of a posterior distribution. Another promising direction would be to transfer information between different $\epsilon$-LS-SVR models using hierarchical models. All those investigations are left for future work.