Source code of the Bayesian SVM described in the paper by Wenzel et al. "Bayesian Nonlinear Support Vector Machines for Big Data"
We propose a fast inference method for Bayesian nonlinear support vector machines that leverages stochastic variational inference and inducing points. Our experiments show that the proposed method is faster than competing Bayesian approaches and scales easily to millions of data points. It provides additional features over frequentist competitors such as accurate predictive uncertainty estimates and automatic hyperparameter search.READ FULL TEXT VIEW PDF
Support vector machines have attracted much attention in theoretical and...
We propose an algorithm for exploring the entire regularization path of
A mean field variational Bayes approach to support vector machines (SVMs...
Tackling pattern recognition problems in areas such as computer vision,
We introduce a new Bayesian multi-class support vector machine by formul...
In this paper, we consider modelling interaction between a set of variab...
Support Vector Machines (SVMs) are one of the most popular supervised
Source code of the Bayesian SVM described in the paper by Wenzel et al. "Bayesian Nonlinear Support Vector Machines for Big Data"
Statistical machine learning branches into two classic strands of research: Bayesian and frequentist. In the classic supervised learning setting, both paradigms aim to find, based on training data, a functionthat predicts well on yet unseen test data. The difference in the Bayesian and frequentist approach lies in the treatment of the parameter vector of this function. In the frequentist setting, we select the parameter that minimizes a certain loss given the training data, from a restricted set of limited complexity. In the Bayesian
school of thinking, we express our prior belief about the parameter, in the form of a probability distribution over the parameter vector. When we observe data, we adapt our belief, resulting in a posterior distribution over
Advantages of the Bayesian approach include automatic treatment of hyperparameters and direct
quantification of the uncertainty111 Note that frequentist approaches can also lead to other forms of uncertainty estimates, e.g. in form of confidence intervals. But since the classic SVM does not exhibit a probabilistic formulation these uncertainty estimates cannot be directly computed.
Note that frequentist approaches can also lead to other forms of uncertainty estimates, e.g. in form of confidence intervals. But since the classic SVM does not exhibit a probabilistic formulation these uncertainty estimates cannot be directly computed.of the prediction in the form of class membership probabilities which can be of tremendous importance in practice. As examples consider the following. (1) We have collected blood samples of cancer patients and controls. The aim is to screen individuals that have increased likelihood of developing cancer. The knowledge of the uncertainty in those predictions is invaluable to clinicians. (2) In the domain of physics it is important to have a sense about the certainty level of predictions since it is mandatory to assert the statistical confidence in any physical variable measurement. (3) In the general context of decision making, it is crucial that the uncertainty of the estimated outcome of an action can be reliably determined.
Recently, it was shown that the support vector machine (SVM) —which is a classic supervised classification algorithm— admits a Bayesian interpretation through the technique of data augmentation [2, 3]. This so-called Bayesian nonlinear SVM
combines the best of both worlds: it inherits the geometric interpretation, its robustness against outliers, state-of-the-art accuracy, and theoretical error guarantees  from the frequentist formulation of the SVM, but like Bayesian methods it also allows for flexible feature modeling, automatic hyperparameter tuning, and predictive uncertainty quantification.
However, existing inference methods for the Bayesian support vector machine (such as the expectation conditional maximization method introduced in ) scale rather poorly with the number of samples and are limited in application to datasets with thousands of data points . Based on stochastic variational inference  and inducing points , we develop in this paper a fast and scalable inference method for the nonlinear Bayesian SVM.
Our experiments show superior performance of our method over competing methods for uncertainty quantification of SVMs such as Platt’s method . Furthermore, we show that our approach is faster (by one to three orders of magnitude) than the following competitors: expectation conditional maximization (ECM) for nonlinear Bayesian SVM by , Gaussian process classification , and the recently proposed scalable variational Gaussian process classification method . We apply our method to the domain of particle physics, namely on the SUSY dataset  (a standard benchmark in particle physics containing 5 million data points) where our method takes only 10 minutes to train on a single CPU machine.
Our experiments demonstrate that Bayesian inference techniques are mature enough to compete with corresponding frequentist approaches (such as nonlinear SVMs) in terms of scalability to big data, yet they offer additional benefits such as uncertainty estimation and automated hyperparameter search.
Our paper is structured as follows. In section 2 we discuss related work and review the Bayesian nonlinear SVM model in section 3. In section 4 we propose our novel scalable inference algorithm, show how to optimize hyperparameters and obtain an approximate predictive distribution. We discuss also the special case of the linear SVM, for which we propose a specially tailored fast inference algorithm. Section 5 concludes with experimental results.
There has recently been significant interest in utilizing max-margin based discriminative Bayesian models for various applications. For example,  employs a max-margin based Bayesian classification to discover latent semantic structures for topic models,  uses a max-margin approach for efficient Bayesian matrix factorization, and 
develops a new max-margin approach to Hidden Markov models.
All these approaches apply the Bayesian reformulation of the classic SVM introduced by . This model is extended by  to the nonlinear case. The authors show improved accuracy compared to standard methods such as (non-Bayesian) SVMs and Gaussian process (GP) classification.
However, the inference methods proposed in  and  have the drawback that they partially rely on point estimates of the latent variables and do not scale well to large datasets. In  the authors apply mean field variational inference to the linear case of the model, but their proposed technique does not lead to substantial performance improvements and neglects the nonlinear model.
Uncertainty estimation for SVMs is usually done via Platt’s technique 
, which consists of applying a logistic regression on the function scores produced by the SVM. In contrast, our technique directly yields a sound predictive distribution instead of using a heuristically motivated transformation. We make use of the idea of inducing point GPs to develop a scalable inference method for the Bayesian nonlinear SVM. Sparse GPs using pseudo-inputs were first introduced in. Building on this idea Hensman et al. developed a stochastic variational inference scheme for GP regression and GP classification [7, 10]. We further extend this ideas to the setting of Bayesian nonlinear SVM.
Let be observations where is a feature vector with corresponding labels . The SVM aims to find an optimal score function by solving the following regularized risk minimization objective:
where is a regularizer function controlling the complexity of the decision function , and is a hyperparameter to adjust the trade-off between training error and the complexity of . The loss
is called hinge loss. The classifier is then defined as.
For the case of a linear decision function, i.e. , the SVM optimization problem (1) is equivalent to estimating the mode of a pseudo-posterior
Here denotes a prior such that . In the following we use the prior , where is a positive definite matrix. From a frequentist SVM view, this choice generalizes the usual -regularization to non-isotropic regularizers. Note that our proposed framework can be easily extended to other regularization techniques by adjusting the prior on (e.g. block -norm regularization which is known as multiple kernel learning ). In order to obtain a Bayesian interpretation of the SVM, we need to define a pseudolikelihood such that the following holds,
By introducing latent variables (data augmentation) and making use of integral identities stemming from function theory,  show that the specification of in terms of the following marginal distribution satisfies (2):
Writing for the matrix of data points and , the full conditional distributions of this model are
with , , and where
denotes a generalized inverse Gaussian distribution. Thelatent variables
of the model scale the variance of the full posteriors locally. The model thus constitutes a special case of a normal variance-mean mixture, where we implicitly impose the improper prioron . This could be generalized by using a generalized inverse Gaussian prior on , leading to a conjugate model for . Henao et al. show that in the case of an exponential prior on
, this leads to a skewed Laplace full conditional for. Note that this, however, destroys the equivalency to the frequentist linear SVM.
By using the ideas of Gaussian processes , Henao et al. develop a nonlinear (kernelized) version of this model . They assume a continuous decision function to be drawn from a zero-mean Gaussian process , where is a kernel function. The random Gaussian vector corresponds to evaluated at the data points. They substitute the linear function by in (3) and obtain the conditional posteriors
with . For a test point the conditional predictive distribution for under this model is
where , , . The conditional class membership probability is
where is the probit link function.
Note that the conditional posteriors as well as the class membership probability still depend on the local latent variables . We are interested in the marginal predictive distributions, but unfortunately the latent variables cannot be integrated out analytically. Both  and  propose MCMC-algorithms and stepwise inference schemes similar to EM-algorithms to overcome this problem. These methods do not scale well to big data problems and the probability estimation still relies on point estimates of the -dimensional . We overcome these problems proposing a scalable inference method and obtaining approximate marginal predictive distributions (that are not conditioned on ).
In the following we develop a fast and reliable inference method for the Bayesian nonlinear SVM. Our method builds on the idea of using inducing points for Gaussian Processes in a stochastic variational inference setting  that scales easily to millions of data points. We proceed by first discussing a standard batch variational scheme in section 4.1 and then in section 4.2 we develop our fast and scalable inference method. We show how to automatically tune hyperparameters in section 4.3 and obtain uncertainty estimates for predictions in section 4.4. Finally, we discuss the special case of the Bayesian linear SVM in section 4.5.
The idea of variational inference is to approximate the typically intractable posterior of a probabilistic model by a variational (typically factorized) distribution. We find the optimal approximating distribution by maximizing a lower bound on the evidence (the so-called ELBO) with respect to the parameters of the variational distribution, which is equivalent to minimizing the Kullback-Leibler divergence between the variational distribution and the posterior[18, 19].
In this section we first develop a batch variational inference scheme [18, 19], which uses the full dataset in every iteration. We follow the structured mean field approach and choose the variational distributions within the same families as the full conditional distributions , with . The coordinate ascent updates can be computed by the expected natural parameters of the corresponding full conditionals (5) leading to
This concludes the batch variational inference scheme.
The downside of this approach is that it does not scale to big datasets. The covariance matrix of the variational distribution has dimension and has to be updated and inverted at every inference step. This operation exhibits the computational complexity , where
is the number of data points. Furthermore, in this setup we cannot apply stochastic gradient descent. We show how to overcome both problems in the next section paving the way to perform inference on big datasets.
We aim to develop a stochastic variational inference (SVI) scheme using only minibatches of the data in each iteration. The Bayesian nonlinear SVM model does not exhibit a set of global variables. Both the number of latent variables and the observations of the latent GP grow with number of data points (c.f. eq.5), i.e. they are local variables. This hinders us from directly developing a SVI scheme. We make use of the concept of inducing points  imposing a sparse GP acting as global variable. This allows us to apply SVI and reduces the complexity to , where is the number of inducing points, which is independent of the number of data points.
We augment our original model (5) with inducing points. Let be pseudo observations at inducing locations . We employ a prior on the inducing points, and connect and setting
where is the kernel matrix resulting from evaluating the kernel function between all inducing points locations, is the cross-covariance between the data points and the inducing points and is given by
. The augmented model exhibits the joint distribution
Note that we can recover the original joint distribution by marginalizing over . We now aim to apply the methodology of variational inference to the marginal joint distribution . We impose a variational distribution on the inducing points . We follow  and apply Jensen’s inequality to obtain a lower bond on the following intractable conditional probability,
Plugging the lower bound into the standard evidence lower bound (ELBO)  leads to the new variational objective
The expectations can be computed analytically (details are given in the appendix) and we obtain in closed form,
where and is the modified Bessel function with parameter . This objective is amenable to stochastic optimization where we subsample from the sum to obtain a noisy gradient estimate. We develop a stochastic variational inference scheme by following noisy natural gradients of the variational objective . Using the natural gradient over the standard euclidean gradient is often favorable since natural gradients are invariant to reparameterization of the variational family [21, 22] and provide effective second-order optimization updates [23, 6]. The natural gradients of w.r.t. the Gaussian natural parameters , are
with . Details can be found in the appendix. The natural gradient updates always lead to a positive definite covariance matrix222This follows directly since and are positive definite. and in our implementation has not to be parametrized in any way to ensure positive-definiteness. The derivative of w.r.t. is
Setting it to zero gives the coordinate ascent update for ,
Details can be found in the appendix. The inducing point locations can be either treated as hyperparameters and optimized while training  or can be fixed before optimizing the variational objective. We follow the first approach which is often preferred in a stochastic variational inference setup [7, 10]. The inducing point locations can be either randomly chosen as subset of the training set or via a density estimator. In our experiments we have observed that the -means clustering algorithm (kMeans)  yields the best results. Combining our results, we obtain a fast stochastic variational inference algorithm for the Bayesian nonlinear SVM which is outlined in alg. 1. We apply the adaptive learning rate method described in .
The probabilistic formulation of the SVM lets us directly learn the hyperparameters while training. To this end we maximize the marginal likelihood , where denotes the set of hyperparameters (this approach is called empirical Bayes ). We follow an approximate approach and optimize the fitted variational lower bound over by alternating between optimization steps w.r.t. the variational parameters and the hyperparameters . We include a gradient ascent step w.r.t. after multiple variational updates in the SVI scheme, this is commonly known as Type II maximum likelihood (ML-II) 
Since the standard SVM does not exhibit a probabilistic formulation, the hyperparameters have to be tuned via computationally very expensive methods as grid search and cross validation. Our approach allows us to estimate the hyperparameters during training time and lets us follow gradients instead of only evaluating single hyperparameters.
In the appendix we provide the gradient of the variational objective w.r.t. to a general kernel and show how to optimize arbitrary differentiable hyperparameters. Our experiments exemplify our automated hyperparameter tuning approach by optimizing the hyper parameter of an RBF kernel.
Besides the advantage of automated hyperparameter tuning, the probabilistic formulation of the SVM leads directly to uncertainty estimates of the predictions. The standard SVM lacks this capability, and only heuristic approaches as e.g. Platt  exist. Using the approximate posterior obtained by our stochastic variational inference method (alg. 1) we compute the class membership probability for a test point ,
where denotes the kernel matrix between test and inducing points and the kernel matrix between test points. This leads to the approximate class membership distribution
where is the probit link function. Note that we already computed inverse for the training procedure leading to a computational overhead stemming only from simple matrix multiplication. Our experiments show that (13) leads to reasonable uncertainty estimates.
We now consider the special case of using a linear kernel. If we are interested in this case we may consider the Bayesian model for the linear SVM proposed by Polson et al. (c.f. eq. 4). This can be favorable over using the nonlinear version since this model is formulated in primal space and, therefore, the computational complexity depends on the dimension and not on the number of data points . Furthermore, focusing directly on the linear model allows us to optimize the true ELBO, , without the need of relying on a lower bound (as in eq. 0.A.1). This typically leads to a better approximate posterior.
We again follow the structured mean field approach and chose our variational distributions to be in the same families as the full conditionals (4),
We use again the fact that the coordinate updates of the variational parameters can be obtained by computing the expected natural parameters of the corresponding full conditionals (4) and obtain
where , and . Since the Bayesian Linear SVM model exhibits global and local variables we can directly employ stochastic variational inference by subsampling the data and only updating minibatches of . Note that for the linear case the covariance matrices have size , i.e. being independent of the number of data points. Therefore, the SVI algorithm (14) for the Bayesian Linear SVM exhibits the computational complexity . Luts et. al develop a batch variational inference scheme for the Bayesian linear SVM but do not scale to big datasets.
The hyperparameter can be tuned analogously to (12). The class membership probabilities are
where are the test points and the approximate posterior obtained by the above described SVI scheme.
We compare our approach against the expectation conditional maximization (ECM) method proposed by Henao et. al , Gaussian process classification (GPC) , its recently proposed scalable stochastic variational inference version (S-GPC) , and libSVM with Platt scaling [29, 8] (SVM + Platt).
For all experiments we use an RBF kernel333The RBF kernel is defined as , where is the length scale parameter. with length-scale parameter .
We perform all experiments using only one CPU core with
2.9 GHz and 386 GB RAM.
Code is available at github.com/theogf/BayesianSVM.
We experiment on seven real-world datasets and compare the prediction performance, the quality of the uncertainty estimates and run time of the methods. The results are presented in table 1. We show that our method (S-BSVM) is up to 22 times faster than the direct competitor ECM and up to 700 times faster than Gaussian process classification444For a comparison with the stochastic variational inference version of GPC, see section 5.3. while outperforming the competitors in terms of prediction performance and quality of uncertainty estimates in most cases. The non-probabilistic SVM is naturally the fastest method. Combined with the heuristic Platt scaling approach it leads to class membership probabilities but, however, still lacks the advantages of a probabilistic model (as e.g. uncertainty quantification of the learned parameters and automatic hyperparameter tuning).
To evaluate the quality of the uncertainty estimates we compute the Brier score which is considered as a good performance measure for probabilistic predictions  being defined as , where is the observed output and is the predicted class membership probability. Note that smaller Brier score indicates better performance.
The datasets are all from the Rätsch benchmark datasets  commonly used to test the accuracy of binary nonlinear classifiers. We perform a 10-fold cross-validation and use an RBF kernel with fixed parameters for all methods. For S-BSVM we choose the number of inducing points as of the training set size, except for the datasets Splice, German and Waveform where we use 100 inducing points. For each dataset minibatches of 10 samples are used.
|Dataset||n||dim.||S-BSVM||ECM||GPC||SVM + Platt|
Average prediction error and Brier score with one standard deviation.
We demonstrate the scalability of our method on the SUSY dataset  containing 5 million points with 17 features. This dataset size is very common in particle physics due to the simplicity of artificially generating new events as well as the quantity of data coming from particle detectors. Since it is important to have a sense of the confidence of the predictions for such datasets the Bayesian SVM is an appropriate choice. We use an RBF kernel555The length scale parameter tuning is not included in the training time. We found by our proposed automatic tuning approach., 64 inducing points and minibatches of 100 points. The training of our model takes only 10 minutes without any parallelization. We use the the area under the receiver operating characteristic (ROC) curve (AUC) as performance measure since it is a standard evaluation measure on this dataset .
We examine the run time of our methods and the competitors. We include both the batch variational inference method (B-BSVM) described in section 4.1 and our fast and scalable inference method (S-BSVM) described in section 4.2 in the experiments. For each method we iteratively evaluate the prediction performance on a held-out dataset given a certain training time budget. The prediction error as function of the training time is shown in fig. 1. We experiment on the Waveform dataset from the Rätsch benchmark dataset (). We use an RBF kernel with fixed length-scale parameter and for the stochastic variational inference methods, S-BSVM and S-GPC, we use a batch size of 10 and 100 inducing points.
Our scalable method (S-BSVM) is around 10 times faster than the direct competitor ECM while having slightly better prediction performance. The batch variational inference version (B-BSVM) is the slowest of the Bayesian SVM inference methods. The related probabilistic model, Gaussian process classification, is around 5000 times slower than S-BSVM. Its stochastic inducing point version (S-GPC) has comparable run time to S-BSVM but is very unstable leading to bad prediction performance. S-GPC showed these instabilities for multiple settings of the hyperparameters. The classic SVM (libSVM) has a similar run time as our method. The speed and prediction performance of S-BSVM depend on the number of inducing points. See section 5.5 for an empirical study. Note that the run time in table 1 is determined after the methods have converged.
In section 4.3 we show that our inference method possesses the ability of automatic hyperparameter tunning. In this experiment we demonstrate that our method, indeed, finds the optimal length-scale hyperparameter of the RBF kernel. We use the optimizing scheme (12) and alternate between 10 variational parameter updates and one hyperparameter update. We compute the true validation loss of the length-scale parameter by a grid search approach which consists of training our model (S-BSVM) for each and measuring the prediction performance using 10-fold cross validation. In fig. 2 we plot the validation loss and the length-scale parameter found by our method. We find the true optimum by only using 5 hyperparameter optimization steps. Training and hyperparameter optimization takes only 0.3 seconds for our method, whereas grid search takes 188 seconds (with a grid size of 1000 points).
The sparse GP model used in our inference scheme builds on a set of inducing points where both the number and the locations of the inducing points are free parameters. We investigate three different inducing point selection methods: random subset selection from the training set, the Gaussian Mixture Model (GMM), and the-means clustering algorithm with an improved -means++ seeding (kMeans) . Furthermore we show how the number of inducing points affects the prediction accuracy and the run time. We test the three inducing point selection methods on the USPS dataset  which we reduced to a binary problem using only the digits 3 and 5 (N=1350 and d=256). For all methods we progressively increase the number of inducing points and compute the prediction error by 10-fold cross validation. We present our results in fig. 3.
The GMM is unable to fit large numbers of samples and dimensions and fails to converge for almost all datasets tried, therefore, we do not include it in the plot. Using the -means selection algorithm leads for small numbers of inducing points to much better prediction performance than random subset selection. Furthermore, we show that using only a small fraction of inducing points (around 1% of the original dataset) leads to a nearly optimal prediction performance by simultaneously significantly decreasing the run time. We observe similar results on all datasets we considered.
We presented a fast, scalable and reliable approximate inference method for the Bayesian nonlinear SVM. While previous methods were restricted to rather small datasets our method enables the application of the Bayesian nonlinear SVM to large real world datasets containing millions of samples. Our experiments showed that our method is orders of magnitudes faster than the state-of-the-art while still yielding comparable prediction accuracies. We showed how to automatically tune the hyperparameters and obtain prediction uncertainties which is important in many real world scenarios.
In future work we plan to further extend the Bayesian nonlinear SVM model to deal with missing data and account for correlations between data points building on ideas from . Furthermore, we want to develop Bayesian formulations of important variants of the SVM as for instance one-class SVMs .
We thank Stephan Mandt, Manfred Opper and Patrick Jähnichen for fruitful discussions. This work was partly funded by the German Research Foundation (DFG) award KL 2698/2-1.
Searching for exotic particles in high-energy physics with deep learning.Nature communications (2014) 4308
In: In Artificial Intelligence and Statistics 12. (2009) 567–574
Fast and Provably Good Seedings for k-Means.NIPS (2016)
Sparse Probit Linear Mixed Model.Machine Learning Journal (2017)
Using an Ensemble of One-Class SVM Classifiers to H. P.-based Anomaly Detection Systems.Data Mining (2006)
Using the abbreviation the first expectation term simplifies to
The entropy of is
where is the modified Bessel function with parameter .
By summing the terms the remaining expectations cancel out and we obtain
First, we compute the standard euclidean gradients of . The derivative w.r.t. the mean and covariance matrix are