Improving the Speed of Support Vector Regression Using Regularized Least Square Regression

One of the most commonly used eager learning method is least-square method (LS) [1], which first was proposed by Carl Friedesh Gauss. In the LS method, the best estimated function is a function that the sum of square of the difference between the observed values and the corresponding estimated values with the estimated function is minimum. Later, on the basis of the LS method, the Regulation Least Square (RLS) method was proposed [2, 3]. The RLS uses a regularization term to reduce the complexity of estimated function and to increase its generalization ability.

ε-insensitive Support Vector Regression (ε-SVR) [4] supposes that the noise in the output variable is at most ε. Therefore, the center of a tube with radius ε, which is used as the estimated function, is determined in a way that the training data are located in that tube. Therefore, this method is robust to such noisy data. In v-SVR [5], the maximum noise level or the tube radius is considered as a model variable, which is determined in the process of solving the model. In both ε-SVR and v-SVR, the maximum noise level or the radius of the tube is considered to be constant or independent in terms of the training data characteristics. To alleviate this shortcoming, Support Vector Interval Regression Networks (SVIRNs) was proposed [6]. In SVIRNs, the upper and lower bounds of the tube are automatically determined by two independent RBF networks. In this model, the radius of the tube can change by changing the features of the training data. Therefore, SVIRNs is robust to a noise whose maximum level depends on training data features. These two networks are initialized using ε-SVR. v-Support Vector Interval Regression Networks (v-SVIRN) [7] determines the upper and lower bounds of the tube by only one neural network or one model. In Twin SVR (TSVR) [8], estimated function is determined by solving two small Quadratic Programming Problems (QPPs), while ε-SVR solves one large QPP. Solving the TSVR model needs less time than solving ε-SVR model. Later, in ε-TSVR [9], a regularization term was used to reduce the complexity of estimated function and to increase its generalization ability. The regularization term of ε-TSVR is L2-norm of the estimated function weight vector, while in L1-TWSVR model [10], L1-norm of the estimate function weight vector was used as the regularization term. L1-TWSVR is a linear programming, while ε-TSVR is a QPP. As well, the estimated function of L1-TWSVR is sparser than that of ε-TSVR. Therefore, the testing time of L1-TWSVR model is less than that of ε-TSVR. Pair v-SVR is extended version of ε-TSVR which has additional advantage of using parameter v for controlling the bounds on fraction of support vectors and error [11].

All the mentioned SVR-based methods assume that the noise distribution is symmetrical, therefore, the center of a tube which contain noisy data is used as the estimated function. In other words, the estimated function is determined in such a way that approximately half of the training data are positioned above it, and the remaining data are positioned below it. In the asymmetric v-SVR method [12], the estimated function is determined in such a way that an arbitrary number of training data are positioned above it and the remaining data are positioned below it. The size of each of these two parts of data is determined based on the prior knowledge about the noise distribution. Next, twin version of asymmetric v-SVR, i.e. asymmetric v-TSVR [13], was proposed to improve the runtime of the asymmetric v-SVR model.

Among the mentioned methods, the LS and RLS methods have the lowest runtime because each of them is an n + 1 variable problem with a close-form solution which is obtained by inverting an n × n matrix, but they are too sensitive to noisy data. Against this, SVR-based models have a good robustness against noise but do not have a good runtime because for example to solve ε-SVR problem, a quadratic programming problem with 2×n variables and linear constraints must be solved which is too time consuming. In this paper, to improve the runtime of ε-SVR, first, an initial estimated function is obtained using the RLS method. Then, unlike the ε-SVR model, which uses all training data to determine the lower and upper limits of the tube, our proposed method uses the initial estimated function for determining the tube and the final estimated function. In other words, the data below and above the initial estimated function are used to estimate the upper limit and the lower limit of the tube, respectively. Thus, the number of the model constraints and, consequently, the model runtime are reduced. The experiments carried out on 15 benchmark data sets confirm that our proposed method is faster than ε-SVR, ε-TSVR and pair v-SVR, and its accuracy are comparable with that of ε-SVR, ε-TSVR and pair v-SVR.

In continue, in Section 2, RLS and ε-SVR are introduced. Then, in Section 3, our proposed method is presented. In Section 4, the results of the experiments are presented and Section 5 provides conclusion.

The RLS have the good runtime, but it is sensitive to noisy data. Against this, ε-SVR model-based models are robust to noise but do not have a good runtime. In this paper, to improve the runtime of ε-SVR, first, an initial estimated function is obtained using the RLS method. Then, unlike the ε-SVR model, which uses all training data to determine the lower and upper limits of the tube, our proposed method uses the initial estimated function for determining the tube and the final estimated function. In other words, the data below and above the initial estimated function are used to estimate the upper and the lower limits of the tube, respectively. Thus, the number of the model constraints and, consequently, the model runtime are reduced.

Let $\mathrm{f}(\mathrm{x})=\sum_{\mathrm{i}=1}^{\mathrm{n}} \alpha_{\mathrm{i}} \mathrm{k}\left(\mathrm{x}_{\mathrm{i}}, \mathrm{x}\right)+\mathrm{b}$ be initial estimated function obtained by using the RLS model based on the training data {(x₁,y₁), (x₂,y₂), …, (x_n,y_n)}. If f(x_i )≤0 then (x_i,y_i) is a point below the function f(.); otherwise, (x_i,y_i) is a point above the function f(.). Let I_a be the index of the training data above the function f, and I_b be the index of the training data below the function f, and |I_a|+|I_b|=n. Our proposed model for determining a tube with the radius ε and the center $\tilde{f}(x)=w^{T} \varphi(x)+b$ is as follows:

$\operatorname{Min}_{\mathbf{w}, \mathbf{b}, \bar{\xi}} \frac{1}{2}\|\mathbf{w}\|^{2}+\mathrm{C}_{2} \sum_{i=1}^{\mathrm{n}} \xi_{\mathrm{i}}$

subject to $\left\{\begin{array}{l}\mathrm{y}_{\mathrm{i}}-\tilde{\mathrm{f}}\left(\mathrm{x}_{\mathrm{i}}\right) \leq \varepsilon+\xi_{\mathrm{i}}, \mathrm{i} \in \mathrm{I}_{\mathrm{a}} \\ \tilde{\mathrm{f}}\left(\mathrm{x}_{\mathrm{i}}\right)-\mathrm{y}_{\mathrm{i}} \leq \varepsilon+\xi_{\mathrm{i}}, \mathrm{i} \in \mathrm{I}_{\mathrm{b}} \\ \xi_{\mathrm{i}} \geq 0, \quad \mathrm{i}=1,2, \ldots, \mathrm{n}\end{array}\right.$     (15)

where, ξ allows outliers are not positioned outside the tube. The number of outliers positioned outside the tube is controlled by the parameter C₂. Unlike ε-SVR which uses all training data to determine the lower and upper bounds of the tube, in our proposed model, only the upper half of the training data is used to determine the upper limit of the tube, and also only the lower half of the training data is used to determine the lower limit of the tube. The upper half and lower half of training data are determined based on the initial estimated function obtained by the RLS. The Lagrange function of the model (15) is as follows:

$\mathrm{L}=\frac{1}{2}\|\mathrm{w}\|^{2}+\mathrm{C} \sum_{\mathrm{i}=1}^{\mathrm{n}} \xi_{\mathrm{i}}$

$-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{a}}} \alpha_{\mathrm{i}}\left(\varepsilon+\xi_{\mathrm{i}}-\mathrm{y}_{\mathrm{i}}+\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)+\mathrm{b}\right)$

$-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{b}}} \alpha_{\mathrm{i}}\left(\varepsilon+\xi_{\mathrm{i}}+y_{\mathrm{i}}-\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)-\mathrm{b}\right)-\sum_{\mathrm{i}=1}^{\mathrm{n}} \gamma_{\mathrm{i}} \xi_{\mathrm{i}}$     (16)

where, γ_i and α_i are Lagrange coefficients. At the optimal point of the model (15) we have:

$\frac{\partial \mathrm{L}}{\partial \mathrm{w}}=0 \rightarrow \mathrm{w}=\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{a}}} \alpha_{\mathrm{i}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{b}}} \alpha_{\mathrm{i}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)$     (17)

$\frac{\partial \mathrm{L}}{\partial \mathrm{b}}=0 \rightarrow \sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{a}}} \alpha_{\mathrm{i}}-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{b}}} \alpha_{\mathrm{i}}=0$     (18)

$\frac{\partial \mathrm{L}}{\partial \xi}=0 \rightarrow \alpha_{\mathrm{i}}=\mathrm{C}_{2}-\gamma_{\mathrm{i}}, \quad \mathrm{i}=1,2, \ldots, \mathrm{n}$     (19)

$\mathrm{y}_{\mathrm{i}}-\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)-\mathrm{b} \leq \varepsilon+\xi_{\mathrm{i}}, \quad \mathrm{i} \in \mathrm{I}_{\mathrm{a}}$     (20)

$\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)+\mathrm{b}-\mathrm{y}_{\mathrm{i}} \leq \varepsilon+\xi_{\mathrm{i}}, \quad \mathrm{i} \in \mathrm{I}_{\mathrm{b}}$     (21)

$\alpha_{\mathrm{i}}\left(\varepsilon+\xi_{\mathrm{i}}-y_{\mathrm{i}}+\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)+\mathrm{b}\right)=0, \quad \mathrm{i} \in \mathrm{I}_{\mathrm{a}}$     (22)

$\alpha_{i}\left(\varepsilon+\xi_{i}+y_{i}-w^{T} \varphi\left(x_{i}\right)-b\right)=0, \quad i \in I_{b}$     (23)

$\gamma_{i} \xi_{i}=0, \quad i=1,2, \ldots, n$     (24)

$\xi_{\mathrm{i}}, \alpha_{\mathrm{i}}, \gamma_{\mathrm{i}} \geq 0, \quad \mathrm{i}=1,2, \ldots, \mathrm{n}$     (25)

By substituting Eq. (17) to Eq. (24) into the Lagrange function, this function is simplified as follows:

$\mathrm{L}=\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{o}}} \alpha_{\mathrm{i}} \mathrm{y}_{\mathrm{i}}-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{u}}} \alpha_{\mathrm{i}} \mathrm{y}_{\mathrm{i}}-\varepsilon \sum_{\mathrm{i}=1}^{\mathrm{n}} \alpha_{\mathrm{i}}$

$-\frac{1}{2}\left(\sum_{i \in I_{0}} \sum_{j \in I_{0}} \alpha_{i} \alpha_{i} k\left(x_{i}, x_{j}\right)+\sum_{i \in I_{u}} \sum_{j \in I_{u}} \alpha_{i} \alpha_{j} k\left(x_{i}, x_{j}\right)-2 \sum_{i \in I_{0}} \sum_{j \in I_{u}} \alpha_{i} \alpha_{j} k\left(x_{i}, x_{i}\right)\right)$     (26)

In addition, according to Eq. (19), and since α_i,γ_i≥0, we have:

$0 \leq \alpha_{\mathrm{i}} \leq \mathrm{C}_{2}, \quad \mathrm{i}=1,2, \ldots, \mathrm{n}$     (27)

Therefore, the dual of the model (15) is as follows:

$\max _{\alpha} \sum_{i \in I_{a}} \alpha_{i} y_{i}-\sum_{i \in I_{b}} \alpha_{i} y_{i}-\varepsilon \sum_{i=1}^{n} \alpha_{i}$

$-\frac{1}{2}\left(\sum_{i \in I_{a}} \sum_{j \in I_{a}} \alpha_{i} \alpha_{j} k\left(x_{i}, x_{j}\right)+\sum_{i \in I_{b}} \sum_{j \in I_{b}} \alpha_{i} \alpha_{j} k\left(x_{i}, x_{j}\right)\right.$

$\left.-2 \sum_{i \in I_{a}} \sum_{j \in I_{b}} \alpha_{i} \alpha_{j} k\left(x_{i}, x_{j}\right)\right)$

subject to $\left\{\begin{array}{l}\sum_{i \in I_{a}} \alpha_{i}-\sum_{i \in I_{b}} \alpha_{i}=0 \\ 0 \leq a_{i} \leq C_{2}, \quad i=1,2, \ldots, n\end{array}\right.$     (28)

If we define:

$\beta_{i}=\left\{\begin{array}{ll}\alpha_{i}, & i \in I_{0} \\ -\alpha_{i}, & i \in I_{u}\end{array}\right.$     (29)

then, the model (28) can be written as follows:

$\max _{\beta} \sum_{i=1}^{n} \beta_{i} y_{i}-\varepsilon\left(\sum_{i \in I_{a}} \beta_{i}-\sum_{i \in I_{b}} \beta_{i}\right)$$-\frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} \beta_{i} \beta_{j} k\left(x_{i}, x_{j}\right)$

subject to $\left\{\begin{array}{l}\sum_{i=1}^{n} \beta_{i}=0 \\ 0 \leq \beta_{i} \leq C_{2}, \quad i \in I_{a} \\ -C_{2} \leq \beta_{i} \leq 0, \quad i \in I_{b}\end{array}\right.$     (30)

which is a convex quadratic programming problem. After solving the model (30) and obtaining the optimal value of $\beta$, the optimal values of $\alpha$ can be determined using Eq. (29). Given Eq. (22), for each $\mathrm{i} \in \mathrm{I}_{\mathrm{a}}$, if α_i>0 then

$y_{i}-w^{T} \varphi\left(x_{i}\right)-b=\varepsilon+\xi_{i}$     (31)

and according to Eq. (19), if α_i<C₂ then

$\gamma_{i}>0$     (32)

Then, according to Eq. (24), we have:

$\xi_{i}=0$     (33)

Thus, for each $\mathrm{i} \in \mathrm{I}_{\mathrm{a}}$, if 0<α_i<C₂, the bias can be obtained from the following equation:

$\mathrm{b}=y_{\mathrm{i}}-\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)-\varepsilon$     (34)

According to Eq. (23), for each i $\mathrm{i} \in \mathrm{I}_{b}$, if α_i>0, then

$-y_{i}+w^{T} \varphi\left(x_{i}\right)+b=\varepsilon+\xi_{i}$     (35)

and according to Eq. (19), if α_i<C₂, then

$\gamma_{i}>0$     (36)

Then, according to Eq. (24), we have:

$\xi_{\mathrm{i}}=0$     (37)

Thus, for each $\mathrm{i} \in \mathrm{I}_{\mathrm{u}}$, if 0<α_i<C₂, the bias can be obtained from the following equation:

$\mathrm{b}=\varepsilon+\mathrm{y}_{\mathrm{i}}-\mathrm{w}^{\mathrm{T}} \varphi\left(\mathrm{x}_{\mathrm{i}}\right)$      (38)

According to Eq. (17), the optimal hyperplane or estimated function is as follows:

$f(x)=w^{T} \varphi(x)+b$$=\mathrm{b}+\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{a}}} \alpha_{\mathrm{i}} \mathrm{k}\left(\mathrm{x}, \mathrm{x}_{\mathrm{i}}\right)-\sum_{\mathrm{i} \in \mathrm{I}_{\mathrm{b}}} \alpha_{\mathrm{i}} \mathrm{k}\left(\mathrm{x}, \mathrm{x}_{\mathrm{i}}\right)$$=\mathrm{b}+\sum_{i=1}^{\mathrm{n}} \beta_{\mathrm{i}} \mathrm{k}\left(\mathrm{x}, \mathrm{x}_{\mathrm{i}}\right)=\mathrm{b}+\sum_{\mathrm{i} \in \mathrm{SV}} \beta_{\mathrm{i}} \mathrm{k}\left(\mathrm{x}, \mathrm{x}_{\mathrm{i}}\right)$     (39)

where, $\mathrm{SV}=\left\{\mathrm{i} \mid \beta_{\mathrm{i}} \neq 0\right\}$ is called support vector set.

The RLS method is one of the fastest methods for function estimation because it is an n + 1 variable problem with a close-form solution which is obtained by inverting an n × n matrix. The RLS has a high sensitivity to noisy data. Against this, ε-SVR has a good resistance against noise, but its runtime is not as good as the RLS runtime, because in order to solve the dual ε-SVR, we need to solve a quadratic programming problem with 2×n variables and linear constraints. ε-SVR supposes that the noise in the output variable is at most ε. Therefore, the center of a tube with radius ε, which is used as the estimated function, is determined in a way that the training data are located in that tube. Therefore, this method is robust to such noisy data. In this paper, to improve the runtime of ε-SVR, first, an initial estimated function was obtained using the RLS method. Then, unlike ε-SVR which uses all data to determine the lower and upper limits of the tube, our proposed method used the initial estimated function for determining the tube and the final estimated function. Strictly speaking, the data below and above the initial estimated function were used to estimate the upper and lower limits of the tube, respectively. Thus, the number of the model constraints and, consequently, the model runtime were reduced.

Our proposed method runtime is also lower than that of ε-TSVR and pair v-SVR because in each of ε-TSVR and pair v-SVR, two quadratic programming problems with n-variables and linear constraints are solved, while in our proposed method, an unconstrained quadratic programming problem with n+1 variables and a close form solution, and then, a quadratic programming problem with n-variables and linear constraints are solved. Solving an unconstrained quadratic programming problem is faster than a constrained quadratic programming problem of the same size. Our experiments on 15 benchmark dataset confirm that our proposed method is faster than ε-SVR, ε-TSVR and pair v-SVR, and its accuracy is comparable with the ε-SVR, ε-TSVR and pair v-SVR.