Stochastic Subgradient Descent + Linear SVM example in R

ssd_svm.R

# Setting up random matrix
set.seed(11111)
x <- data.frame(a = rnorm(n = 100) * 5,
                b = rnorm(n = 100) * 3 + 1,
                c = rnorm(n = 100) * 2 + 2)

# Setting up an imperfect relationship (non-linear problem) with +1 or -1 for classes
y <- 2 * as.integer((2 + (x[, 1] * 2) + (x[, 2] * 3) + (x[, 3] * 4) + (x[, 3] ^ 2) + (x[, 1] * x[, 2])) > 20) - 1

# Setting up polynomial features
columns <- ncol(x)
for (i in 1:columns) {
  x[, paste0(colnames(x)[i], "X", colnames(x)[i])] <- x[, i] * x[, i]
  for (j in i:columns) {
    x[, paste0(colnames(x)[i], "X", colnames(x)[j])] <- x[, i] * x[, j]
  }
}

# Add column names and intercept
colnames(x) <- c("a*2", "b*3", "c*4", "aXa", "aXb*1", "aXc", "bXb", "bXc", "cXc*1")
x <- as.matrix(x)

cost <- function(x, y, a, b){
  threshold <- x %*% a + b
  predicted <- rep(0,length(y))
  predicted[threshold < 0] <- -1
  predicted[threshold >= 0] <- 1
  return(sum(predicted == y) / length(y))
}
grad <- function(x, y, a, b, lambda) {
  
  hard_margin <- y * (x %*% a + b) # define hard-margin SVM
  
  if (hard_margin >= 1) {
    # point is not within the hard-margin, so increase the margin
    gradient <- c(b, a * lambda)
  } else {
    # point is within the hard-margin, so decrease the margin
    gradient <- c(-y, (a * lambda) - (y * x))
  }
  
  return(gradient)
}

SVM_Regularization <- function(x, y, init, epochs, eta, cost, grad, lambda) {
  
  iters <- epochs * length(y)
  param <- data.frame(matrix(nrow = iters + 1, ncol = length(init) + 2))
  colnames(param) <- c("Intercept", colnames(x), "Loss")
  param[1, ] <- c(c(0, init), cost(x, y, init, 0))
  
  for (i in 1:iters) {
    set.seed(i)
    sample_used <- sample.int(length(y), 1)
    param[i + 1, 1:(length(init) + 1)] <- as.numeric(param[i, 1:(length(init) + 1)]) - eta * grad(x[sample_used, ], y[sample_used], as.numeric(param[i, 2:(length(init) + 1)]), param[i, 1], lambda)
    param[i + 1, length(init) + 2] <- cost(x, y, as.numeric(param[i + 1, 2:(length(init) + 1)]), param[i + 1, 1])
  }
  
  cat("Final cost: ", sprintf("%10.07f", param[nrow(param), ncol(param)]), "\n", sep = "")
  cat("Parameters:", as.numeric(param[nrow(param), 1:(length(init) + 1)]), sep = " ")
  
  param <- cbind(Iteration = 0:(nrow(param) - 1), param)
  
  return(param)
  
}

param <- SVM_Regularization(x = x,
                            y = y,
                            init = rep(0, 9),
                            epochs = 200,
                            eta = 0.0001,
                            cost = cost,
                            grad = grad,
                            lambda = 0.01)

melted_param <- reshape::melt(param[(1:((nrow(param) - 1) / 10)) * 10 + 1, 1:(ncol(param) - 1)], id = c("Iteration"))
#car::scatterplot(value ~ Iteration | variable, data = melted_param)
lattice::xyplot(value ~ Iteration | variable, data = melted_param, type = "l", panel = function(...) {panel.xyplot(...); panel.abline(h = 0)})

melted_loss <- reshape::melt(param[(1:((nrow(param) - 1) / 10)) * 10 + 1, c(1,ncol(param))], id = c("Iteration"))
ggplot2::ggplot(melted_loss, ggplot2::aes(x = Iteration, y = value)) + ggplot2::geom_line() + ggplot2::theme_bw() + ggplot2::labs(title = "Evolution of loss by number of boosting iterations\n(no iterations skipped)", x = "Iteration", y = "Accuracy")