Hai-Xiang Zhao
email: haixiang.zhao@ecp.fr

Frédéric Magoulès
email: frederic.magoules@hotmail.com

Parallel Support Vector Machines on Multi-core and Multiprocessor Systems

Keywords: support vector machine, parallel, multi-core, Map-Reduce

This paper proposes a new and efficient parallel implementation of support vector machines based on decomposition method for handling large scale datasets. The parallelizing is performed on the most time-and-memory consuming work of training, i.e., to update the vector f . The inner problems are dealt by sequential minimal optimization solver.

Since the underlying parallelism is realized by the shared memory version of Map-Reduce paradigm, our system is easy to build and particularly suitable to apply to multi-core and multiprocessor systems. Experimental results show that on most of the tested datasets, our system offers higher than fourfold increase in speed compared to Libsvm, and it is also far more efficient than the MPI implementation Pisvm.

I. INTRODUCTION

Support vector machine (SVM) is a popular supervised learning method of solving classification and regression problems [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. It shows robust generalization ability in numerous applications such as image processing, text mining, neural analysis and energy efficiency modeling [START_REF] Lai | Vapnik's learning theory applied to energy consumption forecasts in residential buildings[END_REF], [START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF]. The training of SVM is essentially a quadratic optimization problem which is both time and memory costly while running on computers, making it a challenge to apply SVM on large scale problems. Several optimizing or heuristic methods have been proposed to accelerate the training and reduce the memory occupation, such as shrinking, chunking [START_REF] Osuna | Training support vector machines: an application to face detection[END_REF], kernel caching [START_REF] Joachims | Making large-scale support vector machine learning practical[END_REF], approximation of kernel matrix [START_REF] Chang | Psvm: Parallelizing support vector machines on distributed computers[END_REF]. In addition, certain scalable solvers can be used such as Sequential Minimal Optimization (SMO) [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF], mixture SVMs [START_REF] Collobert | A parallel mixture of svms for very large scale problems[END_REF], primal estimated sub-gradient solver [START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF]. Despite these efforts, however, a more sophisticated and satisfactory solution is still expected for this challenging research problem.

Thanks to the modern chip manufacturing, we are entering the multi-core era. Computers with multi-cores or multiprocessors are becoming more available and affordable. This paper aims to investigate and demonstrate how SVM -a popular machine learning algorithm can benefit from this modern platform. A new parallel SVM that is particularly suitable to shared memory system is proposed. Decomposition method, caching and SMO inner quadratic problem solver are composed in the implementation as the key techniques. For the purpose of achieving easy implementation without sacrificing performance, the state-of-the-art parallel programming framework Map-Reduce is chosen to perform the underlying parallelism. The system is therefore called MRPsvm, which stands for "Map-Reduce parallel SVM". Comparative system analysis and experimental results on benchmark datasets show significant memory saving and overwhelming speed increase in our system.

The following sections are organized as follows. Section II introduces the basic problem of SVM training and how the decomposition method can solve the problem. Section III explains the some key points of system implementation. Section IV describes the related work. Section V presents the numerical experiments on benchmark datasets and corresponding results. Conclusions are drawn in section VI.

II. SUPPORT VECTOR MACHINE

SVM for classification purpose aims at finding a hyperplane to separate two classes with maximum margin. Let x i denotes the ith training sample, y i denotes the corresponding label with the value either -1 or 1, i = 1, 2, ..., l. l is the total number of training samples. The dual form of the SVM can be written as the following convex quadratic function.

min α 1 2 α T Qα - l i=1 α i (1)
subject to

y T α = 0 (2) 0 ≤ α i ≤ C ∀i = 1, 2, ..., l (3)
Where Q is a l by l positive semi-definite matrix. Each element of Q has the form Q ij = y i y j K(x i , x j ) where K(x i , x j ) is called kernel function which can be substituted by any Mercer kernel. C > 0 is the upper bound used to trade off between classifier performance on training data and its generalization ability. α is a vector of l variables, where each element α i is the weight of the corresponding training sample (x i , y i ). The objective of the problem is to find the solution of α which makes (1) minimized and the constraints (2) (3) fulfilled. After we find out the optimal α, we can use the following decision function to predict the labels with the new input of x:

sgn l i=1 y i α i K(x i , x) -b
where b is a constant value which can be easily calculated in the training step, and

sgn(u) = -1, for u < 0 1, for u > 0 A. Decomposition method
Since the size of the kernel matrix Q is the square of the number of samples, i.e., l 2 , so that it is difficult to store the whole matrix in memory when l is very big. Osuna et al. [START_REF] Osuna | Training support vector machines: an application to face detection[END_REF] proposed a method to decompose the problem into smaller tasks. In each task, certain parts of α are chosen to be optimized, while the rest of α remains in constant value. The selected part is called working set. Then the program repeats the select-optimize process until global optimality conditions are satisfied. Let B denote the working set with n variables and N denote the non-working set with (ln) variables. Then, α, y and Q can be correspondingly written as:

α = α B α N , y = y B y N , Q = Q BB Q BN Q N B Q N N
Accordingly, the small task can be written as:

min 1 2 α T B Q BB α B -α T B (1 -Q BN α N ) + 1 2 α T N Q N N α N -α T N (4) subject to α T B y B + α T N y N = 0 (5) 0 ≤ α B ≤ C (6)
Since 4) remains constant in an iteration, this term can be omitted while calculating. So that function (4) has essentially the same form as function [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. The advantage of this problem decomposition method is that the newly generated task is small enough to be solved by most off-the-shelf methods. In our work, we choose SMO as the inner small task solver due to its relative simplicity, yet high performance characteristics. Interestingly, SMO is itself an extreme case of the decomposition method where the working set contains only two variables. This kind of binary sub-problem can be easily solved analytically [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF], [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. As stated in [START_REF] Joachims | Making large-scale support vector machine learning practical[END_REF], the solution of the sub-problem (4) is strictly feasible towards the optimum solution of global problem [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. This feature guarantees the global convergence of this decomposition method.

( 1 2 α T N Q N N α N -α T N ) in (

B. Optimality condition

When to stop the optimizing process is evaluated by the Karush-Kuhn-Tucker condition iteratively. The evaluating procedure can be summarized as follows. Firstly, we classify the training samples into two categories:

I up (α) = {t|α t < C, y t = 1 or α t > 0, y t = -1} I low (α) = {t|α t < C, y t = -1 or α t > 0, y t = 1}
Secondly we search two extreme values m(α) and M (α) by m(α) = max i∈Iup(α)

-y i f i M (α) = min i∈I low (α) -y i f i where f i = l j=1 α j Q ij -1
is the gradient of (1). Since f i should be calculated for all l variables and Q ij may not be stored in the memory, the calculation of f i becomes the main task of the whole optimizing work.

After we find out m(α) and M (α), we define the stop condition as:

m(α) -M (α) ≤ ǫ (7)
where ǫ is a pre-defined value.

C. Working set selection

The selection of working set B directly influences the speed of convergence. For inner SMO solver, the maximal violating pair is selected to be the binary working set according to the second order information [START_REF] Fan | Working set selection using second order information for training support vector machines[END_REF]. We do not state here how inner SMO solver works since it has been discussed in detail in [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF] [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. For the selection of working set B, we simply consider the first order information and select, in some sense, the maximal violating pairs as proposed by Zanni et al. [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF]. Suppose the required size of B is n, we choose q (q < n) variables from α by sequentially selecting pair of variables which satisfy [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF]. The remaining (nq) variables are chosen as those who entered B in the last iteration but not yet selected in current B. The selection of these (n-q) variables follows the sequence: free variables which satisfy 0 < α i < C, lower bound variables which are equal to 0, upper bound variables which have the value C. The reason for putting restraint on the number of new variables entering the working set is to avoid frequent entering-leaving of certain variables. Otherwise, the speed of convergence would considerably slow down [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF].

After the working set is optimized, f is updated by:

f * i = f i + j∈B ∆α j Q ij i = 1, 2, ..., l (8)
where ∆α j is the newly optimized α j minus the old α j . This procedure is crucial as it prepares f for the next iteration to do optimality condition evaluation and working set selection. In fact, this is the most time-consuming step in SVM training due to the computation of Q ij . Therefore, the main parallelizing work of MRPsvm is based on the updating of f . The whole parallel decomposition method is briefly outlined in algorithm 1.

Algorithm 1 Parallel decomposition solver INPUT: data set (x i , y i ), ∀i ∈ {1, ..., l}, n, ǫ

INITIALIZE: α i = 0, f i = -1, ∀i ∈ {1, ..., l} CALCULATE: I up , I low , m(α), M (α) REPEAT select working set B, | B |= n optimize α i , ∀i ∈ B, by SMO solver update f i , ∀i ∈ {1, ..., l} in parallel calculate I up , I low , m(α), M (α) UNTIL m(α) -M (α) ≤ ǫ III. SYSTEM IMPLEMENTATION
In this section, we introduce some key points of our system implementation, i.e., why and how Map-Reduce is used to do the parallelism, caching technique and data representation.

A. Map-Reduce on multi-core system

Map-Reduce is a parallel programming framework originally proposed by Google [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF]. It can help us extract parallelism of computations on large data sets by taking advantage of distributed systems. It allows users to write code in a functional style: map computations on separated data, generate intermediate key/value pairs and then reduce the summation of intermediate values assigned to the same key. The runtime system automatically handles low-level mapping, scheduling, parallel processing and fault tolerance. It is a simple but very useful framework.

Except the plenty researches of Map-Reduce in distributed environments, such as [START_REF] Zaharia | Improving mapreduce performance in heterogeneous environments[END_REF]- [START_REF] Yang | simplified relational data processing on large clusters[END_REF], Chu et al. [START_REF] Chu | Map-reduce for machine learning on multicore[END_REF] has attempted to use this technology to develop a general programming framework on multi-core systems for machine learning applications. For SVM, they have parallelized linear SVM with primal problem. To conduct further development, our work attempts to develop a parallel SVM for general classification problem, and the solved SVM problem is in dual form.

According to the mechanism of Map-Reduce, this framework is naturally suitable to deal with our problem (8) in parallel. Phoenix designed by Ranger et al. [START_REF] Ranger | Evaluating mapreduce for multi-core and multiprocessor systems[END_REF] implements Map-Reduce on shared memory systems. It supplies a common API for users to easily parallelize their applications without conducting concurrency management. The map and reduce tasks are performed in threads. An efficient integrated runtime system is supposed to handle the parallelism, resource management and fault recovery by itself. Their experimental results show that this system works just as efficient as pure Pthreads implementation. Therefore, we choose this system as the underlying Map-Reduce handler to implement our MRPsvm.

The parallel computation works as follows. Variables in working set B is divided into several parts, the calculation of f * is also divided into several parts in the same manner as for B. Each part is then assigned to a map process. The major work of each map process is thus to calculate Q ij if they are not saved in the memory, and then calculate j ∆α j Q ij for all i = 1, 2, ..., l. After the distributed calculations of these maps, final f * is summed up by the reduction. Figure 1 explicitly depicts the detailed problem and the partitioning manner. It is obvious that the problem is divided on kernel columns. (j k , k = 1, ..., n) are the variable index of working set in kernel matrix, which gives the kth variable in B with its index in Q as j k . In practice, it is not necessary to update f on all the n variables since some of them are so marginal that can be omitted.

B. Caching technique

As stated in the previous sections, for large scale problems, kernel matrix Q is too large to be stored in memory, and the calculation of kernel elements Q ij is the dominant work that slows down the training. It is an effective technique to cache the kernel elements in memory as much as possible. MRPsvm maintains a fix-sized cache which stores recently accessed or generated kernel columns. The cache replacement policy is a simple least-recent-use strategy, as same as that of Libsvm. Only the column currently needed but not hit in the cache will be calculated. All parallel maps share an unique copy of cache in the shared memory. In consequence, the operation of inserting a new column into the cache performed by whichever map should be synchronized.

For inner SMO solver, the kernel matrix size is dependent on the size of working set B which is normally set as 1024 according to the knowledge of experience, it is practical to pre-compute and cache the full version of this small kernel matrix.

C. Sparse data representation

To reduce the storage requirements, the sample vectors x i are stored by sparse representation. When calculating a kernel column (Q ij , i = 1, 2, ..., l), we need to unroll the jth sample vector to dense format and then calculate the dot products of this vector with the other l sample vectors.

IV. RELATED WORK

Since the essential convex quadratic problem can be solved by several methods, different kinds of parallel SVMs

1 2 2 * 1 1 1 1 1 1 1 k k m m n n j j j j j j j j j j f f Q Q Q Q Q D D D D D ' ' ' ' ' "" "" 1 1 2 2 * 2 2 2 2 2 2 2 k k m m n n j j j j j j j j j j f f Q Q Q Q Q D D D D D ' ' ' ' ' ""
""

1 1 2 2 * k k m m n n i i j ij j ij j ij j ij j ij f f Q Q Q Q Q D D D D D ' ' ' ' ' ""
"" were proposed according to the particular quadratic problem solver. Based on stochastic gradient descent method, P-packSVM optimizes SVM training directly on the primal form for arbitrary kernels [START_REF] Zhu | P-packsvm: Parallel primal gradient descent kernel svm[END_REF]. Very high efficiency and competitive accuracy have been achieved by the parallel implementation. Psvm proposed in [START_REF] Chang | Psvm: Parallelizing support vector machines on distributed computers[END_REF] is based on interior point solver. It approximates the kernel matrix by incomplete Cholesky factorization. Memory requirement is reduced and scalable performance has been achieved. Bickson et al. [START_REF] Bickson | A gaussian belief propagation solver for large scale support vector machines[END_REF] solve the problem by Gaussian belief propagation, a method from complex system domain. The parallel solver brings competitive speedup on large scale problems.

1 1 2 2 * k k m m n n l l j lj j lj j lj j lj j lj f f Q Q Q Q Q D D D D D ' ' ' ' ' "" "" … … … … … … … map.1 map.r
The decomposition method attracts more attentions than the above solvers. Graf et al. [START_REF] Graf | Parallel support vector machines: The cascade svm[END_REF] train several SVMs on small data partitions, then they aggregate support vectors from two pair SVMs to form new training samples on which another training is performed. The aggregation is repeated until only one SVM remains. The similar idea is used by Dong et al. [START_REF] Dong | Fast svm training algorithm with decomposition on very large data sets[END_REF], in this work, sub-SVMs are performed on block diagonal matrices which are the approximation to the original kernel matrix. Consequently, nonsupport vectors are removed when dealing these sub-problems. Zanni et al. [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF] parallelize SVM-light with improved working set selection and inner quadratic problem solver. Hazan et al. [START_REF] Hazan | A parallel decomposition solver for svm: Distributed dual ascend using fenchel duality[END_REF] propose a parallel decomposition solver using Fenchel Duality.

Similar to our implementation, Cao et al. [START_REF] Cao | Parallel sequential minimal optimization for the training of support vector machines[END_REF] and Catanzaro et al. [START_REF] Catanzaro | Fast support vector machine training and classification on graphics processors[END_REF] parallelized SMO solver for training SVM. Both work mainly focuses on updating gradient for Karush-Kuhn-Tucker condition evaluation and the working set selection. The difference between them is on the implementation details and the programming models. Specifically speaking, the first work is conducted by using MPI on clusters while the second one by Map-Reduce threads on modern GPU platform. In our work, we also adopt SMO algorithm. But we use it as the inner quadratic problem solver without any parallel computation, instead, we perform the parallelization on external decomposition procedure. The main advantage of our coarse-grained parallelism is that it can significantly reduce the burden of overheads since the number of iterations in global decomposition procedure (where n ≫ 2) is extremely smaller than that of pure SMO algorithm (where n = 2).

Pisvm [START_REF] Brugger | Parallel support vector machines[END_REF] uses the same decomposition method as ours to train SVM in parallel. But its implementation is different from our MRPsvm in many aspects. First, it is based on MPI implementation and aims at extracting parallelism from distributed memory systems, while our parallel algorithm is conducted by Map-Reduce threads on shared memory system. They are based on totally different models, and the number of threads, the granularity are different as well. Second, in its implementation, each process stores a copy of data samples. In the contrary, MRPsvm stores only one copy in the shared memory. Third, Pisvm adopts a distributed cache strategy in order to share the saved kernel elements across all of the processes. Each process stores locally a piece of the cache. Consequently, the work for updating gradients is divided and assigned globally to proper processors according to the cache locality. In contrast, our MRPsvm has only one copy of the cache, and each processor accesses the cache equally, so that the overhead of global assignment is avoided. However, we have to note that synchronization on cache write is required. In next section, we will compare the performance of Pisvm with that of MRPsvm on some benchmark datasets, providing direct evidence that our system is more efficient and suitable than the MPI implementation on multi-core systems.

V. EXPERIMENTS AND RESULTS

We test MRPsvm by comparing it with the parallel implementation Pisvm and the serial implementation Libsvm on five widely used benchmark datasets. Although this comparison may not based on systems especially designed for multi-core architecture, we still have good reasons for doing so. Firstly, as the best knowledge as we know, there is no existing parallel implementation of general SVM that is specially developed for multi-core systems. Therefore there is a strong need to verify if our system could outperform the state-of-the-art parallel implementation. Secondly, most of the systems surveyed in section IV are not available to the public, while Pisvm , as a typical parallel implementation of SVM, is easy to obtain. Thirdly, the quadratic problem solver of Pisvm is the same as MRPsvm, hence, if we compare our system with Pisvm, the advantage of Map-Reduce framework is more convincing. Two computers with different hardware architectures are adopted to check hardware effects. As shown in table I, the first computer has 4 cores with a shared L2 cache and memory. The second one is a dual-processor system with 4 cores in each processor. The cores in the same processor share one cache, and the main memory is shared among all of the cores. Both of the two computers are running Linux 2.6.27-7.

The five datasets are shown in table II. They vary in sample size and dimension. We train all SVMs with Gaussian kernel. The tolerance of termination criterion is set to 0.01. Since we focus on comparing of three systems, the outputs of these classifiers may not be optimal. In other words, we do not guarantee the parameters of SVM, i.e., C and γ, to reach optimal values. They are just chosen from the literature as shown in the last two columns of table II.

Since the caching technique is crucial for performance, for a reliable comparison, we set the cache size to be the same for all three systems. Furthermore, we restrict the cache size to be far smaller than the memory size in order to minimize page faults in runtime. Here we have to emphasize that the following reported performance might not be optimal for all three systems, only serving for comparison purpose.

Table III shows the results of the three implementations performed on 4 processors. The time columns represent the whole training time, i.e., from reading the problem to writing the outputs. Here we use "Times" to denote how many times faster of parallel implementation over sequential implementation:

Times = Time of Libsvm --------------

Time of parallel implementation

By analyzing the results, we can see that MRPsvm has successfully parallelized SVM training. For all five datasets, much more time is saved when running MRPsvm than Libsvm. Especially in the first four cases, the speed of MRPsvm is more than 4 times higher than that of Libsvm. In all of the cases, MRPsvm achieves outstanding higher Times than Pisvm, indicating that MRPsvm is more suitable than Pisvm on multi-core systems. We note that in these experiments, the accuracy of the three classifiers is almost the same. In fact, the numbers of support vectors generated by these classifiers are also quite close. Actually, these three implementations have essentially the same mechanism in quadratic problem solving, i.e., to iteratively optimize one pair of variables until achieving global optimization. The difference in runtime mainly caused by the selected working set. Selecting different variables to perform optimization may induce totally different results in an iteration, but generally speaking, as long as global convergence is reached, the influence is marginal.

We show in figure 2 the times up of the two parallel solvers over the sequential solver when running on the second computer. MRPsvm again outperforms Pisvm on all of the datasets. Among them the best Times is achieved on Adult, while the worst Times is found on Kddcup99. This indicates that MRPsvm performs better on smaller datasets. The main reasons for worse performance on larger problem are due to locality and overheads of reduction. In each map, the updating of f requires accessing the whole data samples and several temporal vectors with the size close to l. Therefore, for large datasets, it is difficult to guarantee the locality for using L2 cache, especially when the cache is shared by several threads. Since we partition the global problem by columns, each map generates l intermediate f i , so that the reduction is costly when l is very large.

In fact, the parallel performance on 8 cores only slightly outperforms that on 4 cores. As explained at the beginning of this section, this is because we did not make full use of the memory on computer II. Far more time can be saved if we increase the cache size with caution to the maximum. In this optimal case, the cache size of MRPsvm and Libsvm is larger than that of Pisvm, since the former two systems generally require less memory. Therefore, it is implied that more improvements can be achieved for MRPsvm and Libsvm than Pisvm.

VI. CONCLUSION

This paper proposes a parallel implementation of SVM for multi-core and multiprocessor systems. It implements decomposition method and utilize SMO as inner solver. The parallelism is conducted to update the vector f in the decomposition step and is programmed in the simple, yet pragmatic programming framework Map-Reduce. A shared cache is designed to save the kernel matrix columns when the data size is very large. Extensive experimental results show that MRPsvm is very efficient in solving large scale problems. For instance, the speed on 4 processors can increase more than 4 times than Libsvm for most of the applications. It overwhelms the state-of-the-art Pisvm in all benchmark tests in the sense of both speed and memory requirement.

The multi-core system is already highly available in modern market, and is dominating the trend of processor development. This makes MRPsvm to be practical since it requires only easily available, affordable and single computer. These performance improvements brought by MRPsvm can extend the potential feasibility of SVM in solving increasingly challenging problems. Furthermore, the success of MRPsvm indicates that Map-Reduce is a possible option to parallelize machine learning algorithms.

However, several aspects are still worth considering for the purpose of further improvements. For instance, how to define the best granularity of map work for a particular dataset, how to improve data locality, shrinking, further parallelizing other slow operations such as the computation of inner kernel matrix.

Figure 1 .

 1 Figure 1. Problem partitioning and tasks assignment to mappers for parallel computing.

Figure 2 .

 2 Figure 2. Times up of Pisvm and MRPsvm over Libsvm when running on computer II.

Table I THE

 I PHYSICAL FEATURES OF THE MULTI-CORE SYSTEMS.

	Features	Computer I	Computer II
	# of CPUs	1	2
	# of cores	4	8
	Frequency	1600MHz*4	2327MHz*8
	Memory	2G	4G
	L2 cache	4M	6M*2

Table II DESCRIPTION

 II OF THE FIVE DATASETS AND THE TWO PARAMETERS OF SVM ON EACH DATASET.

	Dataset	# training samples	# testing samples	# Classes	# Dimensions	C	γ
	Web	24,692	25,075	2	300	64	7.8152
	Adult	32,561	16,281	2	123	100	0.5
	Mnist	60,000	10,000	2	576	10	1.667
	Covtype	435,759	145,253	8	54	10	2e-5
	Kddcup99	898,430	311,029	2	122	2	0.6

Table III THE

 III TRAINING TIME AND ACCURACY OF THE THREE SYSTEMS ON FIVE DATASETS PERFORMED ON COMPUTER I. THE UNIT OF TIME IS SECOND.

	Dataset	Time	Libsvm Accuracy	Time	Pisvm Accuracy	Times	Time	MRPsvm Accuracy	Times
	Web	306.4		97.6%	117.5	97.6%	2.6	65.8	97.6%	4.7
	Adult	311.6		82.7%	91.4	82.7%	3.4	59.2	82.7%	5.3
	Mnist	517.8		99.8%	148.7	99.8%	3.5	123.2	99.8%	4.2
	Covtype	20260.7	51.0%	5612.6	51.0%	3.6	3895.1	51.0%	5.2
	Kddcup99	726.8		92.0%	415.5	92.0%	1.7	351.9	92.0%	2.1


Hai-Xiang Zhao
email: haixiang.zhao@ecp.fr

Frédéric Magoulès
email: frederic.magoules@hotmail.com

New Parallel Support Vector Regression for Predicting Building Energy Consumption

Keywords: Support Vector Regression (SVR), Building energy consumption, Parallel computing, Multi-core, Map-Reduce

One challenge of predicting building energy consumption is to accelerate model training when the dataset is very large. This paper proposes an efficient parallel implementation of support vector regression based on decomposition method for solving such problems. The parallelization is performed on the most time-consuming work of training, i.e., to update the gradient vector f . The inner problems are dealt by sequential minimal optimization solver. The underlying parallelism is conducted by the shared memory version of Map-Reduce paradigm, making the system particularly suitable to be applied to multi-core and multiprocessor systems. Experimental results show that our implementation offers a high speed increase compared to Libsvm, and it is superior to the state-of-the-art MPI implementation Pisvm in both speed and storage requirement.

I. INTRODUCTION

Since proposed in 1990s by V.N. Vapnik, support vector machine (SVM) has been a popular supervised learning method of solving classification and regression problems [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. The support vector regression (SVR) has shown robust generalization ability in the application of predicting building energy consumption [START_REF] Dong | Applying support vector machines to predict building energy consumption in tropical region[END_REF]- [START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF]. The essential computation of SVR is to solve a quadratic problem (QP) which is both time and memory costly, causing it a challenge to solve large scale problems. Despite of several optimizing or heuristic methods such as shrinking, chunking [START_REF] Osuna | Training support vector machines: an application to face detection[END_REF], kernel caching [START_REF] Joachims | Making large-scale support vector machine learning practical[END_REF], approximation of kernel matrix [START_REF] Chang | Psvm: Parallelizing support vector machines on distributed computers[END_REF], sequential minimal optimization (SMO) [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF], primal estimated sub-gradient solver [START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF], a more sophisticated and satisfactory resolution is always expected for this challenging problem. As stated in [START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF], the building's energy system is extremely complex involving large number of influence factors, making it a common situation that we tackle large scale datasets.

With the development of chip technologies, computers with multi-core or multiprocessor are becoming more available and affordable in modern market. This paper therefore attempts to investigate and demonstrate how SVR can benefit from this modern platform when solving the problem of predicting building energy consumption. A new parallel SVR that is particularly suitable to this platform is proposed. Decomposition method and inner SMO solver compose the main procedure of training. A shared cache is designed to store the kernel columns. For the purpose of achieving easy implementation without sacrificing performance, the new parallel programming framework Map-Reduce is chosen to perform the underlying parallelism. The proposed system is therefore named as MRPsvm (abbreviation of "Map-Reduce parallel SVM"). Comparative experiments are conducted on three simulated energy consumption datasets, showing significant performance improvement in our system compared to Libsvm [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF] and Pisvm [START_REF] Brugger | Parallel support vector machines[END_REF].

The following sections are organized as follows. Section II introduces the principle of SVR training and the decomposition method. Section III explains how MRPsvm is implemented. Section IV states the related work. Section V presents how to prepare the energy consumption datasets and the numerical experiments. Conclusions are drawn in section VI.

II. SUPPORT VECTOR REGRESSION

A. Principle of ǫ-SVR

Let's present the training data as (x 1 , z 1 ), ..., (x l , z l ),where vector x i is the ith sample, z i is the ith target value corresponding to x i , l is the number of samples. The goal of ǫ-SVR is to find a decision function g(x) that makes the deviation between the calculated targets and the actual targets is at most ǫ [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. The dual form of the SVR can be written in the following quadratic form:

min α 1 2 α T Qα - 2l i=1 p i α i (1)
subject to

y T α =0 (2) 0 ≤ α i ≤ C ∀i = 1, ..., 2l (3)
where α is a vector of 2l variables, Q is a 2l by 2l positive semidefinite matrix. Each element of Q has the following form:

Q ij = K(x m , x n ) m = i if i ≤ l i -l if i > l n = j if j ≤ l j -l if j > l i, j = 1, ..., 2l
where K(x i , x j ) is called kernel function which can be substituted by any Mercer kernel, the most common ones are radial basis function (RBF), polynomial function, sigmoid function and linear function. The parameter p in (1) is defined as:

p i = ǫ + z i if i = 1, ..., l ǫ -z i if i = l + 1, ..., 2l (4)
y in the constraint ( 2) is defined as:

y i = 1 if i = 1, ..., l -1 if i = l + 1, ..., 2l (5)
In the constraint (3), C is the upper bound used to trade off between model performance on training data and its generalization ability. The objective of the problem is to find the solution of α which makes (1) minimized and the constraints (2) (3) fulfilled. After the optimum α is found, the decision function can be formulated as:

g(x) = l i=1 (-α i + α i+l )K(x i , x) + b
where b is a constant value which can be easily calculated in the training step. The training samples that satisfy (-α i + α i+l = 0) are called support vectors. It is obvious that only support vectors have contribution to the decision function, this is the reason why the algorithm is called support vector regression.

B. Decomposition Approach

In building energy application, large amounts of datasets are highly available for model analysis. Since the size of the kernel matrix Q is 2l * 2l, it is difficult to store the whole matrix in memory when l is very large. For the example of 20 buildings which will be described in section V, to store the whole Q for this dataset we need more than 74GB memory which is not affordable by ordinary users. And for the case of 50 buildings, the memory requirement is even larger. Osuna et al. [START_REF] Osuna | Training support vector machines: an application to face detection[END_REF] proposed a method to decompose the problem into smaller tasks. In each task, a working set which contains certain parts of α is chosen to be optimized, while the rest of α remains in constant value. The program repeats the selectoptimize process iteratively until global optimality conditions are satisfied. In each iteration, only the involved partition of kernel matrix needs to stay in the memory. Let B denote the working set which has n variables and N denote the nonworking set which has (2ln) variables. Then, α, y, Q and p are correspondingly written as:

α = α B α N , y = y B y N , Q = Q BB Q BN Q N B Q N N , p = p B p N
Accordingly, the small task in this case can be written as:

min 1 2 α T B Q BB α B -α T B (p B -Q BN α N ) + 1 2 α T N Q N N α N -α T N p N (6) subject to α T B y B + α T N y N = 0 (7) 0 ≤ α B ≤ C (8)
Since the last term 6) remains constant in each iteration, it can be omitted while calculating, so that function [START_REF] Osuna | Training support vector machines: an application to face detection[END_REF] has the same form as function [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. One of the advantages of this decomposition method is that the newly generated task is small enough to be solved by most off-theshelf methods. In our work, we choose SMO as the inner small task solver due to its relative simplicity, yet high performance characteristics. In fact, SMO is itself an extreme case of the decomposition method where the working set contains only two variables. This kind of binary sub-problem can be easily solved analytically [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF] [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. As stated in [START_REF] Joachims | Making large-scale support vector machine learning practical[END_REF], the solution of ( 6) is strictly feasible towards the optimum solution of global problem [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. This feature guarantees the global convergence of the decomposition method.

( 1 2 α T N Q N N α N -α T N p N ) of (
The Karush-Kuhn-Tucker (KKT) optimality condition is verified through evaluating the gradient of (1), i.e., f i = 2l j=1 α j Q ij + p i for all i = 1, ..., 2l. This procedure can be summarized as follows. First, we classify the training samples into two categories

I up (α) = {i|α i < C, y i = 1 or α i > 0, y i = -1} I low (α) = {i|α i < C, y i = -1 or α i > 0, y i = 1}
Then we search two extreme values m(α) and M (α):

m(α) = max i∈Iup(α) -y i f i (9)
M (α) = min i∈I low (α) -y i f i (10)
And then, we define the stop condition as:

m(α) -M (α) ≤ ǫ (11)
The selection of working set directly influences the speed of convergence. For inner SMO solver, maximal violating pair is selected to be the binary working set according to the second order information [START_REF] Fan | Working set selection using second order information for training support vector machines[END_REF]. We do not state here how inner SMO solver works since it has been discussed in detail in [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF] and [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. For the selection of working set B, we simply consider the first order information and select, in some sense, the maximal violating pairs as proposed by [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF]. Suppose the required size of B is n, we choose q (q < n) variables from α by sequentially selecting pair of variables which satisfy [START_REF] Platt | Fast training of support vector machines using sequential minimal optimization[END_REF] and [START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF]. The remaining (nq) variables are chosen as those who entered B in the last iteration but not yet selected in current B. The selection of these (nq) variables follows the sequence: free variables (0 < α i < C), lower bound variables (α i = 0), upper bound variables (α i = C). The reason for putting restraint on the number of new variables entering the working set is to avoid frequent entering-leaving of certain variables. Otherwise, the speed of convergence would considerably slow down [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF].

After the working set is optimized, f is updated by the newly optimized α j , ∀j ∈ B. This procedure is crucial as it prepares f for the next iteration to do optimality condition evaluation and working set selection. In fact, this is the most computational expensive step in SVR training due to the heavy work of computing Q ij . The updating procedure can be written as follows:

f * i = f i + j∈B ∆α j Q ij i = 1, ..., 2l (12)
where ∆α j is the newly optimized α j minus the old α j . The whole decomposition method is summarized in algorithm 1.

Algorithm 1 Decomposition solver of SVR Input: data set (x i , z i ), ∀i ∈ 1, ..., l Initialize: α i = 0, p i by ( 4), y i by ( 5),

f i = p i , ∀i ∈ 1, ..., 2l Calculate: I up , I low , m(α), M (α) Repeat select working set B until | B |= n update α i by SMO solver, ∀i ∈ B update f i , ∀i ∈ 1, ..., 2l calculate I up , I low , m(α), M (α) Until m(α) -M (α) ≤ ǫ III. SYSTEM IMPLEMENTATION
Next, we will introduce some key points of our system implementation, including why and how Map-Reduce is to do the parallelism, shared caching technique and the difference of our system compared to Pisvm.

A. Parallelizing the QP Solver

Map-Reduce is a new parallel programming framework originally proposed in [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF]. It allows users to write code in a functional style: map computations on separated data, generate intermediate key-value pairs and then reduce the summation of intermediate values assigned to the same key. A runtime system is designed to automatically handle low-level mapping, scheduling, parallel processing and fault tolerance. It is a simple, yet very useful framework. It can help people extract parallelism of computations on large datasets by taking advantage of distributed systems.

Problem [START_REF] Brugger | Parallel support vector machines[END_REF] can be regarded as a summation of several computational expensive terms as shown in the top right corner in Figure 1. Therefore, Map-Reduce is naturally suitable to deal with this problem. The working set B is uniformly decomposed into several small pieces, the calculation of f * is also divided into several parts in the same manner as for B. Each part is then assigned to a mapper. After the parallel calculations of these mappers, final f * is summed up by the reducer. Here, (j k , k = 1, ..., n) is the variable index of working set in kernel matrix, which gives the kth variable in B with its index in Q as j k . In practice, since some of ∆α i are so marginal that they can be omitted, it is not necessary to update f on all of the n variables.

In some recent work, Map-Reduce is proved to be an effective parallel computing framework on multi-core systems. Chu et al. [START_REF] Chu | Map-reduce for machine learning on multicore[END_REF] have developed Map-Reduce as a general programming framework on multi-core systems for machine learning applications. Phoenix designed in [START_REF] Ranger | Evaluating mapreduce for multi-core and multiprocessor systems[END_REF] implements a common API for Map-Reduce. It allows users to easily parallelize their applications without conducting concurrency management. The Map-Reduce tasks are performed in threads on multi-core systems. An efficient integrated runtime system is supposed to handle the parallelization, resource management and fault recovery by itself. This system is adopted as the underlying Map-Reduce handler in our MRPsvm. We show the parallel architecture of one iteration in Figure 1. The small tasks are uniformly distributed to mappers which have the same number of processors. Phoenix serves as the role of creating and managing mappers and reducers, making the system is easy to be implemented.

B. Kernel Matrix Caching

Since the calculation of kernel elements dominates the training work, it is an effective consideration to cache the kernel elements in memory as much as possible. MRPsvm maintains a fix-sized cache which stores recently accessed or generated kernel columns. The cache replacement policy is a simple least-recent-use strategy, the same as that of Libsvm. Only the column currently needed but not hit in the cache will be calculated. All parallel mappers share an unique copy of cache in the shared memory. As a result, the operation of inserting a new column into the cache should be synchronized.

For inner SMO solver, since the kernel matrix size is small enough and is dependent on the size of working set B which is normally set to 1024 according to the knowledge of experience, it is practical to cache the full version of this small matrix.

To reduce the storage requirements, the sample vectors x i are stored by sparse representation. When calculating a kernel column (Q ij , i = 1, ..., 2l), we need to unroll the jth sample vector to dense format and then calculate the dot products of this vector with all 2l sample vectors. There is only one copy of the whole dataset (x, z) the shared memory.

C. Comparison of MRPsvm with Pisvm

Pisvm also uses decomposition method to train SVM in parallel. It is an efficient tool to analyze multiple buildings' energy behaviors as stated in our previous work [START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF]. However its implementation is different from our new implementation MRPsvm in many aspects. First, Pisvm is based on MPI implementation and aims at extracting parallelism from distributed memory systems, while our parallel algorithm is conducted

Select B Update aa Split task mapper mapper mapper 1 1 2 2 * 1 1 1 1 1 1 1 k k m m n n j j j j j j j j j j f f Q Q Q Q Q D D D D D ' ' ' ' ' "" "" 1 1 2 2 * 2 2 2 2 2 2 2 k k m m n n j j j j j j j j j j f f Q Q Q Q Q D D D D D ' ' ' ' ' ""
""

1 1 2 2 * k k m m n n d d j dj j dj j dj j dj j dj f f Q Q Q Q Q D D D D D ' ' ' ' ' ""
"" by Map-Reduce threads on shared memory system. The two implementations are based on totally different models. Second, in Pisvm, each process stores a copy of data samples, while in the contrary, MRPsvm stores only one copy in the shared memory. This means MRPsvm can save large amount of storage when the dataset is huge. The saved space can be used to cache more kernel matrix in order to further improve training speed. Third, Pisvm adopts a distributed cache strategy in order to share the saved kernel elements among all of the processes. Each process stores locally a piece of the cache. Consequently, the work for updating gradients is divided and assigned globally to proper processors according to the cache locality. In contrast, MRPsvm has only one copy of the cache, and each processor accesses the cache equally, so that the overhead of global assignment is avoided. However, we have to note that synchronization on cache write is required.

Phoenix

In section V, we will compare the performance of Pisvm with that of MRPsvm in the application of building energy prediction, providing direct evidence that our system is more efficient and suitable than the MPI implementation on multicore systems.

IV. RELATED WORK

Several approaches have been proposed to parallelize SVM, mostly for solving classification problem. They can be classified into several categories according to the type of QP solver. Based on stochastic gradient descent method, P-packSVM optimizes SVM training directly on the primal form of SVM for arbitrary kernels [START_REF] Zhu | P-packsvm: Parallel primal gradient descent kernel svm[END_REF]. Very high efficiency and competitive accuracy have been achieved. Psvm proposed in [START_REF] Chang | Psvm: Parallelizing support vector machines on distributed computers[END_REF] is based on interior point QP solver. It approximates the kernel matrix by incomplete Cholesky Factorization. Memory requirement is reduced and scalable performance has been achieved. Bickson et al. [START_REF] Bickson | A gaussian belief propagation solver for large scale support vector machines[END_REF] solve the problem by Gaussian belief propagation which is a method from complex system domain. The parallel solver brings competitive speedup on large scale problems. The decomposition method attracts more attention than the above solvers. Graf et al. [START_REF] Graf | Parallel support vector machines: The cascade svm[END_REF] train several SVMs on small data partitions, then they aggregate support vectors from two pair SVMs to form new training samples on which another training is performed. The aggregation is repeated until only one SVM remains. The similar idea is adopted by Dong et al. [START_REF] Dong | Fast svm training algorithm with decomposition on very large data sets[END_REF], in their work, sub-SVMs are performed on block diagonal matrices which are regarded as the approximation to the original kernel matrix. Consequently, nonsupport vectors are removed when dealing these sub-problems. Zanni et al. [START_REF] Zanni | Parallel software for training large scale support vector machines on multiprocessor systems[END_REF] parallelize SVM-light with improved working set selection and inner QP solver. Hazan et al. [START_REF] Hazan | A parallel decomposition solver for svm: Distributed dual ascend using fenchel duality[END_REF] propose a parallel decomposition solver using Fenchel Duality. Lu et al. [START_REF] Lu | Parallel randomized sampling for support vector machine (svm) and support vector regression (svr)[END_REF] parallelize randomized sampling algorithms for SVM and SVR.

Cao et al. [START_REF] Cao | Parallel sequential minimal optimization for the training of support vector machines[END_REF] and Catanzaro et al. [START_REF] Catanzaro | Fast support vector machine training and classification on graphics processors[END_REF] parallelize SMO solver for training SVM for classification. Both work mainly focuses on updating gradient for KKT condition evaluation and the working set selection. The difference between them lies on the implementation details and the programming models. Specifically speaking, the first work is conducted by using MPI on clusters while the second one by Map-Reduce threads on modern GPU platform. In our work, we also adopt SMO algorithm. But we use it as the inner QP solver without any parallel computation, in fact, we perform the parallelization on external decomposition procedure. The main advantage of our coarse-grained parallelism is that it can significantly reduce the burden of overheads since the number of iterations in global decomposition procedure (where n ≫ 2) is extremely smaller than that of pure SMO algorithm (where n = 2). Although GPU SVM [START_REF] Catanzaro | Fast support vector machine training and classification on graphics processors[END_REF] and our system are implemented in threads, we will not compare them in experiments since they are designed for different platforms and GPU SVM is specially used to solve classification problems.

V. PREDICTING BUILDING ENERGY CONSUMPTION

In this part, we present the application of SVR on predicting building energy consumption. The proposed parallel SVR is compared with Libsvm and Pisvm on three datasets to demonstrate how much performance improvement can be achieved by the new implementation. Libsvm is a widely used sequential implementation while Pisvm is the state-of-the-art parallel implementation. Although Pisvm is not especially designed for multi-core architecture, we still have good reasons for doing this comparison. Firstly, as the best knowledge as we know, there is no existing parallel implementation of SVR that is specially developed for multi-core systems. Therefore there is a strong need to verify if our system could outperform the parallel implementation designed for distributed memory. Secondly, most of the systems surveyed in section IV are not available to the public, while Pisvm, as a typical parallel implementation of SVM, is easy to obtain. Thirdly, the QP solver of Pisvm is the same as MRPsvm, hence, if we compare our system with Pisvm, the advantage of Map-Reduce framework is more convincing.

A. Energy Consumption Datasets

Three datasets are prepared for the model training. They denote the historical energy consumption of one building, 20 buildings and 50 buildings respectively. The datasets are simulated in Energyplus which is a widely applied, state-ofthe-art and effective building energy simulation tool [START_REF] Energyplus | [END_REF].

The buildings are designed as follows. All of them are for office use and located in urban area. We evenly distribute them into five typical cities of France which are Paris-Orly, Marseilles, Strasbourg, Bordeaux and Lyon. As outlined in Figure 2, the five cities vary remarkably in ambient dry bulb temperatures, making the datasets represent energy requirements under five typical weather conditions. Each building has similar structures, i.e., single-story, mass-built, one rectangle room with attic roof and four windows without shading. Electrical equipments including lighting system, fans, water heaters, are scheduled as common office use. In winter season (from November 1st to March 31st), district heating is applied in order to keep the room temperature at a constant level. Ventilation is adopted for indoor thermal comfort. The number of occupants depends on the housing space and people density, with the average of 0.2 people per zone floor area. During the simulation, some input variables such as size, orientation, window area, scheduling, are set differently to achieve diversity among multiple buildings. The dataset of one building is hourly energy dynamics in a period of one year. We record 8 important features [START_REF] Zhao | Feature selection for predicting building energy consumption based on statistical learning method[END_REF] as shown in Table I. For other two datasets, the recording period is from November to March which is winter season in France, and we record 4 more features which generate the building diversity, i.e., height, length, width and window/wall area ratio. Since people usually do not work in weekends or holidays, the energy requirement in these days is quite small compared to normal working days. This means weekends and holidays have totally different energy behaviors from working days. Therefore, to simplify the model in our practice, we only use the consumption data of working days in the experiments. One more building is simulated for model evaluation purpose. The attributes of the three datasets are shown in Table II.

B. Experiments and Results

We perform the experiments on two shared memory computers which represent two different hardware architectures. As outlined in Table III, the first computer has 2 cores while the second one has 4. Both of them have a shared L2 cache and memory and running 64 bit Linux system (kernel version 2.6.27-7). We train all SVRs with Gaussian kernel, i.e., K(x i , x j ) = exp{-γ x ix j 2 }. The parameters C and γ are evaluated by 5-fold cross validation and ǫ is set to 0.01, as shown in the last three columns of Table II. Since the caching technique is crucial for performance, for a reliable comparison, we set the cache size to be the same for all three SVR implementations. Furthermore, we restrict the cache size to be far smaller than the memory size in order to minimize page faults in runtime. Here we have to emphasize that the following reported performance might not be the best for all three implementations, only serving for comparison purpose. Two performance evaluation methods are applied in this work. One is the mean squared error (MSE) which gives the average deviation of the predicted values to the real ones. Another is the squared correlation coefficient (SCC) which lies in [0, 1] and gives the ratio of successfully predicted number of target values on total number of target values.

On each dataset, we train the sequential implementation Libsvm, and the parallel implementations Pisvm and MRPsvm on both computer-I and computer-II. Table IV shows the results of the three implementations performed on dual-core processor, including the number of support vectors (nSVs), MSE, SCC and training time. We show the training time on quad-core system in Table V. Since nSVs, MSE and SCC in quad-core case are the same as those in dual-core case, we omit them in Table V. Note that the time columns represent the whole training time, i.e., from reading the problem to writing the outputs.

The results demonstrate that MRPsvm has successfully parallelized SVR training. For all three datasets, the accuracy and the nSVs of MRPsvm are quite close to that of Libsvm and Pisvm. Much more time is saved when running MRPsvm than Libsvm and Pisvm in all tests. Here we use speedup to demonstrate how many times faster is the parallel implementation than sequential implementation: ------------------Training time of parallel implementation We show the speedup of MRPsvm and that of Pisvm comparatively in Figures 3,[START_REF] Lai | Vapnik's learning theory applied to energy consumption forecasts in residential buildings[END_REF][START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF] For the case of one building, the speed of MRPsvm is twice higher than that of Libsvm. For the case of 20 and 50 buildings, MRPsvm performs more than 16 times faster than Libsvm. On both computers and on all three datasets, MRPsvm achieves remarkable higher speedup than Pisvm, indicating that MRPsvm is more suitable than Pisvm on multi-core systems. The speed improvement by MRPsvm is particularly obvious in multiple buildings cases, indicating that MRPsvm performs better than Pisvm on large scale datasets. However, the better performance is not guaranteed on even larger problems due to locality and overheads of reduction. In each mapper, the updating of f requires the access of the whole data samples and several temporal vectors which have the size close to 2l. Therefore, for large datasets, it is difficult to guarantee the locality for using L2 cache which is shared among several threads. Since we partition the global problem by columns, each mapper generates 2l intermediate f i , the reduction is costly when l goes to a large scale.

Speedup = Training time of Libsvm

It is not surprising that the parallel performance on quadcore only slightly outperforms that on dual-core. One reason is that the two computers have different architectures, e.g., the single processor frequency of computer-II is half slower than that of computer-I. The other reason is that we set the cache size to be the same for these two architectures and did not make full use of memory. This means far more time could be saved if we increase the cache size with caution to the maximum. In this optimum case, the cache size of MRPsvm and Libsvm is larger than that of Pisvm, since the former two systems generally require less storage space. It is implied that more improvements can be achieved for MRPsvm than Pisvm. We note that the accuracy of MRPsvm is almost the same as that of Libsvm. In fact, these three implementations have essentially the same mechanism in QP solving, i.e., to iteratively optimize one pair of variables until achieving global optimization. Their difference comes from the selection of variable pair which may induce totally different results in a sub-task. But in the global point of view, as long as convergence is reached, the influence of internal difference is marginal.

VI. CONCLUSION

We focus in parallel SVR for solving large scale problems since it is suitable to the building energy analysis. The proposed MRPsvm implements decomposition method with SMO as the sub-problem solver. The simple, but pragmatic programming framework Map-Reduce is adopted to conduct parallel updating of the gradient vector f . The system is relatively easy to be implemented and very high performance is achieved.

Experimental results show that the new system provides the same accuracy as Libsvm does, yet performs far more efficiently than the sequential implementation in solving stated problems. On the smallest dataset, MRPsvm achieves a morethan-twice speed times up than Libsvm while on the largest dataset, the speed times up can reach higher than 16-fold and 19-fold on dual-core system and on quad-core system, respectively. The proposed implementation is superior to the state-of-the-art Pisvm in all tests in the sense of both speed and memory requirement. Since the multi-core system is dominating the trend of processor development and is highly available in modern market, MRPsvm is potentially very practical and feasible in solving large scale regression problems. Furthermore, the success of MRPsvm indicates that Map-Reduce is a possible option to parallelize machine learning algorithms.

However, the proposed system is not yet mature, there are still several aspects worth considering for further improvements, for instance, shrinking, finding the best granularity of parallel work for a particular dataset.

Fig. 1 :

 1 Fig. 1: Architecture of the parallelization in one (d = 2l).

Fig. 2 :

 2 Fig. 2: Dry bulb temperature in the first 11 days of January.

Fig. 4 :

 4 Fig. 4: Speedup of Pisvm and MRPsvm on 20 buildings' data.

Fig. 5 :

 5 Fig. 5: Speedup of Pisvm and MRPsvm on 50 buildings' data.

TABLE I :

 I The features of one building's consumption data.

	Features	Unit
	Outdoor Dry Bulb	C
	Outdoor Air Density	kg/m 3
	Water Mains Temperature	C
	Number of Occupants	-
	Lights Total Heat Gain	J
	Electric Equipment Total Heat Gain	J
	Zone Mean Air Temperature	C
	Zone Infiltration Volume	m 3

TABLE III :

 III The physical features of the experimental environment.

	Features	Computer-I	Computer-II
	# of CPUs	1	1
	# of cores	2	4
	Frequency	3.4GHz*2	1.6GHz*4
	Memory	2G	2G
	L2 cache	2M	4M

TABLE II :

 II Description of the three datasets and the three parameters of SVR on each dataset.

	Dataset	# training samples	# testing samples	# Dimensions	C	γ	ǫ
	One building	5064	1008	8	32	0.594	0.01
	20 buildings	49940	2498	12	16	0.4193	0.01
	50 buildings	124850	2498	12	16	0.4357	0.01

TABLE V :

 V The training time of the three predictors performed on computer-II. Time unit is second.

	Dataset	Libsvm	Pisvm	MRPsvm
	1 building	18.0	7.3	6.9
	20 buildings	2532.1	214.1	133.5
	50 buildings	32952.2	2325.5	1699.0

Fig. 3: Speedup of Pisvm and MRPsvm on one building's data.

TABLE IV :

 IV The training time and performance of the three predictors on three datasets performed on computer-I. bd: building, nSVs: number of support vectors, MSE: mean squared error, SCC: squared correlation coefficient. The unit of time is second.

	Data	nSVs	Libsvm MSE SCC	Time	nSVs	Pisvm MSE	SCC	Time	nSVs	MRPsvm MSE SCC	Time
	1 bd	2150	6.16e-4	0.97	22.3	2162	6.10e-4	0.97	9.2	2168	6.14e-4	0.97	9.0
	20 bd	9014	2.11e-3	0.96	3407.0	8967	2.12e-3	0.96	339.5	8970	2.08e-3	0.96	212.7
	50 bd	22826	3.73e-4	0.97	44212.5	22823	3.74e-4	0.97	4179.8	22799	3.73e-4	0.97	2745.8


Frédéric Hai-Xiang Zhao
email: haixiang.zhao@ecp.fr

Magoul Ès

Feature Selection for Support Vector Regression in the Application of Building Energy Prediction

Keywords: support vector regression, feature selection, building, energy consumption

When using support vector regression to predict building energy consumption, since the energy influence factors are quite abundant and complex, the features associated with the statistical model could be in large quantity. This paper focuses in feature selection for the purpose of reducing model complexity without sacrificing performance. The optimal features are selected by their feasibility of obtaining and the evaluation of two filter methods. We test the selected subset on three datasets and train support vector regression with two different kernels: radial basis function and polynomial function. Extensive experiments show that the proposed method can select valid feature subset which guarantees the model accuracy and reduces the computational time.

I. INTRODUCTION

Support vector regression (SVR) has shown remarkable generalization abilities in the prediction of energy load [START_REF] Dong | Applying support vector machines to predict building energy consumption in tropical region[END_REF][START_REF] Espinoza | Fixed-size least squares support vector machines: a large scale application in electrical load forecasting[END_REF][START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF]. In our application, the statistic model is trained from the historical consumption data and then the derived model is used to predict the energy behaviors of one building in the future or for some new buildings. The energy behavior of buildings is influenced by many factors, such as ambient weather conditions, building's structures, occupants' behaviors, inner facilities, etc, therefore, the possible features of the samples could be in large quantity. However, not all of the features that are important to the statistic model. How to reasonably choose the useful features to be used in model learning is one of the key issues for machine learning methods on such applications.

To the best of our knowledge, there is little work concerning the feature selection for machine learning methods in the application of predicting building energy consumptions. Most of the existing works derive their own model using previously established sets of features. Azadeh et al. [START_REF] Azadeh | Integration of genetic algorithm, computer simulation and design of experiments for forecasting electrical energy consumption[END_REF] and Maia et al. [START_REF] Maia | A methodology for shortterm electric load forecasting based on specialized recursive digital filters[END_REF] forecast electrical energy consumption through analyzing the varying inner targets without any contributory variable involved. Neto et al. [START_REF] Neto | Comparison between detailed model simulation and artificial neural network for forecasting building energy consumption[END_REF] built their neural network based on the input of daily average values of dry bulb temperature, relative humidity, global solar radiation and diffuse solar radiation. Tso et al. [START_REF] Tso | Predicting electricity energy consumption: a comparison of regression analysis, decision tree and neural networks[END_REF] used more than 15 features in their assessment of traditional regression analysis, decision tree and neural networks. Yokoyama et al. [START_REF] Yokoyama | Prediction of energy demands using neural network with model identification by global optimization[END_REF] considered only two features, air temperature and relative humidity in their neural network model. Other similar approaches can be found in [START_REF] Dong | Applying support vector machines to predict building energy consumption in tropical region[END_REF][START_REF] Zhao | Parallel support vector machines applied to the prediction of multiple buildings energy consumption[END_REF][START_REF] Wong | Artificial neural networks for energy analysis of office buildings with daylighting[END_REF][START_REF] Lai | Vapnik's learning theory applied to energy consumption forecasts in residential buildings[END_REF].

This paper proposes an heuristic approach for selecting subset of features, and systematically analyzes how it will influence the model performance. The models are trained by support vector regression (SVR) with two different kernel methods based on three datasets. The feature selection method is evaluated by comparing the models' performances before and after feature selection is performed. Next, section II briefly introduces the principles of SVR model and feature selection algorithms, section III will introduce the extensive experiments and evaluate the proposed feature selection method, and finally, section IV concludes this paper.

II. SVR AND FEATURE SELECTION

A. The Principles of SVR

The basic idea of SVR is to find a model function f (x) to represent the relationship between the features and the target [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. We represent the whole training data as (x 1 , y 1 ), (x 2 , y 2 ), ..., (x l , y l ), where vector x i is the ith sample of the features, y i is the corresponding target, l is the number of samples. If there are n features in the training data, then the dimensionality of x i is n. The model function has the form f

(x) = l i=1 (α * i -α i )K(x i • x) + b
, where K is kernel function which implicitly maps the problem from lower dimensional space to higher feature space, α * i and α i are variables which are optimized by maximizing the following quadratic function:

W (α * i , α i ) = - 1 2 l i,j=1 (α * i -α i )(α * j -α j )K(x i • x j ) + l i=1 y i (α * i -α i ) -ε l i=1 (α * i + α i ) under the constraints l i=1 α * i = l i=1 α i and 0 ≤ α * i , α i ≤ C, i = 1, 2, ..., l.
Then the problem can be solved by many off-the-shelf quadratic problem solvers. Usually, only certain parts of the samples can satisfy the property: α * iα i = 0, they are called support vectors (SVs) since only them have contributions to the model. The popular choices of the kernel are linear function, sigmoid function, radial basis function (RBF) and polynomial function. They represent different mappings from lower space to higher feature space. In our data training process, we choose the last two kernels for the purpose of solving non-linear problems.

B. Algorithems of Feature Selection

The aim of feature selection is to select the most useful feature set in order to establish a good predictor for the concerned learning algorithm. The irrelevant and unimportant features are discarded in order to reduce the dimensionality. Several advantages will be achieved if we wisely perform the selection. The first one is to simplify the calculation while keeping the dimensionality minimized, which could contribute to avoid the curse of dimensionality. The second one is the possible improvement of accuracy of the developed model. The third one is the improved interpretability of the models. The last one is the feasibility of obtaining accurate feature samples, especially for some time series problems in practice.

In our practice, two feature selection methods are used to rank each feature, then the features are filtered according to their ranks. The first one gives each feature a score depending on the correlation coefficient between this feature and the target. We use CC to stand for this method. It is defined as:

CC(f ) = N XY -( X)( Y ) [N X 2 -( X) 2 ][N Y 2 -( Y ) 2 ]
The other one is called Regression, Gradient guided feature selection (RGS), which is firstly proposed by A. [START_REF] Navot | Nearest neighbor based feature selection for regression and its application to neural activity[END_REF]. The basic idea is to assign a weight to each feature and then optimize the weight vector of all the features simultaneously by gradient ascent. The non-linear function K-Nearest-Neighbor is applied in the procedure of optimization. The estimated target of sample x is defined as:

Navot et al. when they use regression to analyze brain neural activities

fw (x) = 1 Z x ′ ∈N (x) f (x ′ )e -d(x,x ′ )/β where N (x) is the set of K nearest neighbors of sample x. d(x, x ′ ) = n i=1 (x ′ i -x i ) 2 w 2 i
is the distance between sample x and one of its nearest neighbors x ′ , n is the number of features, w i is the specific weight assigned to ith feature and w is the weight vector. Z = x ′ ∈N (x) e -d(x ′ ,x)/β is a normalization factor where β is a Gaussian decay factor. Then the optimal w can be found by maximize the following evaluation function:

e(w) = - 1 2 x∈S (f (x) -fw (x)) 2
where S is the samples for model training. Since e(w) is smooth almost everywhere in a continuous domain, one can solve the extremum seeking problem by gradient ascent.

III. NUMERICAL EXPERIMENTS A. Preparing Datasets of Building Energy Consumption

Hourly electricity consumption together with hourly values of several influencing factors is recorded as the building energy consumption data. This data is generated by Ener-gyPlus which is a widely applied, state-of-the-art building energy simulation tool [START_REF] Energyplus | [END_REF]. The test buildings are singlestory mass-built buildings for office use located in five major cities of France, which are Paris-Orly, Marseilles, Strasbourg, Bordeaux and Lyon. They represent the office energy requirements under five typical weather conditions in France. All of the buildings are located in urban area. Each building has similar structures, e.g., one rectangle room with attic roof and four windows in four directions without shading. Electrical equipments including lighting system, fans, water heaters, are scheduled as common office use. The number of occupants depends on the housing space, with the average of 0.2 people per zone floor area. Since the buildings are designed for office use, the energy consumption in holidays (including weekends) is quite small compared to weekdays. In order to simplify the model in practice, we only record the consumption data in weekdays.

Three datasets are generated for model training and testing. The first dataset is one year's consumption of one building which is located in Pairs-Orly. The hourly electricity demands together with hourly behaviors of 23 variables are recorded through one year. The consumption series from January 1st to October 31st is used for training and the samples in remaining two months, from November 1st to December 31st, is used for model testing. In other two datasets, we record the consumption for multiple buildings in winter season which means from November 1st to March 31st. We design 20 and 50 buildings to be the training data for the second and the third dataset respectively. For model testing purpose, we simulate one more building. The differences among these multiple buildings mainly rely on weather conditions, building structures and number of occupants. The buildings have diverse characteristics with randomly generated parameters: length, width, height and window/wall area ratio. The features of one building are listed in Table I.

B. Model Selection

To fully investigate how feature selection on these three datasets influence SVR models, two kernels are involved. The first one is RBF kernel: exp(-γ x ix j 2 ). Three parameters of the learning algorithm, C, p and γ are chosen by stepwise 5-fold cross validation as described in the technical document of Libsvm [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. The initial searching spaces are {2 -3 , 2 -2 , ..., 2 8 }, {2 -10 , 2 -9 , ..., 2 -5 } and {2 -10 , 2 -9 , ..., 2 2 } for C, p and γ respectively. In order to avoid numerical problems in the calculations, before training the model, we scale the training set into the range [0, 1] and then apply the scaling function to scale the testing set.

We also test the feasibility of polynomial kernel (γ * x ′ i * x j + coef ) d which is also applicable on non-linear problems. The kernel parameters coef is set to zero, γ and d are estimated by 5-fold cross validation in searching spaces {2 -10 , 2 -9 , ..., 2 2 } and {2, 3, ..., 7} respectively. In the evaluation step, Mean Squared Error (MSE) and Squared Correlation Coefficient (SCC) are adopted to evaluate the performance of the predictor.

The efficient, widely applied open source tool Libsvm [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF] is used to practice our SVR training.

C. Selection of Features

In practice, our feature selection is based on two criteria. The first one is that the selected features can be easily obtained in reality, i.e., they could be collected from measurements, surveys, related documents like building plans and so on. The second one is that they should potentially be the most important ones to the predictor. The combination of RGS and CC algorithms is used to evaluate the usefulness of features for SVR model. The scores for each feature are listed in Table I columns 2 and 3. As our extensive experimental results have shown, the features with the highest scores under RGS are generally more useful than those with the highest ranks according to CC. This indicates that RGS method is more applicable to SVR than CC method. However, since the feature subsets with low scores are still possibly useful for the learning algorithms [START_REF] Isabelle | An introduction to variable and feature selection[END_REF], we take both RGS and CC into consideration while choosing the features.

The weather condition is normally considered to be the most important factor which influence the building energy consumption. In our selection, we keep two weather features that have the highest scores under RGS, which are dry bulb temperature and outdoor air density while other weather features like relative humidity are discarded no matter how their variations could contribute to energy requirement as we thought naturally. The water main temperature, electrical equipment heat gain such as lights heat gain are determined by their power and occupants' schedule, they could probably be measured or assessed in actual buildings. Zone mean air temperature and zone infiltration volume could also be measured or estimated in a normally operated building. All of the above selected features have scores not less than 1. A special case we have to consider is the number of occupants. This feature takes a middle place under RGS, but since it can be easily counted in real life and has a very high score under the evaluation of CC, we choose to keep it in the final subset. All other features will be discarded since they get low scores or are hard to be collected in reality. For example, Zone Total Internal Total Heat Gain is difficult to be obtained directly and District Heating Outlet Temp is useless according to CC. The selected features are indicated with stars in column Case1 in Table I.

In order to evaluate if the selected feature set is optimal, we also design three other subsets for comparison, they are I.

Compared with case 1, in case 2, we change three of the selected features to other three unselected ones. In case 3, all of the selected features are substituted with other unselected ones except Zone Total Internal Total Heat Gain which is not easily obtained in practice. In case 4, two features which gain lowest scores are removed from selected subset.

D. Numerical Results

Based on the above considerations, we generate four subdatasets of the first dataset for model training and testing. We show the results of all four cases in Table II. Two conclusions can be reached accordingly, the first one is that the designed feature selection method is valid due to model performance in case 1 outperforms other three cases. The other one is that SVR model with RBF kernel has stable performance since high prediction accuracy is always achieved on all of the four subsets.

We put the model performance before and after feature selection on three datasets in Table III. After feature selection, the accuracy of the prediction on 50 buildings' consumptions improves significantly. With regard to 20 building's consumptions, MSE increases to a certain extent, indicating centain decrease in prediction accuracy. However from the standpoint of SCC, the performance of the model with RBF kernel is quite close to the situation without feature selection performed. With regard to polynomial kernel, when training on the original datasets, the prediction ability of the model is just as good as RBF kernel, indicating that polynomial kernel is also applicable on such problem. After adopting feature selection, the performance of the model becomes better for the case of 50 buildings than for 20 buildings. Same trend can also be found for RBF kernel. These phenomena indicate that proposed feature selection approach could give better performance to the models when relatively more training samples are involved.

Another advantage of feature selection we can find in the results is the reduction of training time. For the case of 1 building, 20 buildings and 50 buildings, when we train the model with RBF kernel, before feature selection, the training time is 19.9s, 8208.6s, 40816.6s respectively, while after feature selection, the time decreases to 7s, 2018.4s and 29523.8s respectively. This reduction is remarkable for all of the three cases. We note that for each dataset, we have to assign different model parameters for the learning algorithm,

IV. CONCLUSION

This paper proposes a feature selection method for support vector regression when applied to predict building's energy consumption. The features are selected according to their feasibility in practice and the evaluation under two selection algorithms which can give each feature a score according to its usefulness to the predictor: the gradient guided feature selection and the correlation coefficients. Experimental results show that the selected subset is valid and can provide acceptable predictors. Performance improvement is achieved in some cases, e.g., accuracy enhanced remarkably for the models with either radial basis function or polynomial kernel on fifty buildings' data, the time for model learning decreases to a certain extent. Besides radial basis function kernel, we proved that polynomial kernel is also applicable to this application. However it seems not as stable as radial basis function kernel. Furthermore, it requires more complicated pre-processing work since compared to radial basis function, there are more kernel parameters need to be estimated.

This work serves as the first guide for selecting an optimal subset of features when applying machine learning methods to the prediction of building energy consumption. Future work will be carried out on developing appropriate feature selection algorithms for other modeling methods.

Table II COMPARISON

 II OF MODEL PERFORMANCE ON DIFFERENT FEATURE SETS. (NF: NUMBER OF FEATURES, MSE: MEAN SQUARED ERROR, SCC: SQUARED CORRELATION COEFFICIENT.)

		Case1	Case2	Case3	Case4
	NF	8	8	14	6
	MSE	6.2e-4	7.5e-4	2.1e-3	9.2e-4
	SCC	0.97	0.96	0.90	0.96
	indicated by columns Case2, Case3 and Case4 in Table

Table I THE

 I SCORES OF FEATURES EVALUATED BY RGS AND CC SELECTION METHODS. THE STARS INDICATE SELECTED FEATURES IN THAT CASE.

			Features				RGS	CC	Case1	Case2 Case3	Case4
			Outdoor Dry Bulb			1.61	0.29	*	*	*
			Outdoor Relative Humidity		0.62	0.26	*	*
			Wind Speed			0.52	0.01	*	*
			Direct Solar			0.54	0.47	*
			Ground Temperature			0.99	0.07	*
			Outdoor Air Density			1.26	0.20	*	*
			Water Mains Temperature		1.30	0.07	*	*
			Zone Total Internal Total Heat Gain	1.01	0.67
			People Number Of Occupants		0.93	0.68	*	*
			People Total Heat Gain		0.93	0.68	*
			Lights Total Heat Gain		1.13	0.05	*	*	*
			Electric Equipment Total Heat Gain	1.06	0.69	*	*	*
			Window Heat Gain for each wall	1.03	0.62	*
			Window Heat Loss for each wall	0.93	0.50	*
			Window Heat Gain for each wall	0.82	0.35	*
			Window Heat Loss for each wall	0.82	0.48	*
			Window Heat Gain for each wall	0.73	0.56	*
			Window Heat Loss for each wall	0.82	0.48	*
			Window Heat Gain for each wall	0.89	0.56	*
			Window Heat Loss for each wall	0.95	0.50	*
			Zone Mean Air Temperature		1.14	0.22	*	*
			Zone Infiltration Volume		1.00	0.34	*	*
			District Heating Outlet Temp		0.95	7.35e-4	*	*
			Table III
	PREDICTION RESULTS OF SVR WITH TWO KERNEL METHODS ON THREE
	DATASETS. (BD: BUILDING, BF: BEFORE FEATURE SELECTION, AF:
		AFTER FEATURE SELECTION.)
			One bd	20 bds	50 bds
	RBF	BF	MSE SCC	4.8e-4 0.97	4.3e-4 0.97	4.4e-4 0.97
	kernel	AF	MSE SCC	6.2e-4 0.97	2.1e-3 0.96	3.7e-4 0.97
	Polynomial	BF	MSE SCC	8.0e-4 0.96	5.8e-4 0.96	5.9e-4 0.96
	kernel	AF	MSE SCC	2.1e-3 0.91	0.19 0.85	4.7e-4 0.98
	as explained in section III-B, which always have a great
	influence on the training speed.