Support Vector Machines & Kernel Methods (R) — e1071 · kernlab

Max-margin classification, the kernel trick, and kernel ridge = Gaussian process

e1071::svm (the R interface to LIBSVM) and kernlab are the reference R kernel-machine packages — the cross-check for the from-scratch Pegasos SVM and kernel ridge in the Python notebook. We fit a linear SVM (max-margin, hinge loss) against a row-matched logistic, show the kernel trick turning a linear machine nonlinear on the two-moons data, fit SVR and kernel ridge on California, and verify the bridge to the Bayesian Nonparametrics section: kernel ridge = the posterior mean of a Gaussian process (kernlab::gausspr), the same estimator benchmarked in Bayesian Nonparametric Prediction and developed at length in Gaussian-Process Regression. SVMs are $O(n^2)$, so the credit/California fits use subsamples (as the GP did). ROC-AUC is defined in the CART notebook.

options(repr.plot.width=12.5, repr.plot.height=4.8)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(e1071); library(kernlab)})
auc <- function(y,p){ r<-rank(p); n1<-sum(y==1); n0<-sum(y==0); (sum(r[y==1])-n1*(n1+1)/2)/(n1*n0) }
rmse <- function(a,b) sqrt(mean((a-b)^2))
cat("e1071", as.character(packageVersion("e1071")), "| kernlab", as.character(packageVersion("kernlab")), "\n")

e1071 1.7.17 | kernlab 0.9.33 

1. Max-margin — a linear SVM

The SVM seeks the widest-margin separating hyperplane under the hinge loss; only the support vectors (points on or inside the margin) shape it. The left panel shows the boundary, margins and support vectors on a 2-D toy; then we fit a linear SVM to the credit data and compare its out-of-sample AUC to a logistic baseline fit to the same rows, so that subsampling cannot be mistaken for a difference between the losses.

set.seed(0)
# 2-D toy: two gaussian blobs
n0<-120; A<-cbind(rnorm(n0,-1.3,0.85),rnorm(n0,-1.3,0.85)); B<-cbind(rnorm(n0,1.3,0.85),rnorm(n0,1.3,0.85))
Xt<-rbind(A,B); yt<-factor(rep(c(0,1),each=n0)); dt<-data.frame(x1=Xt[,1],x2=Xt[,2],y=yt)
sv<-svm(y~., data=dt, kernel="linear", cost=1, scale=FALSE)
par(mfrow=c(1,2), mar=c(4,4,3,1))
g<-expand.grid(x1=seq(-4,4,length=200), x2=seq(-4,4,length=200))
dv<-attr(predict(sv,g,decision.values=TRUE),"decision.values")[,1]
image(seq(-4,4,length=200),seq(-4,4,length=200),matrix(dv>0,200), col=c("#2b6cb022","#c5303022"),
      xlab="x1", ylab="x2", main="Max-margin linear SVM (e1071)")
contour(seq(-4,4,length=200),seq(-4,4,length=200),matrix(dv,200), levels=c(-1,0,1), lty=c(2,1,2), lwd=c(1,2,1), add=TRUE)
points(Xt, col=ifelse(yt==0,"#2b6cb0","#c53030"), pch=19, cex=0.6)
points(Xt[sv$index,], col="#2f855a", cex=1.6, lwd=1.5)
legend("topleft", c("support vectors"), col="#2f855a", pch=1, bty="n", cex=0.8)
# credit
d<-read.csv("credit_default.csv"); feat<-setdiff(names(d),"default")
set.seed(0); i<-sample(nrow(d),0.7*nrow(d)); tr<-d[i,]; te<-d[-i,]
sc<-scale(tr[,feat]); ctr<-as.data.frame(sc); ctr$y<-factor(tr$default)
cte<-as.data.frame(scale(te[,feat], center=attr(sc,"scaled:center"), scale=attr(sc,"scaled:scale")))
sub<-sample(nrow(ctr),3000)
lsv<-svm(y~., data=ctr[sub,], kernel="linear", cost=1, probability=TRUE)
plin<-attr(predict(lsv, cte, probability=TRUE),"probabilities")[,"1"]
svm_lin_auc<-auc(te$default, plin)
glm_auc<-auc(te$default, predict(glm(default~., tr, family=binomial), te, type="response"))
gsub_df<-ctr[sub,feat]; gsub_df$default<-tr$default[sub]                 # logistic on the SVM's own 3,000 rows
glm_sub_auc<-auc(te$default, predict(glm(default~., gsub_df, family=binomial), cte, type="response"))
plot.new(); text(0.5,0.5, sprintf("Credit default -- test AUC\n\n linear SVM (e1071, 3k) : %.3f\n logistic (glm, 3k)     : %.3f\n logistic (glm, full)   : %.3f\n\nThe linear SVM lands just under the\nlogistic -- and not for want of rows:\nmatched on the same 3,000, the logistic\nstill wins. Max-margin buys nothing when\nthe classes overlap and the boundary\nis near-linear.", svm_lin_auc, glm_sub_auc, glm_auc), cex=1.05, family="mono")
par(mfrow=c(1,1))
cat(sprintf("linear SVM %.3f against a logistic fit to the same 3,000 rows %.3f (and %.3f on all %d) -- the max-margin\n", svm_lin_auc, glm_sub_auc, glm_auc, nrow(tr)))
cat(sprintf("loss gives up about %.3f AUC here, and subsampling is not the reason: the logistic barely moves between\n", glm_sub_auc-svm_lin_auc))
cat(sprintf("3,000 rows and 21,000. Only the %d support vectors set the boundary; the rest could be moved freely.\n", lsv$tot.nSV))

Warning message:
"glm.fit: fitted probabilities numerically 0 or 1 occurred"

Warning message:
"glm.fit: fitted probabilities numerically 0 or 1 occurred"

linear SVM 0.699 against a logistic fit to the same 3,000 rows 0.713 (and 0.712 on all 21000) -- the max-margin

loss gives up about 0.014 AUC here, and subsampling is not the reason: the logistic barely moves between

3,000 rows and 21,000. Only the 1343 support vectors set the boundary; the rest could be moved freely.

2. The kernel trick — going nonlinear

Swapping the inner product for an RBF kernel $k(x,x')=e^{-\sigma\lVert x-x'\rVert^2}$ lifts the data into an implicit high-dimensional space where a linear margin becomes a curved boundary. On the two-moons data the linear SVM fails and the radial SVM separates cleanly. On the real credit data the RBF kernel adds little — overlapping, imbalanced classes offer scant nonlinear structure.

set.seed(1); n<-200; t<-runif(n,0,pi)
M<-rbind(cbind(cos(t),sin(t)), cbind(1-cos(t),1-sin(t)-0.5)) + matrix(rnorm(4*n,0,0.14),2*n,2)
dm<-data.frame(x1=M[,1],x2=M[,2],y=factor(rep(c(0,1),each=n)))
lin<-svm(y~.,dm,kernel="linear",cost=1,scale=FALSE); rad<-svm(y~.,dm,kernel="radial",cost=1,gamma=1,scale=FALSE)
gm<-expand.grid(x1=seq(-1.6,2.6,length=200), x2=seq(-1.4,1.9,length=200))
par(mfrow=c(1,2), mar=c(4,4,3,1))
for(pair in list(list(lin,"Linear SVM -- cannot separate"), list(rad,"RBF-kernel SVM -- curved margin"))){
  mdl<-pair[[1]]; z<-attr(predict(mdl,gm,decision.values=TRUE),"decision.values")[,1]
  image(seq(-1.6,2.6,length=200),seq(-1.4,1.9,length=200),matrix(z>0,200), col=c("#2b6cb022","#c5303022"),
        xlab="x1",ylab="x2",main=sprintf("%s (acc %.2f)", pair[[2]], mean(predict(mdl,dm)==dm$y)))
  contour(seq(-1.6,2.6,length=200),seq(-1.4,1.9,length=200),matrix(z,200), levels=0, lwd=2, add=TRUE)
  points(M, col=ifelse(dm$y==0,"#2b6cb0","#c53030"), pch=19, cex=0.5) }
par(mfrow=c(1,1))
rsv<-svm(y~., data=ctr[sub,], kernel="radial", cost=1, gamma=0.02, probability=TRUE)
prad<-attr(predict(rsv, cte, probability=TRUE),"probabilities")[,"1"]; svm_rbf_auc<-auc(te$default, prad)
cat(sprintf("RBF SVM on credit (subsample 3000): AUC %.3f, against the linear SVM's %.3f and the logistic's %.3f. The\n", svm_rbf_auc, svm_lin_auc, glm_sub_auc))
cat("kernel recovers part of the gap and finds no more -- overlapping, imbalanced classes hold little smooth structure,\n")
cat("and the curved boundary costs O(n^2) to buy it. The kernel earns its keep on the regression task instead.\n")

RBF SVM on credit (subsample 3000): AUC 0.706, against the linear SVM's 0.699 and the logistic's 0.713. The

kernel recovers part of the gap and finds no more -- overlapping, imbalanced classes hold little smooth structure,

and the curved boundary costs O(n^2) to buy it. The kernel earns its keep on the regression task instead.

3. Kernel ridge & SVR — California, and the GP bridge

For regression the kernel attaches to ridge regression: kernel ridge solves $\alpha=(K+\lambda I)^{-1}y$. We compute it by hand with a kernlab RBF kernel matrix, fit SVR (e1071, $\varepsilon$-insensitive) alongside, and — the key identity — verify that kernel ridge equals kernlab::gausspr, a Gaussian process with matched noise variance. Both beat the linear model; and kernel ridge and the GP produce the same predictions, differing only in that the GP also returns uncertainty. Fit on a 1,500-point subsample.

h<-read.csv("cali_housing.csv"); hf<-setdiff(names(h),"MedHouseVal")
set.seed(0); j<-sample(nrow(h),0.7*nrow(h)); htr<-h[j,]; hte<-h[-j,]
Rtr<-scale(htr[,hf]); Rte<-scale(hte[,hf], center=attr(Rtr,"scaled:center"), scale=attr(Rtr,"scaled:scale"))
Rtr<-matrix(as.numeric(Rtr),nrow(htr)); Rte<-matrix(as.numeric(Rte),nrow(hte))
set.seed(1); s2<-sample(nrow(Rtr),1500); Xs<-Rtr[s2,]; ys<-htr$MedHouseVal[s2]
ym<-mean(ys); ysd<-sd(ys); yz<-(ys-ym)/ysd
sig<-0.125                                             # RBF sigma (1/n_features)
rbf<-rbfdot(sigma=sig); K<-kernelMatrix(rbf, Xs); lambda<-1
alpha<-solve(K + lambda*diag(nrow(Xs)), yz)            # kernel ridge, by hand
Ktest<-kernelMatrix(rbf, Rte, Xs); pred_kr<-as.numeric(Ktest %*% alpha)*ysd + ym
gp<-gausspr(Xs, yz, kernel="rbfdot", kpar=list(sigma=sig), var=lambda, scaled=FALSE)  # GP, matched noise
pred_gp<-as.numeric(predict(gp, Rte))*ysd + ym
svr<-svm(x=Xs, y=ys, type="eps-regression", kernel="radial", gamma=sig, cost=1)
svr_ns<-svm(x=Xs, y=ys, type="eps-regression", kernel="radial", gamma=sig, cost=1, scale=FALSE)
svr_ns_rmse<-rmse(hte$MedHouseVal, as.numeric(predict(svr_ns, Rte)))
pred_svr<-as.numeric(predict(svr, Rte))
kr_rmse<-rmse(hte$MedHouseVal,pred_kr); svr_rmse<-rmse(hte$MedHouseVal,pred_svr)
cat(sprintf("California test RMSE -- kernel ridge %.3f  SVR %.3f  (linear model 0.737)\n", kr_rmse, svr_rmse))
cat(sprintf("kernel ridge vs GP (gausspr) posterior mean: max|diff| %.2e, correlation %.5f -- the SAME estimator.\n",
    max(abs(pred_kr-pred_gp)), cor(pred_kr,pred_gp)))
cat(sprintf("SVR beats kernel ridge here by %.3f, which the Python notebook does not see. The cause is a default: e1071\n", kr_rmse-svr_rmse))
cat(sprintf("scales x and y unless told otherwise, and with scale=FALSE the identical call gives %.3f instead. That is not\n", svr_ns_rmse))
cat("cosmetic. An isotropic RBF kernel measures distance in whatever units it is handed, so rescaling the features IS a\n")
cat("change of kernel. AveOccup is what makes it bite: its full-sample SD is set by a handful of extreme blocks, so the\n")
cat("standardised subsample is nearly constant along that axis (SD 0.11) until e1071 restandardises it back to 1.\n")
par(mfrow=c(1,2), mar=c(4,4,3,1))
plot(pred_kr, hte$MedHouseVal, pch=".", col="#a0aec0", xlab="kernel-ridge predicted ($100k)", ylab="actual",
     main=sprintf("California -- kernel ridge (RMSE %.3f)", kr_rmse)); abline(0,1,lty=2)
qb<-quantile(pred_kr,seq(0,1,length=11)); bb<-cut(pred_kr,unique(qb),include.lowest=TRUE)
points(tapply(pred_kr,bb,mean), tapply(hte$MedHouseVal,bb,mean), type="b", pch=19, col="#c53030", lwd=2)
plot(pred_kr, pred_gp, pch=".", col="#6b46c1", xlab="kernel ridge", ylab="GP posterior mean (gausspr)",
     main=sprintf("Same prediction: max|diff| %.1e", max(abs(pred_kr-pred_gp)))); abline(0,1,lty=2)
par(mfrow=c(1,1))

California test RMSE -- kernel ridge 0.655  SVR 0.601  (linear model 0.737)

kernel ridge vs GP (gausspr) posterior mean: max|diff| 5.46e-14, correlation 1.00000 -- the SAME estimator.

SVR beats kernel ridge here by 0.054, which the Python notebook does not see. The cause is a default: e1071

scales x and y unless told otherwise, and with scale=FALSE the identical call gives 0.656 instead. That is not

cosmetic. An isotropic RBF kernel measures distance in whatever units it is handed, so rescaling the features IS a

change of kernel. AveOccup is what makes it bite: its full-sample SD is set by a handful of extreme blocks, so the

standardised subsample is nearly constant along that axis (SD 0.11) until e1071 restandardises it back to 1.

4. Scoreboard and summary

Support-vector and kernel methods beside the running field. Linear and RBF SVM, kernel ridge and SVR are computed here, and the logistic and linear baselines are refitted; the forest, XGBoost, GP and GAM figures are carried from the R notebooks, which share this notebook's set.seed(0) 70/30 split. Python's corresponding numbers come from a different split and are not interchangeable with these.

lin_rmse<-rmse(hte$MedHouseVal, predict(lm(MedHouseVal~., htr), hte))
# carried figures are the R notebooks' own, on this same set.seed(0) 70/30 split -- not the Python ones,
# which come from a different split and would make the bars incomparable.
credit<-c("logistic"=glm_auc, "linear SVM"=svm_lin_auc, "RBF SVM"=svm_rbf_auc,
          "random forest"=0.774, "XGBoost"=0.776, "GP (R BNP nb)"=0.739, "GAM (R BNP nb)"=0.743)
cali  <-c("linear"=lin_rmse, "kernel ridge"=kr_rmse, "SVR"=svr_rmse,
          "random forest"=0.491, "XGBoost"=0.484, "GP (R BNP nb)"=0.611, "GAM (R BNP nb)"=0.607)
kerc<-c("#dd6b20","#dd6b20");
par(mfrow=c(1,2), mar=c(4,7,3,1))
colc<-ifelse(names(credit)%in%c("linear SVM","RBF SVM"),"#dd6b20", ifelse(grepl("BNP",names(credit)),"#6b46c1","#a0aec0"))
barplot(rev(credit), horiz=TRUE, las=1, xlim=c(0.5,0.8), col=rev(colc), main="Credit -- test AUC (higher better)")
cola<-ifelse(names(cali)%in%c("kernel ridge","SVR"),"#dd6b20", ifelse(grepl("BNP",names(cali)),"#6b46c1","#a0aec0"))
barplot(rev(cali), horiz=TRUE, las=1, col=rev(cola), main="California -- test RMSE (lower better)")
par(mfrow=c(1,1))
cat("Kernel methods (orange) land just under the logistic on credit and clear the linear model on California -- where\n")
cat("kernel ridge IS the GP mean. Trees and boosting still lead raw accuracy at full scale, and kernels pay the O(n^2)\n")
cat("cost that forces subsampling -- the standing reason they do not own large tabular problems.\n")

Kernel methods (orange) land just under the logistic on credit and clear the linear model on California -- where

kernel ridge IS the GP mean. Trees and boosting still lead raw accuracy at full scale, and kernels pay the O(n^2)

cost that forces subsampling -- the standing reason they do not own large tabular problems.

Are the Platt-scaled probabilities honest?

e1071 returns probabilities only when asked (probability=TRUE), and it obtains them by Platt scaling — a logistic fitted to the SVM's decision values by internal cross-validation, since a margin classifier has no probability of its own. The reliability curve tests whether that retro-fit produces something usable: bin by predicted probability, plot the observed default rate, and compare against a logistic that estimates one directly.

rel <- function(p, y, lab, col) {
  q <- unique(quantile(p, seq(0,1,length=11))); b <- cut(p, q, include.lowest=TRUE)
  pp <- tapply(p,b,mean); oo <- tapply(y,b,mean); w <- as.numeric(table(b))/length(y)
  ece <- sum(w*abs(oo-pp), na.rm=TRUE)
  list(pp=as.numeric(pp), oo=as.numeric(oo), ece=ece, lab=lab, col=col)
}
draw <- function(L, main) {
  xm <- max(sapply(L, function(z) max(z$pp))); ym <- max(sapply(L, function(z) max(z$oo)))
  plot(NA, xlim=c(0,xm), ylim=c(0,max(xm,ym)), xlab="predicted P(default)",
       ylab="observed default rate", main=main); abline(0,1,lty=2)
  for (z in L) lines(z$pp, z$oo, type="b", pch=19, col=z$col, lwd=2)
  legend("topleft", sapply(L, function(z) sprintf("%s (ECE %.3f)", z$lab, z$ece)),
         col=sapply(L, function(z) z$col), lwd=2, pch=19, bty="n", cex=0.8)
}
options(repr.plot.width=7, repr.plot.height=5); par(mar=c(4,4,3,1))
pglm <- predict(glm(default~., gsub_df, family=binomial), cte, type="response")
L <- list(rel(plin, te$default, "linear SVM + Platt", "#dd6b20"),
          rel(prad, te$default, "RBF SVM + Platt",    "#2b6cb0"),
          rel(pglm, te$default, "logistic",           "#a0aec0"))
draw(L, "Reliability of the Platt-scaled SVM scores")
for (z in L) cat(sprintf("%-22s ECE %.3f  mean predicted %.3f  base rate %.3f\n",
                         z$lab, z$ece, mean(get(c("plin","prad","pglm")[which(sapply(L,function(q) q$lab)==z$lab)])),
                         mean(te$default)))
cat("\nAn SVM has no probability of its own -- Platt scaling fits one to the decision values afterwards, so the\n")
cat("reliability curve is the check on whether that retro-fit is usable. Compare the ECE figures with the AUC\n")
cat("ranking earlier in the notebook: they measure different things, and a model can rank well while its fitted\n")
cat("probabilities drift. For triage the ranking is enough; for anything that feeds an expected-loss calculation\n")
cat("the diagonal is what matters.\n")

Warning message:
"glm.fit: fitted probabilities numerically 0 or 1 occurred"

linear SVM + Platt     ECE 0.054  mean predicted 0.220  base rate 0.221
RBF SVM + Platt        ECE 0.030  mean predicted 0.219  base rate 0.221
logistic               ECE 0.052  mean predicted 0.217  base rate 0.221

An SVM has no probability of its own -- Platt scaling fits one to the decision values afterwards, so the

reliability curve is the check on whether that retro-fit is usable. Compare the ECE figures with the AUC

ranking earlier in the notebook: they measure different things, and a model can rank well while its fitted

probabilities drift. For triage the ranking is enough; for anything that feeds an expected-loss calculation

the diagonal is what matters.

Summary

e1071 and kernlab reproduce the Python notebook: a max-margin linear SVM that lands just under the logistic on credit even when the two are matched row for row, the RBF kernel bending a linear machine around the two-moons data, and kernel ridge / SVR beating the linear model on California (0.74 → 0.65). The headline identity is verified in R too — kernel ridge equals kernlab::gausspr, a Gaussian process — so the frequentist kernel machine and the Bayesian nonparametric one are the same estimator, differing only in the GP's extra uncertainty. Kernels share the GP's $O(n^2\!-\!n^3)$ cost, the standing reason boosting dominates large tabular data.

This is the package mirror of svm_kernels_python.ipynb, closing the Regularized & Kernel Learning subsection: penalised linear models and their kernelised, max-margin relatives, cross-linked throughout to their Bayesian counterparts (shrinkage priors and Gaussian processes).