Support Vector Machine

Support vector machines (SVMs) provide another method for classifying the clothing categories in the Fashion MNIST data. To better understand what SVMs entail, we’ll have to go through some more complex explanations –mainly summarizing James et. al. (2013)– so please bear with me! The figure below might help you in understanding the different classifiers I will discuss in the next sections (figures taken from here, here and here).

For an \(n \times p\) data matrix and binary outcome variable \(y_i \in \{-1, 1\}\), a hyperplane is a flat affine subspace of dimension \(p - 1\) that divides the \(p\)-dimensional space into two halves, defined by \(\beta_0 + \beta_1 X_1 + \dots + \beta_p X_p\). An observation in the test data is assigned an outcome class depending on which side of the perfectly separating hyperplane it lies, assuming that such a hyperplane exists. Cutoff \(t\) for an observation’s score \(\hat{f}(X) = \hat{\beta}_1 X_1 + \hat{\beta}_2 X_2 + \dots + \hat{\beta}_p X_p\) determines which class it will be assigned to. The further an observation is located from the hyperplane at zero, the more confident the classifier is about the class assignment. If existent, an infinite number of separating hyperplanes can be constructed. A good option in this case would be to use the maximal margin classifier (MMC), which maximizes the margin around the midline of the widest strip that can be inserted between the two outcome classes.

If a perfectly separating hyperplane does not exist, “almost separating” hyperplanes can be used by means of the support vector classifier (SVC). The SVC extends the MMC as it does not require classes to be separable by a linear boundary by including slack variables \(\epsilon_i\) that allow some observations to be on the incorrect side of the margin or hyperplane. The extent to which incorrect placements are done is determined by tuning parameter cost \(C \geq \sum_{i=1}^{n} \epsilon_i\), which thereby controls the bias-variance trade-off. The SVC is preferable over the MMC as it is more confident in class assignments due to the larger margins and ensures greater robustness as merely observations on the margin or violating the margin affect the hyperplane (James et al., 2013).

Both MMCs and SVCs assume a linear boundary between the two classes of the outcome variable. Non-linearity can be addressed by enlarging the feature space using functions of predictors. Support vector machines combine SVCs with non-linear (e.g. radial, polynomial or sigmoid) Kernels \(K(x_i, x_{i'})\) to achieve efficient computations. Kernels are generalizations of inner products that quantify the similarity of two observations (James et al., 2013). Usually, the radial Kernel is selected for non-linear models as it provides a good default Kernel in the absence of prior knowledge of invariances regarding translations. The radial Kernel is defined as \(K(x_i, x_{i'})= \exp{(-\sigma \sum_{j=1}^{p} (x_{ij} - x_{i'j})^2)}\), where \(\sigma\) is a positive constant that makes the fit more non-linear as it increases. Tuning \(C\) and \(\sigma\) is necessary to find the optimal trade-off between reducing the number of training errors and making the decision boundary more irregular (by increasing C). As SVMs only require the computation of \(\bigl(\begin{smallmatrix} n\\ 2 \end{smallmatrix}\bigr)\) Kernels for all distinct observation pairs, they greatly improve efficiency.

As aforementioned, the parameters that need to be tuned are cost C and, in the case of a radial Kernel, non-linearity constant sigma. Let’s start by tuning these parameters using a random search algorithm, again making use of the caret framework. We set the controls to perform 5-fold cross-validation and we use the multiClassSummary() function from the MLmetrics library to perform multi-class classification. We specify a radial Kernel, use accuracy as the performance metric¹ and let the algorithm perform a random search for the cost parameter C over pca.dims (=17) random values. Note that the random search algorithm only searches for values of C while keeping a constant value for sigma. Also, contrary to previous calls to trainControl(), we now set classProbs = FALSE because the base package used for estimating SVMs in caret, kernlab, leads to lower accuracies when specifying classProbs = TRUE due to using a secondary regression model (also check this link for the Github issue).

We begin with training the support vector machine using the PCA reduced training and test data sets train.images.pca and test.images.pca constructed in Part 2.

library(MLmetrics)
svm_control = trainControl(method = "repeatedcv",   
                            number = 5,  
                            repeats = 5, 
                            classProbs = FALSE,
                            allowParallel = TRUE, 
                            summaryFunction = multiClassSummary,
                           savePredictions = TRUE)

set.seed(1234)
svm_rand_radial = train(label ~ ., 
                data = cbind(train.images.pca, label = train.classes),
                method = "svmRadial", 
                trControl = svm_control, 
                tuneLength = pca.dims,
                metric = "Accuracy")
svm_rand_radial$results[, c("sigma", "C", 'Accuracy')]

We can check the model performance on both the training and test sets by means of different metrics using a custom function, model_performance, which can be found in this code on my Github.

mp.svm.rand.radial = model_performance(svm_rand_radial, train.images.pca, test.images.pca, 
                                       train.classes, test.classes, "svm_random_radial")

The results show that the model is achieving relatively high accuracies of 88% and 87% on the training and test sets respectively, selecting sigma = 0.040 and C = 32 as the optimal parameters. Let’s have a look at which clothing categories are best and worst predicted by visualizing the confusion matrix. First, let’s compute the predictions for the training data. We need to use the out-of-bag predictions contained in the model object (svm_rand_radial$pred) rather than the manually computed in-sample (non-out-of-bag) predictions for the training data computed using the predict() function. Object svm_rand_radial$pred contains the predictions for all tuning parameter values specified by the user. However, we only need those predictions belonging to the optimal tuning parameter values. Therefore, we subset svm_rand_radial$pred to only contain those predictions and observations in indices rows. Note that we convert svm_rand_radial$pred to a data.table object to find these indices as computations on data.table objects are much faster for large data (e.g. svm_rand_radial$pred has 4.5 million rows).

library(data.table)
pred_dt = as.data.table(svm_rand_radial$pred[, names(svm_rand_radial$bestTune)]) 
names(pred_dt) = names(svm_rand_radial$bestTune)
index_list = lapply(1:ncol(svm_rand_radial$bestTune), function(x, DT, tune_opt){
  return(which(DT[, Reduce(`&`, lapply(.SD, `==`, tune_opt[, x])), .SDcols = names(tune_opt)[x]]))
}, pred_dt, svm_rand_radial$bestTune)
rows = Reduce(intersect, index_list)
pred_train = svm_rand_radial$pred$pred[rows]
trainY = svm_rand_radial$pred$obs[rows]
conf = table(pred_train, trainY)

Next, we reshape the confusion matrix into a data frame with three columns: one for the true categories (trainY), one for the predicted categories (pred_train), and one for the proportion of correct predictions for the true category (Freq). We plot this as a tile plot with a blue color scale where lighter values indicate larger proportions of matches between a particular combination of true and predicted categories, and darker values indicate a small proportion of matches between them. Note that we use the custom plotting theme my_theme() as defined in the second blog post of this series.

conf = data.frame(conf / rowSums(conf))
ggplot() + 
  geom_tile(data = conf, aes(x = trainY, y = pred_train, fill = Freq)) + 
  labs(x = "Actual", y = "Predicted", fill = "Proportion") +
  my_theme() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
   scale_fill_continuous(breaks = seq(0, 1, 0.25)) +
  coord_fixed()

We observe from this plot that most of the classes are predicted accurately as the light blue (high percentages of correct predictions) are on the diagonal of the tile plot. We can also observe that the categories that are most often mixed up include shirts, tops, pullovers and coats, which makes sense because these are all mostly upper body clothing parts having similar shapes. The model predicts trousers, bags, boots and sneakers well, given that these rows and columns are particularly dark except for the diagonal element. These results are in agreement with those from the random forest and gradient-boosted trees from the previous blog post of this series.

Next, we repeat the above process for fitting a support vector machine but instead of a random search for the optimal parameters, we perform a grid search. As such, we can prespecify values to evaluate the model at, not only for C but also for sigma. We define the grid values in grid_radial.

grid_radial = expand.grid(sigma = c(.01, 0.04, 0.1), C = c(0.01, 10, 32, 70, 150))
set.seed(1234)
svm_grid_radial = train(label ~ ., 
                              data = cbind(train.images.pca, label = train.classes),
                              method = "svmRadial", 
                              trControl = svm_control, 
                              tuneGrid = grid_radial,
                              metric = "Accuracy")
svm_grid_radial$results[, c("sigma", "C", 'Accuracy')]

mp.svm.grid.radial = model_performance(svm_grid_radial, train.images.pca, test.images.pca, 
                                       train.classes, test.classes, "svm_grid_radial")

The grid search selects the same optimal parameter values as the random search (C=32 and sigma = 0.040), therefore also resulting in 88% and 87% training and test accuracies. To get an idea on how C and sigma influence the training set accuracy, we plot the cross-validation accuracy as a function of C, with separate lines for each value of sigma.

ggplot() + 
  my_theme() +
  geom_line(data = svm_grid_radial$results, aes(x = C, y = Accuracy, color = factor(sigma))) +
  geom_point(data = svm_grid_radial$results, aes(x = C, y = Accuracy, color = factor(sigma))) +
  labs(x = "Cost", y = "Cross-Validation Accuracy", color = "Sigma") +
  ggtitle('Relationship between cross-validation accuracy and values of cost and sigma')

The plot shows that the green line (sigma = 0.04) has the highest cross-validation accuracy for all values of C except for smaller values of C such as 0.01 and 10. Although the accuracy at C=10 and sigma = 0.1 (blue line) comes close, the highest overall accuracy achieved is for C=32 and sigma=32 (green line).

Wrapping Up

To compare the models we have estimated throughout this series of blog posts, we can look at the resampled accuracies of the models. We can do this in our case because we set the same seed of 1234 before training each model.² Essentially, resampling is an important tool to validate our models, and to what extent they are generalizeable onto data they have not been trained on. We used five repeats of five-fold cross-validation, which means that the training data was divided into five random subsets, and that throughout five iterations (“folds”) the model was trained on four of these subsets and tested on the remaining subset (changing with every fold), and that this whole process has been repeated five times. The goal of these repetitions of k-fold cross-validation is to reduce the bias in the estimator, given that the folds in non-repeated cross-validation are not independent (as data used for training at one fold is used for testing at another fold). As we performed five repeats of five-fold cross-validation, we can essentially obtain 5*5=25 accuracies per model. Let’s compare these resampled accuracies visually by means of a boxplot. First, we create a list of all models estimated, including the random forests, gradient-boosted trees and support vector machines. We then compute the resampled accuracies using the resamples() function from the caret package. From the resulting object, resamp, we only keep the columns containing the resample unit (e.g. Fold1.Rep1) and the five columns containing the accuracies for each of the five models. We melt this into a long format and from the result, plotdf, we remove the ~Accuracy part from the strings in column Model.

library(reshape2)
model_list = list(rf_rand, rf_grid, xgb_tune, svm_rand_radial, svm_grid_radial)
names(model_list) = c(paste0('Random forest ', c("(random ", "(grid "), "search)"), "Gradient-boosted trees", 
                      paste0('Support vector machine ', c("(random ", "(grid "), "search)"))
resamp = resamples(model_list)
accuracy_variables = names(resamp$values)[grepl("Accuracy", names(resamp$values))]
plotdf = melt(resamp$values[, c('Resample', accuracy_variables)], 
              id = "Resample", value.name = "Accuracy", variable.name = "Model")
plotdf$Model = gsub("~.*","", plotdf$Model)

Next, we create a boxplot with the estimated models on the x-axis and the accuracy on the y-axis.

ggplot() +
  geom_boxplot(data = plotdf, aes(x = Model, y = Accuracy, color = Model)) +
  ggtitle('Resampled accuracy for machine learning models estimated') + 
  my_theme() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  labs(x = NULL, color = NULL) +
  guides(color = FALSE)

We observe from these box plots that the support vector machines perform best, followed by the gradient-boosted trees and the random forests. Let’s also take a look at the other performance metrics from all models we have looked at.

mp.df = rbind(mp.rf.rand, mp.rf.grid, mp.xgb, mp.svm.rand.radial, mp.svm.grid.radial, mp.svm.grid.linear)
mp.df[order(mp.df$accuracy_test, decreasing = TRUE), ]

After taking measures to reduce overfitting, the convolutional neural network from the first blog post of this series achieved training and test set accuracies of 89.4% and 88.8% respectively. The random and grid search for the best value of mtry in the random forests resulted in the selection of mtry=5. The grid search performed better on the training set than the random search on the basis of all metrics except recall (i.e. sensitivity), and better on the test set on all metrics except precision (i.e. positive predictive value). The test set accuracies achieved by the random search and grid search were 84.7% and 84.8% respectively. The gradient-boosted decision trees performed slightly better than the random forests on all metrics and achieved a test set accuracy of 85.5%. Both tree-based models more often misclassified pullovers, shirts and coats, while correctly classifying trousers, boots, bags and sneakers. The random forests and gradient-boosted trees are however outperformed by the support vector machine with radial Kernel specification with tuning parameter values of C=32 and sigma=0.040: this model achieved 86.9% test set accuracy upon a random search for the best parameters. The grid search resulted in slightly worse test set performance, but better training set performance in terms of all metrics except accuracy. Nonetheless, none of the models estimated beats the convolutional neural network from the first blog post of this series, neither in performance nor computational time and feasibility. However, the differences in test set performance are only small: the convolutional neural network achieved 88.8% test set accuracy, compared to 86.9% test set accuracy achieved by the support vector machine with radial Kernel. This shows that we do not always need to resort to deep learning to obtain high accuracies, but that we can also perform image classification to a reasonable standard using basic machine learning models with dimensionality-reduced data.