Improving the Solution of Least Squares Support Vector Machines with Application to a Blast Furnace System

Abstract

The solution of least squares support vector machines (LS-SVMs) is characterized by a specific linear system, that is, a saddle point system. Approaches for its numerical solutions such as conjugate methods Sykens and Vandewalle (1999) and null space methods Chu et al. (2005) have been proposed. To speed up the solution of LS-SVM, this paper employs the minimal residual (MINRES) method to solve the above saddle point system directly. Theoretical analysis indicates that the MINRES method is more efficient than the conjugate gradient method and the null space method for solving the saddle point system. Experiments on benchmark data sets show that compared with mainstream algorithms for LS-SVM, the proposed approach significantly reduces the training time and keeps comparable accuracy. To heel, the LS-SVM based on MINRES method is used to track a practical problem originated from blast furnace iron-making process: changing trend prediction of silicon content in hot metal. The MINRES method-based LS-SVM can effectively perform feature reduction and model selection simultaneously, so it is a practical tool for the silicon trend prediction task.

1. Introduction

As one kernel method, SVM works by embedding the input data 𝐱,𝐳∈𝑋 into a Hilbert space ℋ by a high-dimensional mapping Φ(⋅), and then trying to find a linear relation among the high-dimensional embedded data points [1, 2]. This process is implicitly performed by specifying a kernel function which satisfies 𝑘(𝐱,𝐳)=Φ(𝐱)^TΦ(𝐳), that is, the inner product of the embedded points. Given observed samples {𝐱_𝑖,𝑦_𝑖}^𝑛_𝑖=1 with size 𝑛, SVM formulates the learning problem as a variational problem of finding a decision function 𝑓 that minimizes the regularized risk functional [3, 4] min_𝑓∈ℋ𝑅𝑓=1𝑛^𝑛_𝑖=1𝑉𝑦_𝑖,𝑓𝐱_𝑖+𝜆𝑓²_ℋ, (1.1) where 𝑉(⋅,⋅) is called a loss function, 𝜆 is the so-called regularization parameter to trade off the empirical risk with the complexity of 𝑓, that is, ‖𝑓‖_ℋ, the norm in a reproducing kernel Hilbert space ℋ. By the representer theorem [3, 5], the optimal decision function 𝑓 satisfying (1.1) has the form 𝑓(𝐱)=^𝑛_𝑖=1𝛼_𝑖𝑘𝐱_𝑖,𝐱+𝑏, (1.2) where 𝛼_𝑖∈ℛ for 𝑖=1,…,𝑛, 𝑏∈ℛ. This equation can be easily used to tackle a practical problem if the kernel function is specified. To overcome the high computational complexity of traditional SVM, an interesting variant of the standard SVM, least squares support vector machines, has been proposed by Suykens and Vandewalle [6]. In the case of LS-SVM, the inequality constraints in ℓ₂ soft margin SVM are converted into equality constraints. The model training process of LS-SVM is performed by solving a specific linear equations, that is, a saddle point system which can be efficiently solved by iterative methods instead of a quadratic programming problem. Besides computational superiority extensive empirical studies have shown that LS-SVM is comparable to SVM in terms of generalization performance [7]; these features make LS-SVM an attractive algorithm and also a successful alternative to SVM. For the training of the LS-SVM, Van Gestel et al. [7] proposed to reformulate the 𝑛+1 order saddle point system into two 𝑛 order symmetric positive definite systems which can be solved in turn by the conjugate gradient (CG) algorithm. To speed up the training of LS-SVM, Chu et al. [8] employed the null space method to transform the saddle point system into a reduced 𝑛−1 order symmetric positive definite system which was solved with the CG algorithm also. The minimal residual (MINRES) method proposed by Paige and Saunders is a specialized method for solving a nonsingular symmetric system [9]. This method can avoid the LU factorization and does not suffer from break-down, so it is an efficient numerical method for solving symmetric but indefinite systems. The Karush-Kuhn-Tucker system of LS-SVM is a specified linear system, that is, a saddle point system. Considering the above point, to speed up the solution of LS-SVM model we employ the MINRES method to solve the linear system directly. The main contribution of this paper is to provide a potential alternative to the solution of LS-SVM model. Theoretical analysis of the three numerical algorithms for the solution of LS-SVM model indicates that the MINRES method is the optimal choice. Experiments on benchmark data sets show that compared with the CG method proposed by Suykens et al. and the null space method proposed by Chu et al., the MINRES solver significantly improves the computational efficiency and at the same time keeps almost the same generalized performance with the above two methods. To heel, the MINRES method-based LS-SVM model is constructed and further employed to identify blast furnace (BF) iron-making process, a complex nonlinear system. Practical application to a typical real BF indicates that the established MINRES method-based LS-SVM model is a good candidate to predict the changing trend of the silicon content in BF hot metal with low time cost. The possible application of this work is to aid the BF operators to judge the inner state of BF getting hot or chilling in time properly, which can provide a guide for them to determine the direction of controlling BF in advance. The rest of this paper is organized as follows. In Section 2, we give a review for LS-SVM. Section 3 presents three numerical solutions for LS-SVM. It is followed by extensive experimental validations of the proposed method in Section 4. Section 5 concludes the paper and points out the possible future research.

2. Formulation of LS-SVM

The primal problem of LS-SVM can be formulated following unified format: min_{𝐰,𝑏,𝐞}12𝐰^T𝐶𝐰+2^𝑛_𝑖=1𝑒²_𝑖s.t.𝑦_𝑖=𝐰^TΦ𝐱_𝑖+𝑏+𝑒_𝑖,𝑖=1,…,𝑛, (2.1) for both regression analysis and pattern classification. In (2.1) 𝑛 is the total number of training samples, 𝐱_𝑖 is the 𝑖th input vector, 𝑦_𝑖 is the 𝑖th output value/label for regression/classification problem, 𝑒_𝑖 is the 𝑖th error variable, 𝐶>0 is the regularization parameter, and 𝑏 is the bias term. The Lagrangian of (2.1) is given below: 𝐿()=1𝐰,𝑏,𝐞;𝜶2𝐰^T𝐶𝐰+2^𝑛_𝑖=1𝑒²_𝑖−^𝑛_𝑖=1𝛼_𝑖𝐰^TΦ𝐱_𝑖+𝑏+𝑒_𝑖−𝑦_𝑖, (2.2) where 𝛼_𝑖 is the 𝑖th Lagrange multiplier. For the convex program (2.1), it is obvious that the Slater constraint qualification holds. Therefore, the optimal solution of (2.1) satisfies its Karush-Kuhn-Tucker system 𝜕𝐿𝜕𝐰=0→𝐰=^𝑛_𝑖=1𝛼_𝑖Φ𝐱_𝑖,𝜕𝐿𝜕𝑏=0→^𝑛_𝑖=1𝛼_𝑖=0,𝜕𝐿𝜕𝑒_𝑖=0→𝛼_𝑖=𝐶𝑒_𝑖,𝑖=1,…,𝑛,𝜕𝐿𝜕𝛼_𝑖=0→𝑦_𝑖=𝐰^TΦ𝐱_𝑖+𝑏+𝑒_𝑖,𝑖=1,…,𝑛. (2.3) After eliminating variables 𝐰 and 𝐞 the Karush-Kuhn-Tucker system (2.3) can be reformulated following saddle point system [10]: 1𝐾+𝐶𝐼𝟏_𝑛𝟏^T_𝑛0𝜶𝑏=𝐲0, (2.4) where 𝐾_𝑖𝑗∶=𝑘(𝐱_𝑖,𝐱_𝑗)=Φ(𝐱_𝑖)^TΦ(𝐱_𝑗), 𝐼 stands for unit matrix, 𝟏_𝑛 denotes an 𝑛-dimensional vector of all ones, and 𝐲=(𝑦₁,…,𝑦_𝑛)^T.

3. Solution of LS-SVM

In this section, we give a brief review and some analysis of the three mentioned numerical algorithms for solution of LS-SVM.

3.1. Conjugate Gradient Methods

The kernel matrix 𝐾 is a symmetric positive semidefinite matrix and the diagonal term 1/𝐶 is positive, so the matrix 𝐻∶=𝐾+(1/𝐶)𝐼 is symmetric and positive definite. Through the following matrix transformation 𝑀^T1𝐾+𝐶𝐼𝟏_𝑛𝟏^T_𝑛0𝑀𝑀⁻¹𝜶𝑏=𝑀^T𝐲0, (3.1) where 𝑀=𝐼𝐻⁻¹𝟏_𝑛01, (3.2) the saddle point system (2.4) can be factorized into a positive definite system [11] 𝐻00𝟏^T_𝑛𝐻⁻¹𝟏_𝑛𝜶+𝐻⁻¹𝟏_𝑛𝑏𝑏=𝐲𝟏^T_𝑛𝐻⁻¹𝐲. (3.3)

Suykens et al. suggested the use of the CG method for the solution of (3.3) and proposed to solve two 𝑛 order positive definite systems. More exactly, their algorithm can be described as follows.

Step 1. Employ the CG algorithm to solve the linear equations 𝐻𝜼=𝟏_𝑛 and get the intermediate variable 𝜼.

Step 2. Solve the intermediate variable 𝝁 from 𝐻𝝁=𝐲 by the CG method.

Step 3. Obtain Lagrange dual variables 𝜶=𝝁−𝑏𝜼 and bias term 𝑏=𝟏^T_𝑛𝝁/1^T_𝑛𝜼.
The output of any new data 𝐱 can subsequently be deduced by computing the decision function 𝑓(𝐱)=𝐰^TΦ(𝐱)+𝑏=∑^𝑛_𝑖=1𝛼_𝑖𝑘(𝐱_𝑖,𝐱)+𝑏.

3.2. Null Space Methods

In what was mentioned previously, to get the intermediate variable 𝜼 and 𝝁 two 𝑛 order positive definite systems need to be solved by CG methods. Chu et al. [8] proposed an interesting method to the numerical solution of LS-SVM by solving one 𝑛−1 order reduced system of linear equations. The improved method suggested by Chu et al. can be seen as one kind of null space method. The saddle point system (2.4) can be written as 𝐻𝜶+𝟏_𝑛𝑏=𝐲,𝟏_𝑛𝜶=0. (3.4) Chu et al. specified a particular solution of 𝟏_𝑛𝜶=0 as 𝜶=𝟎 and the null space of 𝟏_𝑛𝜶=0 as 𝐼𝑍=^𝑛−1−𝟏^T_𝑛−1. (3.5) Through solving the following reduced system of order 𝑛−1 for the auxiliary unknown 𝝂, 𝑍^T𝐻Z𝝂=𝑍^T𝐲, (3.6) the solution of the saddle point system (2.4) can be obtained as 𝜶=𝑍𝝂 and 𝑏=(1/𝑛)𝟏^T_𝑛(𝐲−𝐻𝜶).

3.3. Minimal Residual Methods

The vector sequences in the CG method correspond to a factorization of a tridiagonal matrix similar to the coefficient matrix. Therefore, a breakdown of the algorithm can occur corresponding to a zero pivot if the matrix is indefinite. Furthermore, for indefinite matrices the minimization property of the CG method is no longer well defined. The MINRES method proposed by Paige and Saunders [9] is a variant of the CG method that avoids the LU factorization and does not suffer from breakdown. It minimizes the residual in the ℓ₂-norm which is an efficient numerical algorithm for solving symmetric but indefinite systems; the corresponding convergence behavior of the MINRES method for indefinite systems has been analyzed by Van der Vorst [12]. The purpose of this paper is to employ the MINRES method to solve the saddle point system (2.4) directly. Next we gave a brief review of the MINRES algorithm. Let 𝐱₀ be an initial guess for the solution of the symmetric indefinite linear system 𝐴𝐱=𝐛. One can obtain the iterative sequence 𝐱_𝑚,  𝑚=1,2,… such that 𝐫_𝑚2=𝐛−𝐴𝐱_𝑚2=min_{𝐱∈𝐱0+𝒦𝑚(𝐴,𝐫0)}‖𝐛−𝐴𝐱‖₂, (3.7) where 𝐫_𝑚=𝐛−𝐴𝐱_𝑚 is the 𝑚th residual for 𝑚=1,2,…, and 𝒦_𝑚𝐴,𝐫₀=span𝐫₀,𝐴𝐫₀,…,𝐴^𝑚−1𝐫₀ (3.8) is the 𝑚th Krylov subspace. Lanczos methods can be used to generate an orthonormal basis of 𝒦_𝑚(𝐴,𝐫₀), and then only two basis vectors are needed to compute 𝐱_𝑚; see, for example, [12]. The detailed implementation of the MINRES algorithm can be found in [12].

It has been shown that rounding errors are propagated to the approximate solution with a factor proportional to the square of the condition number of coefficient matrix [12]; one should be careful with the MINRES method for ill-conditioned systems.

3.4. Some Analysis on These Three Numerical Algorithms

The properties of short recurrences and optimization [12] make the CG method the first choice for the solution of a symmetric positive definite system. Suykens et al. transformed the 𝑛+1 order saddle point system (2.4) into two 𝑛 order positive definite systems which are solved by CG methods. However, it is time consuming to solve two 𝑛 order positive definite systems with large scales. To overcome this shortcoming, Chu et al. [8] transformed equivalently the original 𝑛+1 order system into an 𝑛 order symmetric positive definite system, and then the CG method can be used. This method can be seen as a null space method. Unfortunately, the transformation may destroy heavily the sparse structure and increase greatly the condition number of the original system. This can hugely slow down the convergence rate of the CG algorithm. Theoretical analysis about the influence of the transformation on the condition number is indispensable, but it is rather difficult. We leave it as an open problem. In this paper, the MINRES method is directly applied to solve the original saddle point problem of 𝑛+1 order. Similar to the CG method, the MINRES method also has properties of short recurrences and optimization.

In light of the analysis mentioned above, the MINRES method should be the first choice for the solution of LS-SVM model, since it avoids solving two linear systems and destroying the sparse structure of the original saddle point system simultaneously.

4. Numerical Implementations

4.1. Experiments on Benchmark Data Sets

In this section we give the experimental test results on the accuracy and efficiency of our method. For comparison purpose, we implement the CG method proposed by Suykens and Vandewalle [6] and the null space method suggested by Chu et al. [8]. All experiments are implemented with MATLAB version 7.8 programming environment running on an IBM compatible PC under Window XP operating system, which is configured with Intel Core 2.1 Ghz CPU and 2 G RAM. The generalized used Gaussian RBF kernel 𝑘(𝐱,𝐳)=exp(−‖𝐱−𝐳‖²/𝜎²) is selected as the kernel function. We use the default setting for kernel width 𝜎², that is, set kernel width as the dimension of inputs.

We first compare three algorithms on three benchmark data sets: Boston, Concrete, and Abalone, which are download from UCI [13]. Each data set is randomly partitioned into 70% training and 30% test sets. We also list the condition numbers of coefficients matrices solved by three methods for the analysis of the computing efficiencies. As shown in Tables 1–3 the condition number for the CG method is the least one and the condition number for the null space method significantly increases.

Table 1 

Experimental results of three methods on Boston data set.

Table 2 

Experimental results of three methods on Concrete data set.

Table 3 

Experimental results of three methods on Abalone data set.

The columns of Cond in Tables 1, 2, and 3 show that compared with the CG method the condition number for the MINRES method increases a bit, but much less than the condition number of the null space method. The orders of linear equations solved by the CG method, the null space method, and the MINRES method are 𝑛, 𝑛−1, and 𝑛+1, respectively. The condition numbers for the CG method and the MINRES method are very close, but we have to solve two systems of 𝑛−1 order using CG methods. Hence, the running time of the MINRES method should be less than that of the CG method. CPU column in Tables 1–3 shows that the MINRES method-based LS-SVM model costs much less running time than the CG method and the null space method-based LS-SVM model in all cases of setting C. So the MINRES method-based LS-SVM model is a preferable algorithm for solving LS-SVM model. In the next subsection, we will employ the MINRES method-based LS-SVM model to solve a practical problem.

4.2. Application on Blast Furnace System

Blast furnace, one kind of metallurgical reactor used for producing pig iron, is often called hot metal. The chemical reactions and heat transport phenomena take place throughout the furnace as the solid materials move downwards and hot combustion gases flow upwards. The main principle involved in the BF iron-making process is the thermochemical reduction of iron oxide ore by carbon monoxide. During the iron-making period, a great deal of heat energy is produced which can heat up the BF temperature approaching 2000°C. The end products consisting of slag and hot metal sink to the bottom and are tapped periodically for the subsequent refining. It will take about 6–8 h for a cycle of iron-making [11]. BF iron-making process is a highly complex nonlinear process with the characteristics of high temperature, high pressure, concurrence of transport phenomena, and chemical reactions. The complexity of the BF and the occurrence of a variety of process disturbances have been obstacles for the adoption of modeling and control in the process. Generally speaking, to control a BF system often means to control the hot metal temperature and components, such as silicon content, sulfur content in hot metal, and carbon content in hot metal within acceptable bounds. Among these indicators, the silicon content often acts as a chief indicator to represent the thermal state of the BF, an increasing silicon content meaning a heating of the BF while a decreasing silicon content indicating a cooling of the BF [11, 14]. Thus, the silicon content is a reliable measure of the thermal state of the BF, and it becomes a key stage to predict the silicon content for regulating the thermal state of the BF. Therefore, it has been the active research issue to build silicon prediction model in the recent decades, including numerical prediction models [15] and trend prediction models [11].

In this subsection, the tendency prediction of silicon content in hot metal is transformed as a binary classification problem. Samples with increasing silicon content are denoted by +1 whereas a decreasing silicon content is denoted by −1. In the present work, the experimental data is collected from a medium-sized BF with the inner volume of about 2500 m³. The variables closely related to the silicon content are measured as the candidate inputs for modeling. Table 4 presents the variables information from the studied BF. There are totally 801 data points collected with the first 601 points as train set and the residual 200 points as testing set. The sampling interval is about 1.5 h for the current BF. Figure 1 illustrates the evolution of the silicon content in hot metal.

Table 4 

A list of input variables.

Figure 1 

Evolution of silicon content in hot metal.

There are in total 15 candidate variables listed in Table 4 from which to select model inputs. Generally, too many input parameters will increase the complexity of model while too little inputs will reduce the accuracy of model. A tradeoff has to be taken between the model complexity and accuracy when selecting the inputs. Therefore, it is necessary to screen out less important variables as inputs from these 15 candidate variables. Here, the inputs are screened out by an integrative way that combines 𝐹-score method [16] for variables ranking and cross-validation method for variables and model parameters selection.

𝐹-score is an effective tool for feature selection in data mining and can give feature ranking by evaluating the discrimination of two sets with real values. For those 15 candidate variables in Table 4, their 𝐹-scores are defined as follows: 𝐹_𝑠(𝑖)=1𝑁₊−1^𝑁+_𝑗=1𝑥^(𝑖)_𝑗,+−𝑥₊^(𝑖)2+1𝑁₋−1^𝑁−_𝑗=1𝑥^(𝑖)_𝑗,−−𝑥₋^(𝑖)2−1⋅𝑥₊^(𝑖)−𝑥^(𝑖)2+𝑥₋^(𝑖)−𝑥^(𝑖)2,𝑖=1,…,𝑑, (4.1) where 𝑥^(𝑖), 𝑥₊^(𝑖) and 𝑥₋^(𝑖) stand for the mean of the 𝑖th attribute of the whole training, positive and negative examples, respectively, while 𝑥^(𝑖)_𝑗,+ and 𝑥^(𝑖)_𝑗,− are the 𝑖th variable of the 𝑗th positive and negative instance, respectively. Hence, a variable ranking can be achieved through 𝐹-score method. Table 4 gives the results of 𝐹-scores of all 15 variables, which are ranked according to the 𝐹-score values. As one kernel-based learning model, the kernel parameter 𝜎², and regularized parameter 𝐶 play an important role in LS-SVM, so one should pay attention to selecting proper parameters. Grid search-based ten-fold cross-validation is executed on the train set for searching the optimal (𝜎²,C). The searching grid for model parameters is set as 2⁻⁵,2⁻⁴,…,2¹⁰×2⁻⁵,2⁻⁴,…,2¹⁰. (4.2)

Mean accuracy in Table 4 stands for the average accuracy under ten-fold cross-validation experiments of LS-SVM model on some grid points with the best performance. In the current work, we first select the variable with highest 𝐹-score as model input and then add variables one by one according to their 𝐹-scores. Mean accuracy under all kinds of input variables can be achieved and the results are shown in Table 4. The following are shown by the mean accuracy column: (1) at the beginning, the mean accuracy increases gradually as more candidate variables are taken as model inputs; (2) the largest mean accuracy appears when CO₂ is included within the input set; (3) when the mean accuracy is beyond the maximum, it will fluctuate as the residual variables are added by turns into the input set. These results indicate that, as the studied BF is concerned, the optimal input set is [Si, S, BI, FS, BV, CO₂] with the model parameters setting (𝜎²,𝐶)=(2⁹,2⁸). Table 5 lists the LS-SVM model accuracy including with/without feature and model selection versions on testing set. In the case of without feature and model selection version, all candidate variables are selected as inputs, and we use the default setting for LS-SVM model; that is, set kernel width 𝜎² equal to the dimension of input variable and set regularized parameter 𝐶 as 1. The information in the second row of this table, such as 34/42, denotes that there are 42 times predicted results that are ascending trend, and 34 times predictions are successful. The confidence level of the LS-SVM model without model and feature selection fluctuates severely between the ascending and descending prediction from 80.95% to 58.86%. The difference of confidence levels of LS-SVM model with model and feature selection between ascending and descending prediction is reduced to 2.19% indicating that model and feature selection procedure enhances the stability of the LS-SVM model obviously. As the last column of Table 5 shows, TSA of LS-SVM model with feature and model selection procedure is significantly improved compared with LS-SVM model without feature and model selection, so the selection procedure is indispensable for the current practical application. Table 6 lists the running time of three mentioned numerical algorithms when performing feature and model selection procedure. The cost time of the MINRES method is reduced significantly compared with the other algorithms. In a word, the feature and model selection procedure can be effectively performed for the MINRES method-based LS-SVM, and it is meaningful for practical using.

Table 5 

Predictive results of LS-SVM model with/without feature and model selection.

Table 6 

Running time of three numerical methods on model identification.

5. Conclusions and Points of Possible Future Research

In this paper, we have proposed an alternative, that is, the MINRES method, to the solution of LS-SVM model which is formulated as a saddle point system. Numerical experiments on UCI benchmark data sets show that the proposed numerical solution method of LS-SVM model is more efficient than the algorithms proposed by Suykens and Vandewalle [6] and Chu et al. [8]. To heel, the MINRES method-based LS-SVM model including feature selection from extensive candidate and model parameter selection is proposed and employed for the silicon content trend prediction task. The practical application to a typical real BF indicates that the proposed MINRES method-based LS-SVM model is a good candidate to predict the trend of silicon content in BF hot metal with low running time.

However, it should be pointed out that despite the MINRES method-based LS-SVM model displaying low running time, lack of metallurgical information may be the root to the limited accuracy of the current prediction model. So there is much work worth investigating in the future to further improve the model accuracy and increase the model transparency, such as constructing predictive model by integrating domain knowledge and extracting rules. The extracted rules can account for the output results with detailed and definite inputs information, which may further serve for the control purpose by linking the output results with controlled variables. These investigations are deemed to be helpful to further improve the efficiency of predictive model.

Acknowledgment

This work was partially supported by National Natural Science Foundation of China under Grant no. 11126084, Natural Science Foundation of Shandong Province under Grant no. ZR2011AQ003, Fundamental Research Funds for the Central Universities under Grant no. 12CX04082A, and Public Benefit Technologies R&D Program of Science and Technology Department of Zhejiang Province under Grant No. 2011C31G2010136.