Functional-bandwidth kernel for Support Vector Machine with Functional Data_An alternating optimization algorithm

European Journal of Operational Research 275 (2019) 195–207 Contents lists available at ScienceDirect European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor Decision Support Functional-bandwidth kernel for Support Vector Machine with Functional Data: An alternating optimization algorithm R. Blanquero a , E. Carrizosa a , A. Jiménez-Cordero a , ∗, B. Martín-Barragán b a Departamento de Estadística e Investigación Operativa, and Instituto de Matemáticas de la Universidad de Sevilla (IMUS), Facultad de Matemáticas. Universidad de Sevilla. C/ Tarfia s/n, 41012 Sevilla, Spain b Business School, 29 Buccleuch Place, University of Edinburgh, EH89JS, Edinburgh, UK a r t i c l e i n f o Article history: Received 8 May 2018 Accepted 8 November 2018 Available online 24 November 2018 Keywords: Data mining Functional Data classification Parameter tuning SVM Functional bandwidth a b s t r a c t Functional Data Analysis (FDA) is devoted to the study of data which are functions. Support Vector Machine (SVM) is a benchmark tool for classification, in particular, of functional data. SVM is frequently used with a kernel (e.g.: Gaussian) which involves a scalar bandwidth parameter. In this paper, we propose to use kernels with functional bandwidths. In this way, accuracy may be improved, and the time intervals critical for classification are identified. Tuning the functional parameters of the new kernel is a challenging task expressed as a continuous optimization problem, solved by means of a heuristic. Our experiments with benchmark data sets show the advantages of using functional parameters and the effectiveness of our approach. ©2018 Elsevier B.V. All rights reserved. 1. Introduction Functional Data Analysis (FDA) has received considerable attention from researchers, ( Ferraty & Vieu, 2006; Ramsay & Silverman, 20 02, 20 05; Wang, Chiou, & Müller, 2016 ) and practitioners in many different fields, such as spectrometry, meteorology, ( MartínBarragán, Lillo, & Romo, 2014 ), client segmentation, ( Laukaitis & Ra ˇ ckauskas, 2005 ), speech recognition, ( Rossi & Villa, 2006 ), or physical, ( Muñoz & González, 2010; Tuddenham & Snyder, 1954 ), and chemical processes, ( Blanquero et al., 2016a; Blanquero, Carrizosa, Jiménez-Cordero, & Rodríguez, 2016b ). FDA can be considered as a generalization of the standard multivariate analysis to address problems in which data have an infinite-dimensional nature. The direct application of classic methods of multivariate analysis on infinite-dimensional data may have dramatic consequences in the obtained results. The curse of dimensionality is a clear example of this situation. Indeed, although theoretically data are described as functions, in practice functional data are represented by high dimensional vectors, yielding problems in which the number of observations is lower than the number of features and which cannot be handled by standard multivariate analysis tools. Furthermore, it is worthwhile to mention ∗Corresponding author. E-mail addresses: rblanq[email protected] (R. Blanquero), [email protected] (E. Carrizosa), [email protected] (A. Jiménez-Cordero), [email protected] (B. MartínBarragán). that the methodologies used for multivariate vectors do not exploit the functional behavior of the data since the high correlations among the different coordinates are not taken into account. In this work, we focus on a challenging problem in FDA: functional binary classification, i.e., how to classify functional data into two predefined classes using the information provided by a training sample ( Baíllo, Cuevas, & Cuesta-Albertos, 2011; Biau, Bunea, & Wegkamp, 2005; Cuevas, Febrero, & Fraiman, 2007; Ferraty & Vieu, 2006; Preda, Saporta, & Lévéder, 2007 ). Support Vector Machine (SVM) ( Carrizosa & Romero Morales, 2013; Cauwenberghs & Poggio, 2001; Cortes & Vapnik, 1995; Cristianini & Shawe-Taylor, 20 0 0; Lessmann & Voß, 20 09; Maldonado, Pérez, & Bravo, 2017; Maldonado, Weber, & Basak, 2011; Suykens & Vandewalle, 1999; Wang, Zheng, Yoon, & Ko, 2018 ) is one of the most used tools in multivariate classification, and it has also been widely applied for functional data. See Blanquero, Carrizosa, Jiménez-Cordero, and Martín-Barragán (2017) , Jiménez-Cordero and Maldonado (2018) , Martín-Barragán et al. (2014) , Muñoz and González (2010) , Rossi and Villa (2006) , and Rossi and Villa (2008) among others. As stated before, solving functional data problems, and more specifically, the functional classification problem, implies the use of specific techniques that take advantage of the intrinsic functional nature of the data. For SVM, ( Rossi & Villa, 2006 ) exploits the functional behavior of the data by adapting the classical kernels to functional kernels through the so-called transformation-based and projection-based kernels. Nevertheless, the whole range of the data is weighted with a single scalar bandwidth. https://doi.org/10.1016/j.ejor.2018.11.024 0377-2217/© 2018 Elsevier B.V. All rights reserved. 196 R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 The functional nature of the data is taken into account in Kästner, Hammer, Biehl, and Villmann (2012) , generalizing the work done in the multivariate case in Hammer and Villmann (2002) ; Sato and Yamada (1996) . Data are classified according to a dissimilarity measure with a functional weight. Such a functional weight is represented in terms of simple basis functions whose parameters are sought via stochastic gradient. To the best of our knowledge, no strategy has been presented in the literature using a supervised tool, e.g., SVM, in which different ranges in the domain of the functions are optimally selected by means of a functional weight in the kernel. Therefore, one of the main contributions of this paper is to define a new functional kernel. Such kernel has a functional bandwidth that optimally weighs the different values of the domain of the function. Similar ideas have been used in references such as Bugeau and Pérez (2007) , Chen, Wynne, Goulding, and Sandoz (20 0 0) , Duong, Cowling, Koch, and Wand (2008) , and Sain (2002) for kernel density estimation, and in Cai, Fan, and Yao (20 0 0) and Wu, Chiang, and Hoover (1998) for functional regression. We propose to embed the new functional kernel into an SVM algorithm. Both the kernel and the SVM parameters are tuned with a surrogate of the accuracy, namely, the correlation between the actual class and the SVM score. See also Berrendero, Cuevas, and Torrecilla (2016) , Székely, Rizzo, Bakirov et al. (2007) , and Torrecilla Noguerales (2015) for more details on the use of surrogate measures for the accuracy. Such parameter tuning yields a continuous optimization problem, allowing us to use gradient methods, known to be more efficient than the optimization methods available for piecewise constant performance measures, such as the misclassification rate. Moreover, the proposed method is enhanced by defining a hierarchy of kernel bandwidths models of increasing complexity, inspired by the nested model previously proposed for Multiple Kernel Learning in Carrizosa, Martín-Barragán, and Romero Morales (2014) . Using this hierarchy provides wide flexibility since complex parameterizations of the functional bandwidth can be efficiently optimized from more simple ones. The remainder of the paper is structured as follows. In Section Section 2 we present the SVM classification model for functional data. Section 3 describes the optimization method used to tune the bandwidth parameters. We focus on the alternating procedure proposed to this purpose, and on the structure of the hierarchy of kernels. Section 4 is devoted to the numerical experiments, showing that our approach outperforms the method in which one single scalar parameter bandwidth is chosen. Finally, some conclusions and extensions are described in Section 5 . 2. Functional bandwidth In this section, we formulate the SVM problem for functional data classification. See Cristianini and Shawe-Taylor (20 0 0) for a broader and more comprehensive presentation of SVM. We have a sample s of observations; each observation i ∈ s has associated a pair ( X i , Y i ), where each X i : [0 , T ] → R belongs to the set X of Riemann integrable functions in the time interval [0, T ]. Furthermore, Y i ∈ {−1 , +1 } denotes the class label for the observation i . Our goal is to find a classification rule to infer the class Y of a new functional observation X ∈ X . The well-known technique SVM considers a kernel K : X ×X → R , ( Cristianini & Shawe-Taylor, 20 0 0; Rossi & Villa, 20 06, 20 08 ), and builds, from a sample s , nonlinear classifiers by means of a score ˆ Y (X) of the form: ˆ Y (X ) =  i ∈ s αi Y i K(X, X i ) , X ∈ X , (1) yielding the following classification rule: a functional observation X ∈ X is assigned to class +1 if and only if ˆ Y (X) > β, where βis a given threshold value. Here the values αi , i ∈ s , are obtained as the optimal solution of the following optimization problem: ⎧ ⎪ ⎨ ⎪ ⎩ max α i ∈ s αi −1 2  i,j∈ s αi αj Y i Y j K(X i , X j ) s.t.  i ∈ s αi Y i = 0 αi ∈ [0 , C] , i ∈ s, (2) for a scalar regularization parameter C to be tuned, usually by k - fold cross-validation with a grid search on a sufficiently large interval. Many types of kernels for data in R d are proposed in the literature, e.g., the linear kernel, the polynomial kernel, or the Gaussian (RBF) kernel, given by: K(X i , X j ) = exp  − d  t=1 (X it −X jt ) 2 ω  , X i , X j ∈ R d (3) where ω is a scalar bandwidth to be tuned, ( Carrizosa et al., 2014; Carrizosa & Romero Morales, 2013; Cristianini & ShaweTaylor, 20 0 0; Hofmann, Schölkopf, & Smola, 2008; Keerthi & Lin, 2003 ). In this paper, for simplicity, we only focus on the Gaussian kernel, one of the most used and effective kernels, which will be used in what follows. The expression (3) of the Gaussian kernel for data in R d has been generalized to a Gaussian kernel for functional data, e.g., Kadri, Duflos, Preux, Canu, and Davy (2010) and Wang and Yao (2015) . Nevertheless, in these papers, the associated bandwidth is always considered to be a scalar value. In our proposal we extend the fixed scalar bandwidth parameter ω in an RBF kernel to a functional bandwidth, ω( t ), that varies along the range of the functional data, (4) : K(X i , X j ) = exp − T 0 (X i (t) −X j (t)) 2 ω(t) dt (4) Throughout this paper, we assume that ω in (4) is a non-negative Riemann integrable function in [0, T ], and thus K is well-defined. It is worth mentioning that the simplest extension from the kernel with vector data (3) to the kernel with functional data (4) would be to consider ω( t ) as a constant function, as in Kadri et al. (2010) and Wang and Yao (2015) . Nevertheless, the main contribution of this paper is to consider such bandwidth as a function which adapts to the structure and shape of the data and may lead to better insight and classification rates. More specifically, making ω depend on t allows us to identify those subintervals in [0, T ] which are critical for classification, namely, those for which ω( t ) takes highest values. Example 2.1. As an illustration, let us study the regions data set ( Martín-Barragán et al., 2014 ), in which the daily temperature has been measured along a year in each of 35 Canadian weather stations. Two groups can be distinguished: Atlantic climate (label -1), with 15 records, versus the rest of climates (label 1), with 20 records. Our objective is to discriminate between both classes. Fig. 1 depicts the 15 curves in the interval [1, 365] corresponding to the Atlantic climate, in solid black line, and the 20 curves corresponding to the rest of climates, in dashed red line, with the data measured every single day. This is, by nature, a Functional Data classification problem. However, it may be considered as a classic classification problem with 15 + 20 records in R d , d = 365 (the number of time instants at which the temperature has been actually recorded), and thus one can apply the classic SVM in the form (3) for some ω to be tuned. Observe that this model is the same as model (4) with ω(t) = ω, ∀ t ∈ [0 , T ] with T = 365 , (5) R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 197 −30 −20 −10 0 10 20 Months Temperature JFMAMJ J ASOND Atlantic rest Fig. 1. regions data set. Table 1 Confusion matrix with ω as in (5) . Label −1 Label 1 Label −1 51.42% 5.71% Label 1 11.42% 31.42% Table 2 Confusion matrix with ω as in (6) . Label −1 Label 1 Label −1 54.28% 2.85% Label 1 8.57% 34.28% and the integral evaluated numerically in the grid of time instants where the temperature is measured. Using SVM with a constant ω( t ) as in (5) leads to a classifier with the out-of-sample confusion matrix shown in Table 1 . Now, let us consider the very same RBF model with a functional bandwidth ω( t ) of the form ω(t) = ω 1 , if 0 ≤t ≤τ1 ω 2 , if τ1 < t ≤365 , (6) where ω 1 , ω 2 , τ1 are parameters to be tuned using the techniques described in this paper. In other words, with the bandwidth in (6) we split the interval [0, T ] into two pieces, giving different weights to each time interval. The SVM classifier obtained this way leads to the out-of-sample confusion matrix in Table 2 . Comparing Tables 1 and 2 we can see that the traditional SVM yields an accuracy of 82.84%. On the other hand, our SVM with the very same RBF kernel but using a functional parameter of the form (6) yields an accuracy of 88.56%, instead. Regarding the interpretability of the results, Figs. 2 and 3 show the boxplots of the values of the bandwidth ω as in (5) , and the values of ω 1 , ω 2 and τ1 , as in (6) . The single-bandwidth approach gives the same importance to all the months of the year with the majority of the bandwidth values between 50 and 150. In contrast, our functional-bandwidth methodology with two different pieces proposed to divide the whole year into two parts, before and after summer (months of June and July), see Fig. 3 . Moreover, according to the values of ω 1 and ω 2 , in order to get good classification predictions, we should focus on the second half year and give more importance to the second part, i.e., the autumn and first months of winter, which coincides to the time instants when the temperature begins to decrease. The previous illustrative example demonstrates that even a simple functional bandwidth such as (6) may yield improvements in accuracy. Such improvement is a consequence of the adequate choice of the parameter τ1 , which combined with good values of ω 1 and ω 2 allow us to identify the suitable intervals for classification. The functional bandwidth parameter ω( t ) gives more flexibility, which should result in greater precision. For instance, it may be chosen in the class of piecewise constant non-negative functions in [0, T ] with H pieces, i.e., one can naturally assume that ω( t ) has the form (7) ω(t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ω 1 , if 0 ≤t ≤τ1 ω 2 , if τ1 < t ≤τ2 ··· ω h , if τh −1 < t ≤τh ··· ω H , if τH−1 < t ≤T (7) 198 R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 Fig. 2. (a) and (b) Show the bandwidth values for the regions data set when ω has the form of (5) and (6) , respectively. Fig. 3. Time instant results for the regions data set with ω as in (6) . where ω 1 , . . . , ω H ≥0 and 0 ≤τ1 ≤···≤τH−1 ≤T are parameters to be tuned. Instead of piecewise constant functions, one could consider ω( t ) belonging to the class of polynomials of degree H which are non-negative in [0, T ], the class of piecewise polynomial functions non-negative in [0, T ], or the non-negative splines, ( De Boor, 1978; Friedman, Hastie, & Tibshirani, 2001b ). The use of functional parameters in the kernel may lead to significant improvements in the accuracy, as demonstrated in our numerical experiments. The price to pay for obtaining such gains in the accuracy is the fact that tuning the functional parameters calls for using more sophisticated optimization procedures. In Section 3 we detail how the underlying optimization problem for tuning ω( t ) is solved. 3. Optimal selection of the functional bandwidth In this Section, a detailed study of the mathematical formulation of the (functional) parameter tuning in SVM classification is presented. Section 3.1 explains how to formulate the optimization problem involved and how to solve it. In Section 3.2 a nested R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 199 heuristic to address the tuning problem more efficiently is described. In this way, we exploit the fact that the bandwidths considered are elements of a nested family of kernels. Section 3.3 details how to choose the number of pieces, H , of the functional bandwidth. 3.1. Problem formulation and optimization Parameter tuning in the classification of functional data with SVM implies the optimal choice of two very different elements: the scalar regularization parameter C in (2) , and the kernel K in (4) through ω( t ). The problem of finding the best function ω( t ) in (4) , is not tractable as a rule in its full generality. Hence, we restrict our attention to certain classes of functions parameterized by a vector θbelonging to a certain set , i.e., ω is expressed as ω( t , θ), and the choice of the function ω is equivalent to choosing the parameters θ. Example 3.1. For the bandwidth given in (7) , one would have that θ= (ω 1 , . . . , ω H , τ1 , . . . , τH−1 ) , and = { (ω 1 , . . . , ω H , τ1 , . . . , τH−1 ) : ω h ≥0 , ∀ h, τh ∈ [0 , T ] , h = 1 , . . . , H −1 , τ1 ≤···, ≤τH−1 } . For convenience, we consider τ0 = 0 and τH = T . In principle, in order to find the optimal values of the parameters, C , and θ, a strategy based on a grid search on both parameters could be applied. Given a set of predefined pairs of values ( C , θ), one first solves (2) to obtain the coefficients αof the score function (1) , and then the corresponding accuracy associated to that pair is computed. However, this approach may become too time consuming, and thus a more sophisticated heuristic is proposed in these lines. We propose to follow the standard grid approach to optimize C . Nevertheless, when seeking the parameters θand α we propose to solve a bilevel problem where some measure of the quality of θis maximized for the αprovided by the SVM classifier, i.e., for the αsolving (2) . Many criteria can be chosen to guide the choice of parameter θ. One may, for instance, minimize the misclassification rate, which is the default approach for tuning the parameter C . However, the misclassification rate has a discrete nature that would prevent us from using continuous optimization techniques, and, in particular, from gradient-based methods. Instead, we propose to maximize the Pearson correlation, R , between the class label Y i of the functional data X i and the score, ˆ Y (X i , θ, α) in (1) , where all the variables, including the time instants, are treated as continuous variables. Other references in the literature, such as Blanquero et al. (2017) and Jiménez-Cordero and Maldonado (2018) , have previously used with excellent results the Pearson correlation coefficient. Despite the fact that, when using the Pearson correlation coefficient as a surrogate of accuracy, a linear relationship between the binary label, Y ∈ {−1 , 1 } , and the real-valued score, ˆ Y ∈ R , is implicitly assumed, this coefficient is very fast to compute, and even more important, it also allows us to use gradientbased methodologies since its optimization amounts to solving a continuous optimization problem. It is very well known that building a classifier and evaluating its performance over the same data set leads to overfitting. In such a case, the model fits the data set too well but performs poorly in unseen data. On top of that, the classifier can depend on parameters that must be tuned, usually done by performing a grid search in a suitable range of values. The usual way to avoid overfitting in this general situation is to split the data set, perhaps within a k -fold cross-validation framework, in three parts, the so-called training, validation, and test samples. For a given choice of the parameters, the first two ones are used to build the model and estimate its performance, respectively; once the best parameters have been chosen, the final model is tested on the last sample. In our case, we take this idea further by creating four independent samples, due to the structure of our resolution method. First, the data set is divided into k folds. Second, k −1 folds are again split into three samples named s 1 , s 2 , and s 3 , while the remaining fold is denoted by s 4 . Samples s 1 and s 2 play the role of training samples, whereas s 3 and s 4 form the validation and testing sets, respectively, as will be detailed next. The first independent sample s 1 is employed for the resolution of Problem (2) , that is the classic SVM formulation, to obtain a classification rule by means of α, for fixed parameters θand C . The second independent sample s 2 is used to measure the quality of parameters θ, i.e., it is used to calculate R ((Y i , ˆ Y (X i , θ, α)) i ∈ s 2 ) , the correlation between the class labels and the scores. To find the regularization parameter C , we measure the accuracy in the sample s 3 for all the different possible values of C in the grid, and we keep the C providing the largest accuracy. Finally, the accuracy in the independent sample s 4 is measured. After all these considerations, for fixed C , the bilevel problem can be expressed as: ⎧ ⎨ ⎩ max θ,αR ((Y i , ˆ Y (X i , θ, α)) i ∈ s 2 ) s.t. αsolves (2) in s 1 θ∈  (8) Note that we have emphasized the dependence of the score ˆ Y on θand αby including them in the notation. In the cases where the values of the parameters in θ, or the classification coefficients α are clear enough, we will omit them for the sake of simplicity. Problem (8) is a nonlinear bilevel optimization problem, which can be handled with off-the-shelf strategies, as those described in Colson, Marcotte, and Savard (2007) . These techniques are, however, rather expensive. Recall that (8) is only a surrogate of our real problem. Hence, instead of the above-mentioned standard methodologies, we next propose an alternating approach for which only a few iterations will be carried out. Firstly, in the first step of our alternating approach, for fixed parameters θand C , a classification rule is obtained by means of αsolving Problem (2) , that is, the classic SVM formulation. Problem (2) is a concave quadratic maximization problem, which can be solved by standard local search optimizers or specific routines, as in Ferris and Munson (2004) ; Richtárik and Takᡠc (2016) . Secondly, in the second step, for fixed αand C , θis chosen by solving: max θ∈ R ((Y i , ˆ Y (X i , θ)) i ∈ s 2 ) (9) Problem (9) is a continuous optimization problem which is solved by using standard local search techniques within a multi-start strategy. The alternating procedure will alternate these two steps until some stopping criterion is met. Suitable values for θand α will be obtained by this procedure for a specific value of the regularization parameter C . The value of C will be chosen by a grid search, as commonly done in standard SVM. This means that, for every value of C in a given grid, we measure the accuracy in s 3 of the classification rule obtained with the best θand αfound as solutions of Problem (9) . The C with the largest accuracy in s 3 will be chosen. Finally, we estimate the correct classification rate using the fourth independent sample, s 4 . A pseudocode of the heuristic is outlined in Algorithm 1 , and in Section 3.2 , we detail an extension to more complex models by means of a nested heuristic, described above. 3.2. Optimization enhancement. A nested optimization When the dimension of θis high, the approach described in Section 3.1 may be time-consuming. The main reason is that, on top of the grid search needed for C , Problem (9) may have many 200 R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 Algorithm 1 Heuristic for parameter tuning. Input: H •Randomly split the sample s into s 1 , s 2 , s 3 and s 4 . for C in the grid do Alternating Procedure repeat 1. Fixed θ, calculate the parameters αof the SVM classifier by solving Problem (2) in s 1 . 2. Fixed α, calculate θby solving Problem (9) in s 2 . until stopping criteria •Evaluate the accuracy in the sample s 3 with C fixed. end for •The optimal value of C is the one with the best accuracy in s 3 . The optimal values of αand θare the ones associated to the optimal parameter C. Output: optimal parameters C and θ, optimal classification coefficients α, and the corresponding accuracy estimated from s 4 . local minima, and therefore multiple local searches are required to find a good solution. The success of the method would be improved if, instead of a random multi-start, a more intelligent search strategy were possible. This is the case, for instance, for models of the bandwidth parameter ω( t , θ) that can be plugged into a sequence of models of increasing complexity. Thus the optimal solution obtained in the simple model can be used as a starting solution in the following more complex case. The above-explained methodology can be easily embedded in a nested heuristic for SVM parameter tuning, Carrizosa et al. (2014) , in which a nested structure of kernels is assumed. More precisely, given a family of kernel functions, we construct a series of nested kernel models with their associated parameters, or equivalently, a series of H nested functional bandwidths ω (1) ( t , θ(1) ) ≺≺ ω ( H ) ( t , θ( H ) ). By ω (h ) (t, θ(h ) ) ≺ω (h +1) (t, θ(h +1) ) we denote that the bandwidth ω ( h ) ( t , θ( h ) ) has parameters which are part of the parameters of the bandwidth ω (h +1) (t, θ(h +1) ) . Example 3.2. Consider the family of piecewise constant functions with 3 pieces, in (7) . We have that ω (1) (t, θ(1) ) = ω 1 , with θ(1) = ω 1 , ω (2) (t, θ(2) ) = ω 1 I [0 ,τ1 ] + ω 2 I (τ1 ,T ] , with θ(2) = (ω 1 , ω 2 , τ1 ) , and finally ω (3) (t, θ(3) ) = ω 1 I [0 ,τ1 ] + ω 2 I (τ1 ,τ2 ] + ω 3 I (τ2 ,T ] , with θ(3) = (ω 1 , ω 2 , ω 3 , τ1 , τ2 ) . Here I [ r,r  ] denotes the indicator function, i.e., the function which is equal to 1 in the interval [ r , r  ] and 0 otherwise. The idea of using nested models is to take advantage of the easy-to-tune structure of the elementary models and consider them as a simplification of the complex models. When solving Problem (8) for ω ( H ) ( t , θ( H ) ) we will use a sequential approach where the (suboptimal) solution obtained when using ω ( h ) ( t , θ( h ) ), will be used as an initial solution of Problem (8) with ω (h +1) (t, θ(h +1) ) . Example 3.3. For the bandwidth given in (7) , once we have obtained the (suboptimal) solution of ω ( h ) ( t , θ( h ) ) by θopt (h ) = (ω opt 1 , . . . , ω opt h , τopt 1 , . . . , τopt h −1 ) , we randomly select an interval [ τ −1 , τ ) and split it into two pieces by its midpoint, assigning the same bandwidth value to such two new pieces. In other words, the initial point of the parameters in the level h + 1 turns out to be θ(h +1) = (ω opt 1 , . . . , ω opt  −1 , ω opt  , ω opt  , ω opt  +1 , . . . , ω opt h , τopt 1 , . . . , τopt  −1 , τopt  + τopt  −1 2 , τopt  , . . . , τopt h ) . The pseudocode of the nested algorithm defined in Section 3.1 , is shown in Algorithm 2 . Algorithm 2 Nested heuristic for parameter tuning. Input: H, nested functional bandwidths ω (1) (t, θ(1) ) ≺···≺ ω (H) (t, θ(H) ) . •Randomly split the sample s into s 1 , s 2 , s 3 and s 4 . for C in the grid do Initialization: •h := 1 . •Randomly select an initial solution θ(h ) ∈ (h ) . •Set θ:= θ(h ) while h ≤H do 1. Using samples s 1 and s 2 , run the Alternating Procedure of Algorithm 1 for ω(t, θ(h ) ) , starting from θand yielding θopt (h ) = (ω opt 1 , . . . , ω opt h , τopt 1 , . . . , τopt h −1 ) as solution. 2. Randomly select  ∈ { 1 , 2 , . . . , h } . 3. Set θ:= (ω opt 1 , . . . , ω opt  −1 , ω opt  , ω opt  , ω opt  +1 , . . . , ω opt h , τopt 1 , . . . , τopt  −1 , τopt  + τopt  −1 2 , τopt  , . . . , τopt h −1 ) and h := h + 1 . 4. Evaluate the accuracy in the sample s 3 with C fixed. end while end for •For h fixed, the optimal value of C is the one with the best accuracy in s 3 . The optimal values of αand θ(h ) are the ones associated to the optimal parameter C. Output: optimal parameters C, θopt (h ) , ∀ h , the associated classification coefficients α, and the accuracy estimated from s 4 . 3.3. Choice of the number of pieces, H Thus far we have assumed that the number H of pieces is given as input in the problem, and hence the results are dependent on H . The larger H is, the better the accuracy (in the training sample) since more flexibility is added to the model. However, if a too large value of H is chosen, the number of parameters involved in the problem increases considerably, and this may deteriorate the accuracy in the test sample. Therefore, it is sensible to define a strategy to determine the best H . In this respect, standard criteria, such as BIC, AIC or ICL, ( Akaike, 1974; Biernacki, Celeux, & Govaert, 20 0 0; Schwarz, 1978 ) can be applied in the SVM context, as done in Claeskens, Croux, and Kerckhoven (2008) for instance. They proposed two new information criteria which are inspired, but not equal to AIC and BIC, with the aim of giving consistent selection criteria without much additional computational costs. In contrast, in this paper, we propose to keep the parameter H with the largest accuracy on the validation sample s 3 . 4. Numerical experiments This section details the experiments performed ( Section 4.1 ) and the main characteristics of the data bases here considered ( Section 4.2 ). Finally, Section 4.3 is devoted to the computational results obtained. 4.1. Description of the experiments In this section, a detailed description of the experiments carried out to test our methodology is made. To obtain stable estimates, k - fold cross-validation has been used to evaluate the performance of the algorithm on different data sets. The number k of folds varies depending on the size of the database. For small databases, k is equal to the number of observations, i.e., we performed leave-oneout, whilst for large databases we take k = 10 . A database is con- R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 201 Table 3 Real data description summary. #Records #Points measurements #Records label −1 #Records label + 1 ECG 200 96 67 133 growth 93 31 54 39 gun 200 96 100 100 MCO 89 360 44 45 phoneme 200 150 100 100 phoneme_large 1717 256 1022 695 rain 35 365 15 20 regions 35 365 20 15 synthetic_magnitude 150 100 75 75 tecator 215 100 77 138 wine 111 234 54 57 yoga 306 426 150 156 sidered small here if and only if it has less than 100 observations. See Table 3 . Algorithm 2 is run k times, one per fold. Each time, the division into four independent samples s 1 , s 2 , s 3 , and s 4 is done as explained in Section 3.1 . The number of runs of the multi-start local search optimization method is set to five. The algorithm is run until the maximum number of iterations reached to ten, or when the difference between the objective values in two consecutive iterations is less than 10 −5 . The functional bandwidth ω( t , θ) is the piecewise constant function in (7) with H = 8 . The regularization parameter C varies in the set { 2 −10 , . . . , 2 10 } . The parameters θ( h ) are in the set (h ) = { (ω 1 , . . . , ω h , τ1 , . . . , τh −1 ) : ω  ≥2 −4 ,  = 1 , . . . , h, 0 ≤τ1 ≤. . . ≤τh −1 ≤T } , ∀ h = 1 , . . . , 8 . For comparison purposes, apart from the standard SVM, i.e., our approach with H = 1 , we have run three supervised classification methods for functional data, available at the fda.usc library of R ( Febrero-Bande & Oviedo de la Fuente, 2012 ), namely classif.depth , classif.kernel , classif.knn with the default parameters. In order to obtain a fair comparison, the accuracy obtained is estimated on the very same testing sample s 4 used in our approach. Our algorithm is coded in R and is carried out on a cluster with 2 terabyte of RAM memory at 6.2 TFlops, running CentOS Linux 7.3. The code is available upon request. 4.2. Description of the data sets Our methodology has been tested in 12 benchmark data sets, widely used in the functional data classification literature, namely, ECG , ( Chen et al., 2015; Xing, Pei, & Philip, 2009 ), growth , ( Cuevas et al., 2007; Muñoz & González, 2010; Torrecilla Noguerales, 2015 ), gun , ( Chen et al., 2015; Xing et al., 2009 ), MCO , ( Baíllo et al., 2011; Cuevas, Febrero, & Fraiman, 2006; Ruiz-Meana et al., 2003 ) and Online companion of ( Carrizosa et al., 2014 ), phoneme , ( Ferraty & Vieu, 2006; Muñoz & González, 2010; Rossi & Villa, 2006; Torrecilla Noguerales, 2015 ), phoneme_large , ( Berrendero et al., 2016; Delaigle & Hall, 2012; Friedman, Hastie, & Tibshirani, 20 01a; 20 01b ), rain , ( Martín-Barragán et al., 2014 ), regions , ( Martín-Barragán et al., 2014 ), synthetic_magnitude , model 3 of ( López-Pintado & Romo, 2009 ), tecator , ( Ferraty & Vieu, 2006; Martín-Barragán et al., 2014; Rossi & Villa, 2006; Torrecilla Noguerales, 2015 ), wine , ( Chen et al., 2015 ) and yoga , ( Wei, 20 06; Wei & Keogh, 20 06 ). Note that the data set phoneme is used as described in the fda.usc library, ( Febrero-Bande & Oviedo de la Fuente, 2012 ), of R . Table 3 summarizes the data sets description, which gives the overall number of records, the number of time measurements, and the number of records of each class. A plot with a sample of 10 instances of each data set is shown in Fig. 4 , depicting in solid black line the observations with label −1 and in dashed red line the records with label 1. The number of folds is determined by leave-one-out in the data sets growth , MCO , rain , and regions , and with 10 −fold cross-validation in the remaining databases. 4.3. Results We provide the boxplots of the accuracy measured on s 4 from h = 1 to h = 8 for the different folds in the k -fold accuracy estimation procedure. Boxplots are not very informative for small data sets, for which leave-one-out is performed. Indeed, for each fold either one obtains an accuracy of 0% or 100%, since either the testing observation is wrongly or correctly classified. For this reason, only the boxplots of the largest data sets are depicted. See Fig. 6 . Moreover, the exact values of the average accuracy and its standard deviation, as well as the corresponding values for the three fda.usc library methods considered in Section 4.1 , are also presented in Table 4 for the sake of comparison. The four gray columns correspond to the four methods we are comparing with, denoted as depth , kernel , knn and classic SVM, h = 1 . Finally, last column of Table 4 gives the best number of pieces chosen, according to the strategy explained in Section 3.3 . We have highlighted in bold in Table 4 the maximum of the accuracy values for h = 2 , . . . , 8 which are equal or greater than any of the four methods. In general, our method for h = 2 , . . . , 8 is better than the four comparative approaches in the data sets growth , MCO , phoneme , phoneme_large , and regions . This improvement may be produced by the shape of the curves. The different class labels seem to be easy to identify depending on the time subinterval, and therefore our strategy makes easier such separation. Observe for instance, the growth data set, in which the two classes have a different pattern around the time instant 15. Moreover, it is seen that the improvement in the accuracy strongly depends on the data set considered. Indeed, no improvement is seen for the databases gun , rain , and tecator when comparing our methodology with h = 1 and h ≥2. However, for some of the values h ≥2 the accuracy obtained in gun is better than that provided by depth . The results of our approach in the database rain are always better than the ones provided by the three fda.usc methods. In contrast, such three methods should be applied if the tecator data set is studied. In the databases ECG , growth , phoneme_large and yoga there is a minor improvement (about a 0.5%) when comparing the classic SVM with our approach for h ≥2. Such improvement also holds in the ECG data set when comparing with the depth method. The accuracy value obtained in phoneme_large with our approach when h = 4 pieces are optimally chosen is better than all the three fda.usc methods. Analogous conclusions are obtained in the yoga data set. A considerably larger accuracy is obtained in databases MCO , phoneme , regions , synthetic_magnitude , and wine when solving the problem with h ≥2 than when solving with h = 1 , i.e., the classic SVM. In some cases such improvement yields 202 R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 Fig. 4. Sample of functional data in the real data sets analyzed. R. Blanquero, E. Carrizosa and A. Jiménez-Cordero et al. / European Journal of Operational Research 275 (2019) 195–207 203 Fig. 4. Continued