title: How to tune the RBF SVM hyperparameters?: An empirical evaluation of 18 search algorithms
authors: Jacques Wainer Pablo Fonseca
text: 
  
  
  
  
  
  
  
  How to tune the RBF SVM hyperparameters?: An empirical evaluation
    of 18 search algorithms
      Jacques Wainer^1Corresponding author: wainer@ic.unicamp.br Pablo Fonseca^2
      ^1Computing Institute
  
   University of Campinas
  
   Campinas, SP, 13083-852, Brazil
  
  ^2Facultad de Ciencias y Filosofía
  
    Universidad Peruana Cayetano Heredia
  
  Lima, Peru
  ============================================================================================================================================================================
  
  
  
  
  
  
  
  SVM with an RBF kernel is usually one of the best classification algorithms for most data sets, but it is important to tune the two hyperparameters C and γ to the data itself. In general, the selection of the hyperparameters is a non-convex optimization problem and thus many algorithms have been proposed to solve it, among them: grid search, random search, Bayesian optimization, simulated annealing, particle swarm optimization, Nelder Mead, and others. There have also been proposals to decouple the selection of γ and C. We empirically compare 18 of these proposed search algorithms (with different parameterizations for a total of 47 combinations) on 115 real-life binary data sets. We find (among other things) that trees of Parzen estimators and particle swarm optimization select better hyperparameters with only a slight increase in computation time with respect to a grid search with the same number of evaluations. We also find that spending too much computational effort searching the hyperparameters will not likely result in better performance for future data and that there are no significant differences among the different procedures to select the best set of hyperparameters when more than one is found by the search algorithms.
  
  Keywords: Hyperparameters; SVM; grid search; random search; non-convex
  optimization algorithms
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Support vector machines (SVM) and in particular SVM with the RBF kernel is a powerful classification technique, ranking amont the best two or three classifiers for a large set of real life datasets <cit.>. SVM with RBF kernel requires one to select the value of two hyperparameters C and γ before computing the SVM itself. There have been many proposals in the literature on how to search for the best pair of hyperparameters for a particular data set, among them general optimization procedures such as grid search, random search, simulated annealing, Bayesian optimization, and many others.  The main goal of this paper is to compare many of these proposals on a large number of real-life binary data sets. The 18 families of search algorithms that are compared in this research are discussed in Section <ref>, including their different parameterizations. In total, we tested 47 different combinations of algorithms and parameterizations.
  
  SVM, in general, is formulated as a minimization problem:
  
      min     1/2 ||w||^2 + C∑^n_i=1ξ_i 
      subject to     y_i(w ·ϕ(x_i)+b)≥ 1-ξ_i 
      with     ξ_i ≥ 0
  
  The parameter C is a regularization factor and it is one of
  the SVM hyperparameters: this constant must be set before solving the
  minimization problem. 
  
  ϕ(x) is a non-linear transformation that takes the data into a high dimensional space, sometimes called a Reproducing Kernel Hilbert Space and sometimes called the feature space. Fortunately, ϕ(x) does not need to be computed explicitly. The dual formulation of the SVM minimization does not need ϕ(x) but only the inner product of the transformation of two data points ϕ(x_i)^T ϕ(x_j). The dual form of the SVM minimization is:
  
  
      min     1/2∑_ij y_i α_i y_j α_j
              K(x_i, x_j) - ∑_i α_i 
      subject to     ∑_i y_i α_i = 0 
      with     0 ≤α_i ≤ C
  
  where
  K(x_i,x_j)=ϕ(x_i)^T ϕ(x_j) is a kernel function
  that represents the ϕ  mapping.
  
  The Gaussian Radial Basis Function (RBF) kernel is the chosen kernel function discussed in this paper, and it is expressed as:
  
      K(x_i,x_j) = exp(-γ||x_i-x_j||^2).
  
  
  The constant γ is also an RBF SVM hyperparameter: a value that must be set before the optimization in Equation <ref> is computed.  To solve a particular classification problem, one must choose a particular pair of C and γ that is the most appropriate to that problem.
  
  Let us define (G' |G , C, γ) as the accuracy of an RBF SVM on the data set G' when it has been trained on the data set G with hyperparameters C and γ.
  
  
  Let us assume a given data set G as an i.i.d. sample from an unknown distribution D. One would like to select the best pairs of hyperparameters as:
  
      C_best, γ_best =   _C, γ _g ∼ D [(g|G,
      C, γ) ]
  
  where _g ∼ D is the expectation taken over a random variable g that represents future data with the same distribution as the data G.
  
  One can consider the right side of Equation <ref> as the definition of a function of two parameters, C and γ, which we will call the true response surface of the SVM for the data set G. We will represent the true response surface as _G(C,γ).
  
  Of course, one does not have access to the new data g or the underlying distribution D, and thus one cannot compute _g [(g|G, C, γ) ]. One solution is to compute an estimate of that value using some form of resampling of the original data set G. Resampling will divide the data set G into one or more pairs of data sets TR_i and TE_i, such that TR_i ∩ TE_i = ∅ and TR_i ∪ TE_i ⊆ G. K-fold cross-validation, bootstrapping, and leave-one-out are examples of resampling procedures.
  
  The estimate of _g [(g|G, C_1, γ_1) ] for some fixed C_1 and γ_1 is the mean of the accuracy of training the SVM (with C_1 and γ_1 as hyperparameters) on each of the TR_i (the training set) and testing on the corresponding TE_i (the test set). In this paper, we will use 5-fold cross-validation as the resampling procedure.
  
  
  We will define an approximation to the function _G (C,γ)
  as:
  
      (G, C, γ) = _G (C,γ) =   _i  ( TE_i | TR_i, C, γ)
  
  
  For simplicity, we will call this the response surface (without the qualifier true) of the SVM for the data set G and for a particular resampling choice of pairs { TR_i, TE_i }.
  
  Selecting the best hyperparameters for the data set G is finding the maximum of the response surface _G. Or, to be more consonant with the optimization terminology, one can think of minimizing the error surface 1- _G.
  
  In general, the response surface is not concave (or the error surface is not convex) and thus there may be many local maxima, which requires optimization algorithms that are not just based on gradient ascent. There have been many proposals of algorithms that will find hyperparameters with higher values of _G. The first goal of this research is to compare many of those algorithms.
  
  One must remember that the aim of hyperparameter search is to find a pair of hyperparameters that maximizes the true response surface _G, that is, the pair that will result in the classifier with the highest accuracy for future data. The _G surface is an approximation of the _G surface, and thus there will likely be differences in the value of _G(C_1,γ_1) and _G(C_1,γ_1). The second goal of this paper is to have some understanding if it is worth spending the computational effort on algorithms that will likely select better hyperparameters in  but may not achieve a better solution in . We will estimate the correlation of improving the selection in  and the expected gain in accuracy on future data.
  
  
  The _G is also a step function - small changes in C and γ will not change the accuracy of the resulting SVM on a finite set of data. Therefore, it is likely that any search algorithm for pairs of C and γ will produce a set of pairs of hyperparameters with the same value of _G.
  
  If more than one pair of hyperparameters is selected there is a reasonable argument to choose the pair with the smallest C and the smallest γ. For example, the tutorial page for hyperparameters  in SVM in sklearn <cit.> makes that argument. Both hyperparameters work as inverse regularization terms. A large C will place emphasis on lowering the number of support vectors since each one of them contributes to the ∑^n_i=1ξ_i cost in the optimization. A lower C will allow more support vectors, resulting in larger margins.
  
  The γ parameter controls how fast the “influence” of a point decreases with distance. The kernel value for two points will decrease as γ increases. As γ increases, the decision surfaces become more "curvy" and fit closely to the training data. A smaller γ will generate decision surfaces that are flatter, and thus a simpler model.
  
  Thus, in the case of search algorithms that return a set of best hyperparameters, one has to perform a further selection procedure. We call this post-search selection.  The third goal of this paper is to compare some different post-search selection procedures regarding the accuracy of future data.  
  
  
  
  
  The fourth goal of the paper is to get some insight into the distributions of the best C and γ selected by the many algorithms, whether these values are concentrated on some region of the hyperparameter space.
  
  
  
  
  § HYPERPARAMETER SEARCH ALGORITHMS
  
  
  Broadly speaking, there are two different approaches to the SVM hyperparameter search. The first approach considers the hyperparameter search as a general non-convex optimization problem and uses general algorithms to solve it.  Among such general optimization algorithms are grid search, random search, uniform design <cit.>, gradient-free minimization algorithms such as Nelder-Mead simplex <cit.>, genetic algorithms, particle swarm optimization, ant colony optimization, Bayesian optimization <cit.>, simulated annealing <cit.>. All these alternatives have been proposed as search algorithms for the SVM hyperparameters.
  
  The second general approach considers the mathematical properties of the SVM formulation itself, and propose particular search strategies to solve specifically the SVM hyperparameter optimization problem. This approach can be further divided into two broad directions. The first one defines smooth approximations to the estimated error of the classifier (usually the leave-one-out error) and uses gradient descent algorithms on these approximations <cit.>. We will not test these algorithms in this research.
  
  The second direction tries to decouple the choice of γ from the choice of C (the more frequent approach) so that γ is selected first (sometimes using a grid search, sometimes not) based on some measure of the data, and only then the C is chosen (usually by grid search).
  
  
  
  
  
  
   §.§ General optimization approaches
  
  
  We will discuss some of the general optimization algorithms for searching for the SVM hyperparameters.
  
  
  
    * grid: In the grid search the log_2 C ×log_2 γ space is divided into an equally spaced grid and only the pairs in the grid points are tested.
  
    We used the following constraints on C and γ not only for the grid search but for all algorithms:
  
      2^-5≤    C ≤ 2^15
      
         2^-15≤   γ≤ 2^3
  
  which are the limits of the grid search suggested by the LibSVM package <cit.>, a very popular SVM solver. The grid search and all other algorithms perform the search in the space (and not in the C ×γ space).
  
  As we will discuss further in Section <ref>, we are interested in comparing the quality of the different algorithms when controlling for the number of times the optimizer evaluates  for some point in the space. We call this number N. Notice that since we are using 5-fold to evaluate , the SVM learning and evaluation methods will be called 5*N times.
  
  We tested the grid procedure for N = 25, N=100 and N=400, and in each case, we divided both the log C and the logγ dimensions into an equal number of intervals (5, 10, and 20 respectively)
  
  
    * ud: Uniform design is a space-filling method <cit.> used in the design of computer experiments. It approximates a uniform distribution over a space of parameters (or hyperparameters in our case), given the number of points one needs to sample from that space. The sampling is deterministic and it chooses points that minimize a discrepancy measure in relation to the uniform distribution. Uniform design has been proposed as a method to select SVM hyperparameters <cit.> .
  
  To generate the uniform design tables we used the  function from the  package <cit.>.  We tested ud with N = { 25, 100, 400}.
  
  
    * rand: N random points are uniformly sampled from the space within the bounds of Equation <ref>.
  
  
  
    * normrand: Randomly select log C from a normal distribution with mean 5 and standard deviation 5 and logγ from a normal distribution with mean -5 and standard deviation 5.
  
    We tested the other algorithms discussed in this section on 115 datasets (see Section <ref>) and collected the C and γ selected by them. The distribution of the selected log C and logγ follow a Gaussian distribution (see Section <ref>). This prompted the realization that a random sampling using these Gaussian distributions could achieve, in general, better results than the rand search discussed above.
   
    We tested normrand with N = { 25,100, 400}. But we must warn that the analysis of the results for this algorithm has a risk of overfitting - the algorithm was tuned using all the datasets available and we are measuring their quality on the same set of datasets.
  
   
  
  
    * gridhier: Hierarchical grid. A grid search is performed to select the point x as the one with the highest accuracy. A second grid search centered on x and bounded by the closest neighbours of x in the first grid is performed, and the highest accuracy point among x and from the second search is selected.
  
    We tested the hierarchical procedure both for grid search with 25+25, that is, with N =25 in the first level, followed by N=25 in the second level, and with 100+100
  
  
    * udhier: Hierarchical uniform design. Similar to the gridhier, but the uniform design search is performed at both levels. Tested with N = 25+25 and 100+100.
  
  
  
    * nelder: Nelder-Mead simplex algorithm <cit.> is a well-known derivate-free optimization algorithm. The algorithm maintains a simplex of points in the parameter space (in our case of 2 dimensions the simplex is a triangle) and updates the simplex point with the worst evaluation with a new point symmetric to it on the other side of the centroid of the remaining points. The algorithm converges when either the simplex is too small or when the improvement on the evaluation on changing the worse simplex vertex is below a threshold. Nelder Mead has been suggested as a method for selecting SVM hyperparameters  <cit.>.
  
    We used the  function of the  package <cit.> with the initial point in the center of the space. There are many algorithms to generate an initial simplex from a single initialization point <cit.>; unfortunately, neither the  R package nor the underlying  C implementation documents which algorithm is used. 
  
  
  
  
    We used the bounded version of the Nelder Mead algorithm with the constraints in Equation <ref> We tested Nelder-Mead with N={25,100, 400}. The Nelder Mead algorithm has some other stopping criterium, besides the number of evaluations. In the case of an early stop, we restarted the nelder algorithm on a random point in , setting the number of evaluations to the balance of evaluations.  
  
  
  
  
  
    * bobyqa: Bound optimization by quadratic approximation (BOBYQA) <cit.> is a derivative-free optimization algorithm that approximates the response function by quadratic interpolation. We are not aware of any publication that proposed using any of the algorithms in this family to select SVM hyperparameters, but we feel it is an interesting alternative to the Nelder-Mead algorithm.
  
    We used the  function from the  package <cit.>, with the initial point in the center of the space.  We tested bobyqa with N={25,100, 400}. Similarly to nelder, the algorithm will restart at some random point on early stopping. 
  
  
    * sa: Simulated annealing <cit.> starts at a state and probabilistically accepts moving to a different random state given the difference of the evaluations. At “high temperature,” at the beginning of the process, even states that worsen the evaluation are accepted with high probability. As the temperature lowers, the probability of accepting a bad transition also lowers, and at very low temperatures only new states that have better evaluation are accepted. Simulated annealing has been suggested as a method to select SVM hyperparameters <cit.>,.
  
    We used the  function from the  package <cit.>. We used the default values of the simulated annealing for this function and stopped the interaction after N calls, for N={25,100, 400}
  
  
    * pso: Particle swarm optimization <cit.> considers a set of solutions to the optimization problem as a set of particles, which have a momentum (and thus will tend to move in the direction it has been moving before) but are also attracted to other particles that found good solutions to the optimization problem. Particle swarm optimization has been used to select the SVM hyperparameter <cit.>.
  
    We used the  function from the  package <cit.>. We tested with N={25,100, 400}.
  
   
    * cma: Covariance matrix approximation - evolutionary strategy <cit.> is an evolution-inspired algorithm for general non-convex problems. CMA-ES maintains a set of K points in the search space with the corresponding function value. It then adjusts a multidimensional Gaussian distribution on the search space so that points with higher accuracy would have higher likelihood, and generate new points based on this distribution. The algorithm then keeps the best K among the old and the new points and repeats the process.
  
      We used the  package <cit.>. We used the default number of the initial population μ and the default number for the subsequent generations (λ) and the number of generations was computed so the search would not exceed the specified N ={100, 400}.  
  
  
  
    * bogp: Bayesian optimization <cit.> models the response surface as a Gaussian process regression. This modeling allows the optimizer to evaluate not only which data points should be queried because they are likely to result in a good evaluation, but also query which data point will improve its knowledge of the response function itself. Thus a Bayesian optimization model finds a particular balance between exploitation (of the data points it believes will result in good evaluations) and exploration (of new data points that will improve its estimate of the response surface itself).
  
    Bayesian optimization has received a lot of attention recently, but our previous experience with it is that it is computationally costly. We, therefore, opted to use a fast, mostly compiled implementation of Bayesian optimization (as opposed to an R implementation). We used the Python RoBO <cit.> implementation.  We tested with N={100, 400}.  
  
  
  
    * tpe: Tree of Parzen Estimators <cit.> is a Bayesian optimization inspired algorithm, but instead of using Gaussian processes to model the distribution of values of the response surface, it models it as two Gaussian mixture models. TPE is usually much faster than Gaussian process based Bayesian optimization.
  
    We used the HyperOpt python implementation <cit.> of the  TPE, with N={100, 400}.
  
  
  
  
  
  We call grid, ud, normrand, and rand as grid-like search algorithms because the evaluation points in the space are defined before any evaluation need to be performed. We call gridhier and udhier the hierarchical algorithms and we call nelder, bobyqa, sa, cma, pso, bogp, and tpe as general optimization algorithms
  
  
  
  
  
   §.§ SVM specific algorithms
  
  
  
  
  
   
   
    * svmpath:  <cit.> discuss an algorithm that for a given γ, computes all values of C which cause a change in the support vectors set, and thus may cause a change in the accuracy of the classifier. This is called the regularization path of an SVM. Furthermore, computing the regularization path has similar computational cost as computing the SVM itself, and thus, in principle, one could use such an algorithm to search independently for the best hyperparameters. We performed a grid search on the γ value (with N evaluations) for each value of γ we computed the regularization path. The svmpath algorithm not only returns the C values but also the number of points in the support vector for each C. We use this number as an approximation to the error rate of the SVM in the training set.
  
     The svmpath algorithm only uses information from the training set and thus there is no need to use a test set to select the best set of hyperparameters. Therefore, it does not require to use the 5-fold cross-validation - one can use the whole training set and discover the combination of values γ and C that has the lowest number of points in the support vector. Unfortunately, usually, there are many different C values with the same number of support vectors (for a fixed γ), since the svmpath concerns with changes in the support vector and not just its size ( a point may get in as another gets out). In these cases, we selected a random C and γ among the ones with minimal size of the support vector.
  
  
     We used the  function from the  package <cit.> and N ∈{ 5,10,20}.  
  
  
    * asymp:  <cit.> discuss that there should be a good solution to the hyperparameters selection for an RBF SVM in the line logγ = logĈ - log C in the space. Thus they propose that a linear SVM should be used to find the Ĉ constant, and then a 1D grid on the line logγ =logĈ - log C to select the best C and γ for the RBF SVM.
  
    The parametrization of this algorithm involves a pair of numbers: the number of grid points for the C used in the linear SVM problem, and the number of grid points in the logγ =logĈ - log C line to select the pair C and γ.  We used 5+5, 10+10 and 20+20 points.
  
  
    * dbtc:  <cit.> propose that the γ parameter should be chosen using a geometric criterium, in particular, it should be chosen by maximizing the distance between the center of the two classes (in the feature space), which is called distance between two classes (DBTC). DBTC is one of a family of techniques called internal metrics that use some metric on the training set (not on the test set) to select γ such as Xi-alpha bound <cit.>, Generalized approximate cross-validation <cit.>, Maximal Discrepancy <cit.>, among others, DBTC achieves significantly higher accuracy and is significantly fasterthan the others internal metrics  <cit.>.
  
    In our solution, a grid in logγ is searched in order to maximize the DBTC. With the selected γ, a grid on log C is used to select the other hyperparameter. We used 5+5 and 10+10 evaluation points.
  
   We also implemented a sampled version of DBTC, where a sample of 50% of data points on both classes is used to compute the distance. We call this version sdbtc after sampled DBTC and tested it also on 5+5 and 10+10 evaluations. 
   
  
  
    * skl: The popular Python package  <cit.> and the popular SVM solver implementation libSVM <cit.> use as the default for the γ hyperparameter the inverse of the number of dimensions of the data. Sklearn also uses 1.0 as the default value for C.  We do not know of any published evidence for that choice of the default value of C.
  
     One must note the importance of default values for hyperparameters since the less sophisticated user of popular packages usually relies on the default values in the hope that there is good empirical or theoretical evidence for that selection. In fact, the authors are unsophisticated users of many of the packages described above (, , , , and so on) and we did rely on the default values of those packages for the hyperparameters of the corresponding search algorithm. Thus, it is important to test whether the default values for SVM implementations are indeed reasonable choices.
  
  
     We call skl1 the choice of using both default values. And sklN the choice of setting γ as the sklearn default, and performing a grid search on the C, with N={5, 10, 20}.  
  
   
    * sigest:  <cit.> argue that γ could be derived from the distribution of the distances between the data points (in the original space). They argue that any 1/γ in the range between the 10% and the 90% percentile of the distribution of | x_i - x_j |^2 should be an acceptable value. We chose the median (50% percentile) of that distribution as a fixed γ and performed a grid search on log C.  We used the  function of the  package <cit.>. The best γ is computed from the formula and a grid search on log C is performed to select the other hyperparameter. We used N = {5, 10, 20 } for the C grid.
  
    The  function is used in the R package  to determine the default value of the γ parameter in their SVM implementation.
  
    sigest like other algorithms discussed above that use information from the training set to define the value of γ (skl and dbtc) used the whole training set to define the selected γ^*. Using that γ^*, the grid on the C value used the 5-fold cross-validation to select the best value of C.  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Post-search selection procedures
  
  
  
  As discussed, there are reasonable arguments for selecting pairs with lower C and lower γ if more than one pair was found as the best set of hyperparameters.
  
  Let us assume that B = { C_i, γ_i } are the set of best hyperparameters found by a search procedure. But to select the one with the lowest C and the lowest γ, one has to make a decision on which hyperparameter to minimize first. One can select the subset with the lowest C first, and among them, select the one with the lowest γ. Or one may perform the dual algorithm, of selecting first the minimal γ and then C.
  
  It is also reasonable to assume that if the set of B solutions is a plateau on _G, a peak (in _G) would be somewhere in the middle of the plateau, and thus a mean of the C_i and the mean of the γ_i would also be a reasonable solution.
  
  In this research, we will compare the following post-search selection procedures:
  
  
  
    * minCg First selects the minimum value of the C_i and then the minimum value of the γ_i among the subset with the selected minimum C_i. Formally:
    
      Ĉ    = min{ C_i | ⟨ C_i, γ_i ⟩∈ B }
      γ̂    = min{γ_j |  ⟨Ĉ, γ_j ⟩∈ B }
  
  where Ĉ is the selected C and γ̂ the selected
  γ.
  
  
    * mingC: First select the minimum γ_i, and then the minimum C_i among the subset with the selected γ_i.
      
  
    * meanCg The mean of the log C_i and the mean of the
    logγ_i
      
      logĈ    = {log C_i | ⟨ C_i, γ_i
            ⟩∈ B }
      logγ̂    = {logγ_i |  ⟨ C_i,
          γ_i ⟩∈ B }
  
  
  
    * randCg: Select a random pair
    ⟨ C_i, γ_i ⟩ from the set B
  
  
    * maxCg: First select the maximum value of the C_i and then the maximum value of the γ_i among the subset with the selected C_i.
  
    
    * maxgC: First select the maximum γ_i, and then the
    maximum C_i among the subset with the selected γ_i.
    
  
  
  
  The last two procedures maxCg and maxgC are there only to test a procedure that our arguments should point out as worse than the other ones.
  
  
  
  
  
  § EXPERIMENTAL PROTOCOL AND STATISTICAL ANALYSIS
  
  
  
  
  
   §.§ Datasets
  
  
  
  
  
  We used 115 datasets from UCI, first processed by <cit.> and further processed by us. In particular, each multi-class dataset was transformed into a binary problem.
  
  
  For each of the original 115 UCI datasets:
  
  
    * the categorical features were converted to an integer value: 1 for the arbitrary first categorical value, 2 for the second value and so on.
  
    * each feature was scaled to mean equal to 0 and standard deviation equal to 1
  
    * the multi-class datasets were converted to binary by assigning the most frequent class to 0, the second most frequent to 1, the third most frequent to 0, and so on.
  
  
  The distribution of the size of the datasets, the number of features, and the ratio of the less frequent class are shown in the histograms in Figure <ref>.  
  
  
  
  
  
  
  
  
   §.§ Experimental procedure
  
  
  
  
  In general terms, for each dataset, we follow a nested cross-validation procedure, 5-fold inside a 2-fold. The internal cross-validation (5-fold) is used to select the hyperparameters (using the different search procedures) and the external 2-fold is used to estimate the accuracy of the classifier (with the different hyperparameters) on future data (section <ref>).
  
  
  
  
  Formally, each dataset D_i is divided into two subsets S_1i and S_2i with the same size and proportion of classes (this is the external 2-fold). Let us also define 1̂ = 2 and 2̂ = 1, that is, the hat operator denotes the other subset.
  
  Let us define θ(x,i,j) as the value of the hyperparameters that are selected when we use the search algorithm x in subset j of dataset i, that is:
  
      θ(x,i,j) = C, γ∈ x (S_ji, C, γ)
  
  where C, γ∈ x indicates that C and γ are selected from the ones “generated” by algorithm x. (S_ji, C, γ) is defined in Equation <ref> and in this research we use the 5-fold cross-validation to compute it.
  
   In the case of multiple best pairs of hyperparameters, we select a random one from the set, that is, we follow the randCg post-search selection. There are two reasons for this option. First, some of the general search algorithms only return a single solution, so even if the algorithm found many equally good solutions, it seems to us that this unique returned solution is a random one among them. Second, as we will discuss in Section <ref>, there are no significant differences between all of the selection procedures, and thus, any of them can be used to select among equally good pairs of hyperparameters.  
  
  Let us denote (x,i,j) as 
  
      (x,i,j) = (S_ji , θ(x,i,j) )
  
  that is, the highest accuracy found in the response surface of S_ji by algorithm x, again using the 5-fold evaluation . That is, θ(x,i,j) is the pair of best hyperparameters selected by x on subset S_ij and (x,i,j) is the value of _S_ji on that point.
  
  Finally, let us denote β(x,i,j) as: 
  
      β(x,i,j) = (S_
        i | S_ji , θ(x,i,j))
  
  that is, the accuracy on the second subset (S_ i) of an SVM trained on the first subset (S_ji) with hyperparameters selected by algorithm x on the first subset.
  
  For each subset S_ji, we recorded the value of θ(x,i,j) and measured (x,i,j) and β(x,i,j) , that is, for each subset we recorded the position (θ), the height () of the highest point in _S_ji, and the accuracy of testing the SVM trained on S_ji with the best hyperparameters discovered by x on a similar dataset (S_ i).
  
  
  
  
  
  
   §.§ Accuracy gain
  
  
  
  Let us define for each algorithm x and dataset i the mean of the best accuracy for both subsets:
  
      (x,i) = (x,i,j) +  (x,i,) /2
  
  This is a 2-fold estimate of the maximal accuracy achieved by the search algorithm x for dataset i.
  
  
  We will consider the 10 × 10 grid search as the base/standard search algorithm. In fact, a popular guide to SVM <cit.> suggests an 11 × 10 grid search, with C and γ limited by the constraints in Equation <ref>. We want to compare the algorithms described in Section <ref> with the 10 × 10 grid search both in terms of the accuracy achieved and the execution time.
  
  
  We will compute the accuracy gain of an algorithm x on dataset i in relation to the 10 × 10 grid search as:
  
      (x,i) = (x,i) - (,i)
  
  
  
  For each algorithm x, we will compute the average accuracy gain (by averaging over all datasets) and the confidence interval for that mean. The confidence interval is computed using bootstrap with 5000 replicates.
  
  
  
   §.§ Expected accuracy on future data
  
  
  
  
  The expected accuracy of an SVM trained on a dataset D_i with hyperparameters θ(x,i) is:
  
      _g [(g|G,θ(x,i)) ]
  
  
  A 2-fold estimate of that value is: 
  
      E_2(x,i) = (S_ i| S_ji , θ(x, i , j )) + (S_ji| S_
        i, θ(x, i,))/2 = β(x,i,j)+β(x,i,)/2
  
  
  The accuracy gain (in relation to the 10x10 grid) for future data, for dataset D_i and algorithm x is
  
      E_2(x,i) - E_2(,i)
  
  We will compute the accuracy gain on future data for all datasets and algorithms.
  
  
  
  
  
  
  
  
   §.§ Time ratio
  
  
  We measured the total execution time of each search algorithm when executing sequentially – we will not discuss parallelization in this research. We computed the ratio of the execution of the search algorithm x on dataset i, and the execution time of the base search, the 10 × 10 grid search. In all cases, the search algorithm ran sequentially.  We computed the mean and confidence interval on the mean time ratio using bootstrap on 5000 replicates.
  
  The algorithms bogp and tpe were implemented in compiled code with a Python interface. For those two algorithms, we compared the execution time with the grid100 also implemented in Python, using sklearn.
  
  Each combination of algorithms and datasets was allowed to run for at most 5 days; the combinations that exceded that time limit were removed from both the accuracy gain and the time ratio the analysis. 
  
  
    
  
   §.§ Post-search selection procedures
  
  
    For the grid-like algorithms grid, ud, rand, and normrand, we collected all the best solutions for each subset. That is, we collected all solutions to:
  
      θ(x,i,j,k) = C, γ∈ x (S_ji, C, γ)
  
  where k is the index of each best solution.
  
  We then selected a particular θ(x,i,j,k̂) using the different post-search selection algorithms discussed in section <ref>. We then computed the accuracy on the other subset of that choice of hyperparameters when tested on the other subset of the dataset i, that is we computed
  
      (S_ i| S_ji , θ(x, i , j, k̂ )
  
  as an estimate of the accuracy on future data when using all of the selection procedures.
  
  
  
  
  
  
  
   §.§ Statistical analysis
  
  
  
  The usual procedure for comparing different algorithms (in this case, the different search algorithms and the different post-search selection procedures) on different datasets is to follow the one proposed by  <cit.>: a broad comparison of the accuracy of all algorithms using the Friedman test (a non-parametric version of a repeated measure ANOVA), and if its p-value is below the usual 0.05 (for a 95% confidence), a pairwise comparison of all algorithms - the Nemenyi test. Those comparisons with a p-value below 0.05 would be deemed as significantly different (at a 95% confidence).
  
  For the post-search selection comparisons, we perform this procedure. We compared the expected accuracy gain on future data for each of the selection procedures (for algorithms grid400, rand400, ud400, and normrand) to find out which one is the procedure that yields significantly better results for the accuracy on future data.
  
  
  For the accuracy gain and time ratio, we will only compute the 95% confidence interval for those measures.  We do not seek to make a broad claim that search algorithm A is “significantly better” than the other algorithms. One of the reasons is that very likely there will be no such algorithm: we are comparing 47 different combinations of algorithms and parameterizations and it is unlikely that one of such combinations is so much better than the others so that it can be statistically distinguished from all other 46 algorithms with just 115 datasets/measurements. For that to be true the algorithm should be “very good” in finding the maximum of the _G response surface.  That is even less likely true given that the _G is a step function, and thus “around” any very good γ and C points there is a neighborhood of equally good points (with the same value of _G). If the dataset is large, the neighbourhood may be small and maybe none of the other search algorithms found that region, but that is increasingly unlikely for smaller datasets.
  
  Finally, we must also point out that one should not use the confidence intervals computed in this research as a substitute for hypothesis testing.  When comparing two sets of values, confidence intervals can subsume hypothesis testing: if the 95% confidence intervals have no intersection, one can claim that at 95% confidence that the two sets are significantly different (the reverse is not always true <cit.>). But with multiple groups, as pointed out by <cit.> one would have to include the effect of multiple comparisons in the confidence interval calculations. This was not done in this research.
  
  
  
  
  
  
   §.§ Reproducibility
  
  
  
  The datasets used in the paper are available at <https://figshare.com/s/d0b30e4ee58b3c9323fb> as described in <cit.>. The program to run the different procedures and the different classifiers, the results of the multiple runs, and the R program to perform the statistical analysis described in this paper are available at <https://figshare.com/s/35f990ae2b7f65062cca>.
  
  
  
  
  
  
  § RESULTS ON THE COMPARISON OF THE SEARCH ALGORITHMS
  
  
  
  Table <ref> is the main result of the comparisons of the search algorithms. The first column is the name of the search algorithm and the N parameterization. The next three columns are the mean accuracy gain and the low and high limits of the 95% confidence interval for the accuracy gain. The last three columns are the mean time ratio between the algorithm and the base search algorithm, and the low and high limits of the 95% confidence interval. The table is ordered by increasing mean accuracy gain and the small gap in the table indicates the start of the range of algorithms which have a mean accuracy gain of 0 or greater.
  
  Figure <ref> displays the results of Table <ref> in graphical form - which relates the mean accuracy gain and mean log_10 of the time ratio for each combination of algorithm and parametrization.
  
  
  Accuracy gain and time ratio for the search algorithms 
        3c     3cTime ratio
   
  
  algorithm     mean     low     high     mean     low     high 
   
    
  6c
                –
                continued from the previous page 
  
  algorithm     mean     low     high     mean     low     high 
   
    
  6rContinues on next page 
   
  
  
  
  lrrr|rrr
    svmpath5     -0.0388     -0.0483     -0.0312     0.15     0.12     0.20 
   
    asymp10     -0.0384     -0.0502     -0.0304     1.86     1.37     2.47 
   
    asymp20     -0.0346     -0.0458     -0.0268     3.15     2.33     4.16 
   
    skl1     -0.0336     -0.0416     -0.0274     0.02     0.01     0.02 
   
    asymp40     -0.0292     -0.0407     -0.0218     5.79     4.27     7.54 
   
    svmpath10     -0.0292     -0.0376     -0.0224     0.28     0.23     0.35 
   
    svmpath20     -0.0250     -0.0325     -0.0188     0.55     0.45     0.68 
   
    sdbtc5     -0.0206     -0.0262     -0.0161     0.07     0.06     0.08 
   
    skl5     -0.0179     -0.0223     -0.0144     0.06     0.05     0.07 
   
    dbtc5     -0.0176     -0.0227     -0.0139     0.15     0.13     0.17 
   
    sigest5     -0.0144     -0.0174     -0.0119     0.06     0.05     0.06 
   
    sdbtc10     -0.0127     -0.0164     -0.0100     0.11     0.10     0.12 
   
    dbtc10     -0.0111     -0.0145     -0.0086     0.22     0.19     0.24 
   
    sigest10     -0.0102     -0.0129     -0.0081     0.09     0.08     0.10 
   
    sdbtc20     -0.0101     -0.0138     -0.0076     0.19     0.17     0.21 
   
    skl10     -0.0098     -0.0136     -0.0054     0.10     0.08     0.11 
   
    skl20     -0.0094     -0.0138     -0.0057     0.16     0.14     0.18 
   
    sigest20     -0.0086     -0.0113     -0.0065     0.16     0.14     0.17 
   
    sa25     -0.0083     -0.0106     -0.0066     0.34     0.31     0.38 
   
    dbtc20     -0.0074     -0.0101     -0.0053     0.34     0.30     0.37 
   
    grid25     -0.0055     -0.0069     -0.0042     0.27     0.26     0.28 
   
    normrand25     -0.0043     -0.0062     -0.0029     0.26     0.25     0.29 
   
    rand25     -0.0041     -0.0053     -0.0030     0.23     0.22     0.24 
   
    ud25     -0.0039     -0.0052     -0.0026     0.31     0.30     0.32 
   
    cma100     -0.0026     -0.0053     -0.0005     0.98     0.84     1.15 
   
    pso25     -0.0024     -0.0037     -0.0011     0.28     0.26     0.31 
   
    bobyqa25     -0.0020     -0.0039     -0.0006     0.30     0.26     0.35 
   
    sa100     -0.0009     -0.0022     0.0001     1.23     1.15     1.33 
   
    bobyqa400     -0.0006     -0.0027     0.0007     4.27     3.70     4.82 
   
    bobyqa100     -0.0004     -0.0022     0.0008     1.14     1.00     1.30 
   
    nelder25     -0.0002     -0.0015     0.0018     0.33     0.29     0.38 
  
    
  
    gridhier50     0.0000     -0.0011     0.0013     0.54     0.51     0.61 
   
    nelder100     0.0007     -0.0006     0.0020     1.14     0.99     1.29 
   
    rand100     0.0007     -0.0004     0.0020     0.89     0.86     0.92 
   
    cma400     0.0008     -0.0017     0.0029     3.77     3.23     4.32 
   
    ud100     0.0009     -0.0001     0.0021     1.15     1.11     1.20 
   
    normrand100     0.0010     0.0001     0.0021     1.01     0.96     1.07 
   
    udhier50     0.0012     0.0001     0.0024     0.55     0.52     0.59 
   
    nelder400     0.0020     0.0006     0.0037     4.33     3.82     4.85 
   
    pso100     0.0027     0.0018     0.0039     1.09     0.99     1.19 
   
    udhier200     0.0030     0.0020     0.0044     2.09     1.96     2.23 
   
    tpe100     0.0033     0.0024     0.0045     1.12     0.97     1.29 
   
    sa400     0.0034     0.0025     0.0046     4.76     4.42     5.12 
   
    rand400     0.0036     0.0027     0.0048     3.51     3.40     3.60 
   
    ud400     0.0038     0.0029     0.0050     4.48     4.20     4.77 
   
    grid400     0.0038     0.0028     0.0059     3.74     3.65     3.83 
   
    bogp100     0.0039     0.0028     0.0053     36.08     24.43     50.69 
   
    gridhier200     0.0041     0.0033     0.0050     1.88     1.78     2.00 
   
    normrand400     0.0042     0.0032     0.0057     4.00     3.83     4.21 
   
    pso400     0.0055     0.0045     0.0070     4.26     3.83     4.74 
   
    tpe400     0.0057     0.0047     0.0070     3.19     2.82     3.62 
   
    bogp400     0.0059     0.0048     0.0074     461.25     307.12     647.53 
   
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Discussion of some of the  search algorithm
  
  
  
  As discussed, one cannot say that an algorithm “performs significantly worse” than the 10x10 grid using the confidence intervals alone even in these cases where the confidence intervals do not include the 0. Thus, we will only use the mean accuracy gain to make statements on whether an algorithm performs better or worse than the 10x10 grid, and the reader should use both the range of the interval and whether it contains or not the 0, as a measure of the strength of that statement.
  
  
  Two groups of algorithms are clearly distinguished from the others. The SVM specific algorithms asymp and svmpath are clearly inferior to the other ones - they achieve a lower accuracy gain and are much most computationally costly than the other algorithms with a similar number of evaluations. On the other side, the bogp - Bayesian optimization with Gaussian process is very costly, although bogp400 did achieve the highest accuracy gain, it runs 300 times slower than the baseline grid100.
  
  Also, in general, as expected, algorithms with a larger number of probes to the response surface will perform better than algorithms with a lower number of probes and will take correspondingly more computational time. But the more interesting aspect is to consider the continuous line in Figure <ref> as the mean behavior of the search procedures (excluding the svmpath, asymp, and bogp), and consider the ones that are above or below the “expected” behaviour. Among the algorithms with 100 and 400 probes, tpe and pso perform better than the grid-like algorithms, which perform better than the other general optimization procedures (sa, cma, nelder, and bobyqa). Among the grid-like algorithms, the normrand seems slightly better than the others, but we must remind the reader that there it is possible that the normrand results due to a kind of “overfitting” because the algorithm's parameters (the mean and standard deviations of the two distributions) are somewhat tunned to this particular set of datasets.
  
  Among the algorithms with N=25, the derivative-free algorithms nelder and bobyqa seem also interesting choices. Unsurprisingly, the hierarchical algorithms perform above the mean in comparison to the other algorithms.
  
  Finally, among the remaining SVM-specific algorithms, dbtc performs below the average, and the sampled version of dbtc is even worse. sigest seems to perform a little better than the skl algorithm.
  
  
  In general terms, for the very fast search (1, 5 or 10 probes), one should use the sigest algorithm to select the γ and perform a 5 or 10 grid search for C. skl is maybe a simpler (implementation-wise) alternative but there are some accuracy costs. dbtc is not a competitive alternative. For the fast searches (20 to 25 probes) one should use nelder with N=25, or any of the grid-like  with N=25. SVM-specific algorithms are not competitive with these options even at N=20.
  
  For the conventional range on the number of evaluations (N=50 or N=100), one should consider the hierarchical uniform design with N=50 which is faster but performs similarly to the grid-like algorithms. For 100 evaluations one should use either tpe or pso. If one is willing to pay the computational cost of 200 to 400 evaluations, one should again use tpe or pso, a hierarchical grid of 100+100, or the grid-like algorithms at N=400. Finally, among grid-like algorithms, the normrand seems to perform slightly better than the others.
  
  
  
  All algorithms limited to 200 and 400 evaluations performed better than
  the 10x10 grid search (grid100). Clearly, they are, on average, finding points a little
  higher in the response surface; the question is whether it is worth to spend the extra
  effort in selecting the hyperparameters.  We will discuss this question below.
  
  
  
  
   §.§ Is it worth to spend more effort searching for better hyperparameters?
  
  
  
  
  This research, so far, has been discussing search algorithms to find maxima, or peaks, in the _G surface. To remind the readers, _G is a 5-fold estimate of _G which is the true response surface of the SVM classifier when G is used as the learning set and the accuracy is the expectation when the testing data is sampled from the same distribution that G is sampled from.  One really wants to find the best set of hyper-parameters on _G, but one is searching it on _G. The point is that, most likely, _G is a noisy approximation of _G, and peaks in _G may not correspond to peaks in _G. If that is the case, it may not be worth it to spend a lot of computational effort on the most costly search algorithms.
  We will provide two forms of evidence that that is the case.
  
  The first evidence is that we do have another estimate of the _G. As we discussed in Section <ref>, we measure what we called the expected accuracy on future data, which is the accuracy of training the SVM on a subset j of dataset i (with hyperparameters selected on subset j) and testing on the other subset  (Equation <ref>). We will display the average gain on the expected accuracy on future data from each search algorithm in comparison to the 10x10 grid.
  
  The second evidence will compare the peaks on two different _G surfaces. The _G is based on a 5-fold split of the known data G. If a different random generator seed is used, different folds would be created. We call each of an instance of running a search algorithm using a different random seed to generate different 5-folds a run. If the peaks in _G are “grounded” on true peaks of _G, there should not be much of a difference between the best-selected set of hyperparameters between the two runs.
  
  We ran the grid400,ud400,rand400, and normrand400 algorithms on all datasets, with a different instance of the 5-fold cross-validation. The grid400 and the ud400 probe the same points in the space for both _G surfaces, and the random seeds for the rand400 and normrand400 were also set to the same number for the two runs, so also they generate the same sequence of points to be probed in both instances. Therefore, for each of these algorithms, the same set of points were probed, but in two different _G surfaces. We computed different measures of the distance between the two sets of best hyperparameters. We believe that the magnitude of the differences indicate that the maximum in the _G surface for both runs are not “grounded” on a real maximum of the _G surface, but are indeed “noise” particular to each _G surface. 
  
  
  
  
  
    §.§.§ Gains in the expected accuracy on future data
  
  
  
  
  Table <ref> displays the mean accuracy gain on future data in relation to grid100 and the 95% confidence interval on the mean. The algorithms are ordered in increasing order of the accuracy gain.
  
  Notice that all algorithms except for the top 10 (indicated by a small gap in the table) have confidence intervals that include the 0. The fact that a confidence interval includes the 0 is the reverse of the situation we warned the reader in Section <ref>. There, we warned that if the confidence interval did not include the 0, or if two intervals did not have an intersection one could not conclude that the difference was significative because of the corrections needed for the multiple comparisons. Here, the fact that the confidence interval does include the 0 indicates that the differences are not significant – the corrections for multiple comparisons would only enlarge the confidence intervals <cit.>.
  
  Therefore, besides the first 10 SVM-specific algorithms in Table <ref>, none of the other search algorithms seems to be significantly different from grid100 in terms of accuracy on future data. The sa100 algorithm seems to present an anomaly, it also has its confidence interval on the accuracy gain in the negative range, that is, it does not include the 0. But, again, one cannot conclude that the algorithm is significantly worse than grid100. We believe that this result is spurious, in the sense that sa100 is not different than the other algorithms on the expected accuracy gain for future data, it was just an accident that its confidence interval does not include the 0.
  
  
  Accuracy gain on future data  
    algorithm     mean     low     high
  
    
  4c
                –
                continued from the previous page 
  
  algorithm     mean     low     high 
  
   
  3rContinues on next page 
   
  
  
  
  lrrr
  svmpath5     -0.0271     -0.0370     -0.0196 
   
    asymp10     -0.0261     -0.0381     -0.0179 
   
    asymp20     -0.0232     -0.0349     -0.0155 
   
    asymp40     -0.0198     -0.0309     -0.0117 
   
    svmpath10     -0.0192     -0.0273     -0.0123 
   
    skl1     -0.0168     -0.0262     -0.0101 
   
    svmpath20     -0.0158     -0.0231     -0.0098 
   
    sdbtc5     -0.0063     -0.0105     -0.0027 
   
    dbtc5     -0.0056     -0.0094     -0.0025 
   
    sigest5     -0.0041     -0.0079     -0.0006 
  
    
  
    skl5     -0.0032     -0.0078     0.0017 
   
    sa100     -0.0028     -0.0059     -0.0004 
   
    skl10     -0.0025     -0.0071     0.0027 
   
    sa25     -0.0025     -0.0057     0.0007 
   
    dbtc10     -0.0023     -0.0057     0.0007 
   
    sdbtc10     -0.0022     -0.0057     0.0012 
   
    sdbtc20     -0.0018     -0.0052     0.0017 
   
    sigest10     -0.0016     -0.0055     0.0021 
   
    sigest20     -0.0015     -0.0055     0.0019 
   
    cma400     -0.0014     -0.0049     0.0024 
   
    rand100     -0.0013     -0.0038     0.0008 
   
    sa400     -0.0012     -0.0039     0.0010 
   
    cma100     -0.0011     -0.0048     0.0025 
   
    skl20     -0.0011     -0.0062     0.0043 
   
    dbtc20     -0.0010     -0.0041     0.0016 
   
    tpe100     -0.0010     -0.0036     0.0012 
   
    pso25     -0.0009     -0.0031     0.0016 
   
    bogp400     -0.0009     -0.0033     0.0012 
   
    ud400     -0.0007     -0.0032     0.0016 
   
    normrand400     -0.0007     -0.0031     0.0013 
   
    bogp100     -0.0006     -0.0034     0.0018 
   
    tpe400     -0.0006     -0.0031     0.0018 
   
    nelder25     -0.0005     -0.0029     0.0020 
   
    bobyqa400     -0.0004     -0.0027     0.0024 
   
    bobyqa100     -0.0004     -0.0026     0.0018 
   
    gridhier200     -0.0004     -0.0025     0.0014 
   
    grid25     -0.0003     -0.0022     0.0018 
   
    nelder400     -0.0000     -0.0022     0.0025 
   
    nelder100     0.0000     -0.0019     0.0024 
   
    gridhier50     0.0001     -0.0020     0.0021 
   
    udhier50     0.0002     -0.0016     0.0020 
   
    udhier200     0.0002     -0.0020     0.0027 
   
    bobyqa25     0.0004     -0.0018     0.0024 
   
    pso100     0.0004     -0.0018     0.0025 
   
    rand400     0.0004     -0.0022     0.0030 
   
    pso400     0.0005     -0.0015     0.0024 
   
    rand25     0.0005     -0.0020     0.0031 
   
    ud100     0.0006     -0.0020     0.0035 
   
    normrand100     0.0006     -0.0016     0.0028 
   
    ud25     0.0008     -0.0016     0.0028 
   
    grid400     0.0009     -0.0012     0.0031 
   
    normrand25     0.0010     -0.0011     0.0032 
   
  
  
  
  
  
  Figure <ref> displays the gain in accuracy on future data for all search algorithms and their log_10 time ratio in relation to the 10x10 grid. The placement of the search algorithms
  on the horizontal axis is the same as in Figure <ref>, for easier comparison.
  
  
  
  
  
  The important conclusion is whatever gains each algorithm may have in relation to grid100 in finding the peaks in _G surface, those gains seem not to yield better accuracy on future data. 
  
  
  
  
  
    §.§.§ Differences on two _G surfaces
  
  
  
  Table <ref> reports the results on comparing two different _G surfaces. The first line is the number of datasets for which there was only one peak (a single best point C' and γ') in at least one of the runs. The second line is the proportion among those datasets that the two runs resulted in the same best selection. The third line is the average rank of the accuracy of the pair C' and γ' that was selected as the single best solution in one run on the other run. The fourth line is the median of the difference of accuracy between the best selection and the second-best selection in one run. The fifth measurement is the median of the difference in the accuracy of the best-selected solution on each run. Finally, the last measurement is the average distance in space of the two best solutions in both runs.
  
  
  
  
  
  Let us inspect the result for the ud400 search algorithm. 45 datasets resulted in at least one run selecting a single best set of hyperparameters. Let us call C_1 and γ_1 the best selection for the first run and C_2 and γ_2 the best selection for the second run. ^1_1 and ^1_2 are the corresponding accuracies of the best (superscript 1) points for each run (subscript 1 and 2).  ^2_1 and ^2_2 are the accuracies of the second-best (superscript 2) solution for each run. For 9% of the datasets, both runs resulted in the same solution, that is C_1=C_2 and γ_1 = γ_2.  The average rank of the accuracy at the point C_2, γ_2 in the first run and the point C_1, γ_1 in the second run was 20.2. That is, on average, the best point in one run ranks as the 20th best point in the other run. This seems the strongest evidence that the ud400 is mostly searching for peaks in the “noise” of the _G surface - if the peaks in the accuracy were really peaks in the _G(C,γ) surface one would not expect such a large difference from one run to the other. This evidence can also be seen considering the accuracy values. The median difference between ^1_1 and ^2_1 (and ^1_2 and ^2_2), that is, the best and second best solutions of a single run is 8.6 × 10^-4. The median difference between ^1_1 and ^1_2, that is, the accuracy of the best solutions across runs, is ten times as large 8.9 × 10^-3. Again, one would not expect such a difference between the two runs. Finally, the average distance between C_1, γ_1 and C_2, γ_2 is 4.7.
  
  The numbers are similar for the rand400 and normrand400 algorithms, but for them, there were fewer datasets that generated a single best solution. The exception is the grid400 algorithm, which has a much higher rate of agreement of the best solution between the two runs, and therefore a much lower difference between the median average between first and second-best in a single run, the difference of the best accuracy between runs, and the average distance between the best solutions between the two runs. We do not know how to explain this discrepancy of the grid search in relation to the others.
  
  
  
  
  
  
   §.§ Results on the different post search selection procedures
  
  
  
  
  There were 995 combinations of grid-like algorithms (grid, ud, rand, normrand) and subsets (half of each dataset) where there was more than one best solution for the hyperparameters. For these combinations, the median number of solutions was 5, the mean number of solutions was 24.57, the minimum 2, and the maximum 293.
  
  Table <ref> displays the mean rank of the accuracy on future data for each selection procedure. As one can see, there seems to be very little difference between the ranks. In fact, the Friedman test on the accuracy on future data of the selection procedures results in a p-value of 0.42 (chi-squared = 4.948, df = 5) which shows that there are no significant differences among the selection procedures.
  
  
  
  
  
    It could be the case that the non-significant differences are due to a large number of comparisons. Table <ref> only compares the “more reasonable” selection procedures: minCg mingC, randomCg, and meanCg. The resulting Friedman test has a p-value of 0.80 ( chi-squared = 1.0188, df = 3). Therefore, there is no evidence of differences among the post-search selection procedures.
    
  
  
  
  
  Finally, even if we restrict the number of solutions to 6 or above (which results in 464 combinations of grid-like algorithms and subsets) the p-value of the Friedman test for all 6 selection procedures is 0.63 (chi-squared = 3.4738, df = 5) and for the 4 selection procedures 0.56 (chi-squared = 2.0767, df = 3).
  
  So, there is ample evidence that despite reasonable arguments to use one or another post-search selection procedure, there are no significant differences in future data accuracy among them.
  
  
  
  
  
   §.§ Hyperparameters selected by the best algorithms
  
  
  
  Figure <ref> plots the best log C and logγ found by all the algorithms with 200 and 400 probes. The mean log C for the best algorithms without the normrand400 is 5.37 and the standard deviation is 5.33. The mean of the logγ is -5.56 and the standard deviation is 3.57.  These numbers are the source for our suggestion of the normrand algorithm.  A mean of 5 and a standard deviation of 5 for the log C and a mean -5 and a standard deviation of 5 for the logγ are mnemonic approximations of those values.
  
  
  
  
  
  
  
  
  
  
  § DISCUSSION
  
  
  
  
  
  
  
  
  An important issue to discuss is whether the results of this research can be generalized to any future dataset. The data sets used in this research are a large sample of those available at UCI in 2014, collected by the authors of <cit.>. We believe they can be considered as an unbiased sample of real-world data sets. But some classes of data sets are not represented in this set: very large data sets and data sets with a high number of dimensions. Fortunately, RBF SVM are usually not applied in these types of data sets: SVM is at least quadratic on the number of training examples and therefore it does not scale well to large data sets, with the possible exception of online SVM solvers; SVM with RBF kernel is usually not used on large dimensionality problems since for these cases the linear SVM is probably a better solution. Finally, all datasets are binary classification problems. We do not believe that for multi-class problems, where the SVMs are used in one-vs-one or one-vs-all configurations the conclusions would be different.
  
  
  
  
  
  The choice of search algorithms was necessarily ad-hoc. We believe that this research covered most of the important proposal for SVM hyperparameter search algorithms. A somewhat non-conventional inclusion was the uniform design (ud) search – we wanted to include at least one algorithm from the research area of Design of Computer Experiments, which has not been very well represented in the hyperparameter search research.  The most notable absentee in this research are the many versions of sequential testing <cit.>.
  
  Also, we must point out that we made no attempt to select hyperparameters for the search procedures themselves. The more complex general optimization algorithms such as pso, sa, cma have themselves hyperparameters: one has to define a function that will reduce the temperature in sa; one has to define the number of particles and the relative strength of the momentum and the attractive force among the particles in pso; and so on. Although we cited research that proposed using general optimization algorithms such as pso, sa, cma and so on, <cit.> we did not attempt to reproduce the algorithms hyperparameters as proposed by the authors or perform any search for those hyperparameters - we relied on the default of the implementations used. As we will discuss below, we do not think that a better choice of search hyperparameters would have resulted in a better selection of the SVM hyperparameters. Similarly, we do not believe that the inclusion of some of the missing search algorithms such as ant colony optimisation would have changed the main conclusions of this research in any significant way.
  
  
  
  
  
  The time measurement results are less generalizable: the results reflect the 2019 implementation in R of the search algorithms, and the relative speeds of the different search algorithms may change as better and faster implementations become available. Also, we must point out that the timing data refer to the sequential execution of the search algorithms. Parallel execution is trivial for grid-like algorithms and there is published research on parallel versions of pso, sa, and bogp.
  
  
  We are not aware of any research published that systematically compares different search algorithms for SVM, with the exception of <cit.> who compares on 70 data sets 6 different search algorithms: grid search, random search, default values for the hyperparameters, genetic algorithms (GA), estimation of distribution algorithms (EDA), and particle swarm optimization (pso). Different to this research, they run experiments of up to 10000 evaluations, but they use a non-standard cross-validation to tune and measure the accuracy of the SVM. What they call the S-CV (single cross-validation step) divides the data set into K subsets, where they use a 2-fold to train with the different combinations of hyperparameters, which are tested on one of the folds remaining. The last fold is only used once to measure the accuracy of the best selection of hyperparameters. The process is repeated K times, changing the last two folds and using the remained in training. We used the more traditional nested cross-validation - the subsets are an external 2-fold CV (to compute the accuracy) and the selection of the hyperparameters is performed by the internal 5-fold CV (one selection per fold). We do not know if the simpler, less costly S-CV, introduces bias to the evaluation, but we do know that nested CV is considered a more standard way of evaluating many classifiers <cit.>.  They performed the standard significance tests on the 6 search algorithms and did not find any significant differences except for the default settings which performed worse than the other 5 algorithms.
  
  There has been some recent research on hyperparameter optimization in general <cit.> but they tend to focus on the more challenging problem of searching high dimensional spaces <cit.>.
  
  
  
  
  
  Finally, this research has some strong similarities to <cit.>. There the authors fixed the search procedure to select the SVM hyperparameters (an 11x10 grid search) and varied the resampling procedure used to compute the _G. That paper finds a similar conclusion that it is not worth it to use very precise and costly resampling procedures for the selection of hyperparameters. In particular, the authors suggest using a 2 or 3-fold resampling to select the hyperparameters. Here, we find that it is not worth it to use costly search procedures. In both cases, it seems that efforts in finding precisely the maxima in _G or that measuring precisely _G are wasteful.  The _G surface does not seem to have narrow, high maxima - narrow so that it is necessary to not only measure it precisely but also to spend a lot of search effort to find it, and high because the difference in accuracy between this high peak and the other not so good points is significant enough.  _G for SVMs seem to be smooth (at least near the maxima) and one should not expect that either measuring it more precisely or probing it many times will yield a “sweet spot” in accuracy. We can only speculate whether this might be true for other classification algorithms.
  
  
  
  § CONCLUSIONS
  
  
  
  
  
  Regarding the main goal of this research, the comparison of different search algorithms for the C and γ hyperparameters of an RBF SVM, we found that:
  
  
    * more evaluations of the cross-validated accuracy will, in general, result in the selection of hyperparameters that have higher cross-validated accuracy, which is not surprising.
  
  
    * if one is willing to allow many evaluations (100 or 400),  tree of Parzen estimators  <cit.> and particle swarm optimization <cit.> seem to be the two algorithms that achieve a good balance of execution time and finding higher accuracy points in the _G surface.
  
    * Bayesian optimization <cit.> does find high accuracy points but at a very high computational cost.
  
  
    * all grid-like algorithms grid, rand, ud, and normrand seem to have similar performance, but the accuracy gains for normrand seems a little higher than the other three.
  
  
    * the hierarchical versions of grid and ud are, not surprisingly, good alternatives to the flat versions of the algorithms.
  
  
    * at lower evaluation points (25) nelder and pso seem to be the best alternative
  
  
    * at even lower number of evaluations (10 and 20) the default value of the γ parameter selected by sigest  <cit.> seems a better solution than the alternatives proposed by Sklearn and distance between two classes. The sigest default value for 1/γ tested in this research is the median distance squared among a 50% sample of the data points.
  
    
    * the SVM-specific algorithms svmpath <cit.> and asym <cit.> are not competitive alternatives.
    
    
  
   Regarding the second goal of this research, efforts in finding better hyperparameters with higher values of the cross-validated accuracy seems not to be correlated to better expected accuracy on future data. There is no important difference in (an estimative of) the expected accuracy on future data among the best performing algorithm and the grid100. 
  
  Regarding the third goal of this research, we could not find any significant difference between the different post-search selection procedures. That is, despite arguments that one should select a lower C and lower γ if multiple best solutions were found, there is no evidence that this selection procedure is any better than a random selection, or other different procedures in terms of accuracy on future data. 
  
  Regarding the fourth goal of this research, at a first approximation, the distribution of the best log C (for datasets whose features are scaled to zero mean and unit standard deviation) is a normal distribution of mean equal to 5 and standard deviation equal to 5. And  the best logγ follows a  normal distribution with mean equal to -5 and standard deviation of 5.
  
  Finally, 
  an important conclusion of this paper is that is in likely not worth it to spend too much computational time searching for better hyperparameters. It is likely that the search algorithms are finding “noisy” peaks in the _G surface which are not representative of “real” peaks in the _G surface. If that is indeed the case, then one should not spend a lot of computational effort searching for the best set of hyperparameters, and our conclusions that pso and tse are, in general, the best search algorithms for a high number of evaluations are only of theoretical interest. In this case, one should choose algorithms with a lower number of evaluations. Given that grid like algorithms can be trivially parallelized, if the user has enough computational instances, we believe that a grid-like search with 25 or a few more evaluations is likely enough to find good hyperparameters for most of the problems. 
  
  
  
  
  
  
  54
  
  
  
  urlstyle
  
  
  [Aarts and Korst(1988)]aarts1988simulated
  Aarts E, Korst J (1988) Simulated annealing and Boltzmann machines. John
    Wiley and Sons Inc.
  
  [Anguita et al(2003)Anguita, Ridella, Rivieccio, and
    Zunino]Anguita2003109
  Anguita D, Ridella S, Rivieccio F, Zunino R (2003) Hyperparameter design
    criteria for support vector classifiers. Neurocomputing 55(1–2):109 – 134
  
  [Bendtsen(2012)]pso
  Bendtsen C (2012) pso: Particle swarm optimization. version 1.0.3 
  
  [Bergstra and Bengio(2012)]bergstra2012random
  Bergstra J, Bengio Y (2012) Random search for hyper-parameter optimization.
    Journal of Machine Learning Research 13:281–305
  
  [Bergstra et al(2011)Bergstra, Bardenet, Bengio, and
    Kégl]bergstra2011algorithms
  Bergstra J, Bardenet R, Bengio Y, Kégl B (2011) Algorithms for
    hyper-parameter optimization. In: Advances in Neural Information Processing
    Systems, pp 2546–2554
  
  [Bergstra et al(2013)Bergstra, Yamins, and Cox]bergstra2013making
  Bergstra J, Yamins D, Cox DD (2013) Making a science of model search:
    Hyperparameter optimization in hundreds of dimensions for vision
    architectures. Journal of Machine Learning Research 28:115–123
  
  [Caputo et al(2002)Caputo, Sim, Furesjo, and
    Smola]caputo2002appearance
  Caputo B, Sim K, Furesjo F, Smola A (2002) Appearance-based object recognition
    using SVMs: which kernel should I use? In: NIPS Workshop on statistical
    methods for computational experiments in visual processing and computer
    vision.
  
  [Cawley and Talbot(2007)]nelder-lssvm
  Cawley G, Talbot N (2007) Preventing over-fitting during model selection via
    bayesian regularisation of the hyper-parameters. Journal of Machine Learning
    Research 8:841–861
  
  [Cawley and Talbot(2010)]cawley2010over
  Cawley GC, Talbot NL (2010) On over-fitting in model selection and subsequent
    selection bias in performance evaluation. Journal of Machine Learning
    Research 11:2079–2107
  
  [Chang and Lin(2011)]libsvm
  Chang CC, Lin CJ (2011) LIBSVM: a library for support vector machines. ACM
    Transactions on Intelligent Systems and Technology 2(3):27
  
  [Cohen et al(2005)Cohen, Ruch, and Hilario]nelder1
  Cohen G, Ruch P, Hilario M (2005) Model selection for support vector
    classifiers via direct simplex search. In: FLAIRS Conference, pp 431–435
  
  [Demšsar(2006)]demsar
  Demšsar J (2006) Statistical comparisons of classifiers over multiple
    data sets. Journal of Machine Learning Research 7:1–30
  
  [Duarte and Wainer(2017)]duarte2017empirical
  Duarte E, Wainer J (2017) Empirical comparison of cross-validation and internal
    metrics for tuning SVM hyperparameters. Pattern Recognition Letters 88:6–11
  
  [Eberhart et al(1995)Eberhart, Kennedy et al]eberhart1995new
  Eberhart R, Kennedy J, et al (1995) A new optimizer using particle swarm
    theory. In: International Symposium on Micro Machine and Human Science, vol 1, pp 39–43
  
  [Fang et al(2000)Fang, Lin, Winker, and Zhang]ud1
  Fang KT, Lin D, Winker P, Zhang Y (2000) Uniform design: theory and
    application. Technometrics 42(3):237–248
  
  [Fernández-Delgado et al(2014)Fernández-Delgado, Cernadas,
    Barro, and Amorim]delgado14
  Fernández-Delgado M, Cernadas E, Barro S, Amorim D (2014) Do we need
    hundreds of classifiers to solve real world classification problems? Journal
    of Machine Learning Research 15:3133–3181
  
  [Hansen et al(2003)Hansen, Müller, and
    Koumoutsakos]hansen2003reducing
  Hansen N, Müller SD, Koumoutsakos P (2003) Reducing the time complexity of
    the derandomized evolution strategy with covariance matrix adaptation
    (cma-es). Evolutionary Computation 11(1):1–18
  
  [Hastie(2012)]svmpath
  Hastie T (2012) svmpath: the SVM Path algorithm. version 0.953 
  
  [Hastie et al(2004)Hastie, Rosset, Tibshirani, and Zhu]path1
  Hastie T, Rosset S, Tibshirani R, Zhu J (2004) The entire regularization path
    for the support vector machine. Journal of Machine Learning Research
    5:1391–1415
  
  [Hsu et al(2003)Hsu, Chang, and Lin]libsvmguide
  Hsu CW, Chang CC, Lin CJ (2003) A practical guide to support vector
    classification. Tech. rep., Department of Computer Science National Taiwan
    University, updated: May 19, 2016
  
  [Huang and Dun(2008)]pso2
  Huang CL, Dun JF (2008) A distributed PSO–SVM hybrid system with feature
    selection and parameter optimization. Applied Soft Computing 8(4):1381–1391
  
  [Huang et al(2007)Huang, Lee, Lin, and Huang]ud-svm
  Huang CM, Lee YJ, Lin D, Huang SY (2007) Model selection for support vector
    machines via uniform design. Computational Statistics & Data Analysis
    52(1):335–346
  
  [Hutter et al(2011)Hutter, Hoos, and
    Leyton-Brown]hutter2011sequential
  Hutter F, Hoos HH, Leyton-Brown K (2011) Sequential model-based optimization
    for general algorithm configuration. In: International Conference on Learning
    and Intelligent Optimization, Springer, pp 507–523
  
  [Hutter et al(2019)Hutter, Kotthoff, and Vanschoren]automl
  Hutter F, Kotthoff L, Vanschoren J (2019) Automated Machine Learning. Springer
  
  [hyp(accessed 2019)]hyperopt
   Hyperopt: Distributed asynchronous hyper-parameter
    optimization. <http://github.com/hyperopt/hyperopt>
  
  
  [Imbault and Lebart(2004)]sa3
  Imbault F, Lebart K (2004) A stochastic optimization approach for parameter
    tuning of support vector machines. In: International Conference on Pattern
    Recognition, IEEE, vol 4, pp 597–600
  
  [Joachims(2000)]joachims2000maximum
  Joachims T (2000) The maximum-margin approach to learning text classifiers:
    Methods, theory, and algorithms. PhD thesis, Doctoral dissertation, Dortmund
    University
  
  [Karatzoglou et al(2004)Karatzoglou, Smola, Hornik, and
    Zeileis]kernlab
  Karatzoglou A, Smola A, Hornik K, Zeileis A (2004) kernlab – an S4 package
    for kernel methods in R. Journal of Statistical Software 11(9):1–20
  
  [Keerthi(2002)]keerthi2002efficient
  Keerthi S (2002) Efficient tuning of SVM hyperparameters using radius/margin
    bound and iterative algorithms. IEEE Transactions on Neural Networks
    13(5):1225–1229
  
  [Keerthi and Lin(2003)]asymptotic
  Keerthi S, Lin CJ (2003) Asymptotic behaviors of support vector machines with
    Gaussian kernel. Neural computation 15(7):1667–1689
  
  [Keerthi et al(2006)Keerthi, Sindhwani, and
    Chapelle]keerthi2006efficient
  Keerthi S, Sindhwani V, Chapelle O (2006) An efficient method for
    gradient-based adaptation of hyperparameters in SVM models. In: Advances in
    neural information processing systems, pp 673–680
  
  [Klein et al(2017)Klein, Falkner, Mansur, and
    Hutter]klein-bayesopt17
  Klein A, Falkner S, Mansur N, Hutter F (2017) Robo: A flexible and robust
    bayesian optimization framework in python. In: NIPS 2017 Bayesian
    Optimization Workshop
  
  [Krueger et al(2015)Krueger, Panknin, and Braun]krueger2015fast
  Krueger T, Panknin D, Braun M (2015) Fast cross-validation via sequential
    testing. Journal of Machine Learning Research 16:1103–1155
  
  [Lin et al(2008a)Lin, Lee, Chen, and Tseng]sa1
  Lin SW, Lee ZJ, Chen SC, Tseng TY (2008a) Parameter determination
    of support vector machine and feature selection using simulated annealing
    approach. Applied soft computing 8(4):1505–1512
  
  [Lin et al(2008b)Lin, Ying, Chen, and Lee]pso1
  Lin SW, Ying KC, Chen SC, Lee ZJ (2008b) Particle swarm
    optimization for parameter determination and feature selection of support
    vector machines. Expert Systems with Applications 35(4):1817–1824
  
  [Ludbrook(2000)]ludbrook2000multiple
  Ludbrook J (2000) Multiple inferences using confidence intervals. Clinical and
    Experimental Pharmacology and Physiology 27(3):212–215
  
  [Mantovani et al(2015)Mantovani, Rossi, Vanschoren, Bischl, and
    de Carvalho]mantovani2015effectiveness
  Mantovani RG, Rossi AL, Vanschoren J, Bischl B, de Carvalho AC (2015)
    Effectiveness of random search in SVM hyper-parameter tuning. In:
    International Joint Conference on Neural Networks (IJCNN), pp 1–8
  
  [Nelder and Mead(1965)]nelder1965simplex
  Nelder J, Mead R (1965) A simplex method for function minimization. The
    Computer Journal 7(4):308–313
  
  [Pai and Hong(2005)]sa2
  Pai PF, Hong WC (2005) Support vector machines with simulated annealing
    algorithms in electricity load forecasting. Energy Conversion and Management
    46(17):2669–2688
  
  [Pedregosa et al(2011)Pedregosa, Varoquaux, Gramfort, Michel, Thirion,
    Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos,
    Cournapeau, Brucher, Perrot, and Duchesnay]scikit-learn
  Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M,
    Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D,
    Brucher M, Perrot M, Duchesnay E (2011) Scikit-learn: Machine learning in
    Python. Journal of Machine Learning Research 12:2825–2830
  
  [Pelikan et al(1999)Pelikan, Goldberg, and Cantú-Paz]bo1
  Pelikan M, Goldberg DE, Cantú-Paz E (1999) BOA: The Bayesian
    optimization algorithm. In: Annual Conference on Genetic and Evolutionary
    Computation, pp 525–532
  
  [Powell(2009)]powell2009bobyqa
  Powell M (2009) The BOBYQA algorithm for bound constrained optimization
    without derivatives. Tech. Rep. Report NA2009/06, Department of Applied
    Mathematics and Theoretical Physics, University of Cambridge
  
  [Schenker and Gentleman(2001)]schenker2001judging
  Schenker N, Gentleman J (2001) On judging the significance of differences by
    examining the overlap between confidence intervals. The American Statistician
    55(3):182–186
  
  [Sklearn(2019)]sklsvm
  Sklearn (2019) RBF SVM parameters.
    <https://scikit-learn.org/stable/auto_examples/svm/plot_rbf_parameters.html>
  
  [Snoek et al(2012)Snoek, Larochelle, and Adams]bo2
  Snoek J, Larochelle H, Adams RP (2012) Practical bayesian optimization of
    machine learning algorithms. In: Advances in Neural Information Processing
    Systems, pp 2951–2959A
  
  [Sun et al(2010)Sun, Zheng, Li, and Zhou]sun2010analysis
  Sun J, Zheng C, Li X, Zhou Y (2010) Analysis of the distance between two
    classes for tuning SVM hyperparameters. IEEE Transactions on Neural
    Networks 21(2):305–318
  
  [Trautmann et al(2011)Trautmann, Mersmann, and Arnu]cmaes
  Trautmann H, Mersmann O, Arnu D (2011) cmaes: Covariance matrix adapting
    evolutionary strategy. version1.0-11 
  
  [Wahba et al(2001)Wahba, Lin, Lee, and Zhang]Wahba01onthe
  Wahba G, Lin Y, Lee Y, Zhang H (2001) On the relation between the gacv and
    Joachims' xi-alpha method for tuning support vector machines, with extensions
    to the non-standard case. Tech. Rep. 1039, Statistics Department University
    of Wisconsin, Madison WI
  
  [Wainer(2016)]arxiv
  Wainer J (2016) Comparison of 14 different families of classification
    algorithms on 115 binary datasets. ArXiv e-prints 1606.00930
  
  [Wainer and Cawley(2017)]wainer2017empirical
  Wainer J, Cawley G (2017) Empirical evaluation of resampling procedures for
    optimising svm hyperparameters. Journal of Machine Learning Research
    18(15):1–35
  
  [Wessing(2018)]wessing2018proper
  Wessing S (2018) Proper initialization is crucial for the Nelder–Mead simplex
    search. Optimization Letters pp 1–10
  
  [Y. Xiang et al(2013)Y. Xiang, Gubian, Suomela, and Hoeng]GenSA
  Y Xiang, Gubian S, Suomela B, Hoeng J (2013) Generalized simulated annealing
    for efficient global optimization: the GenSA package for R. The R Journal
    Volume 5/1, June 2013
  
  [Ypma et al(2014)]nlopt
  Ypma J, et al (2014) nloptr: R interface to NLopt. version 1.0.4 
  
  [Zhu(2016)]mixtox
  Zhu X (2016) mixtox: Curve fitting and mixture toxicity assessment. version 1.3
    
  
  
  
  
  
  
  