{"word_count": 1347, "character_count": 15343}
Seventh post of our series on classification from scratch. The latest one was on the neural nets, and today, we will discuss SVM, support vector machines.A formal introductionHere [latex]y[/latex] takes values in [latex]\{-1,+1\}[/latex]. Our model will be [latex display="true"]m(\mathbf{x})=\text{sign}[\mathbf{\omega}^T\mathbf{x}+b][/latex] Thus, the space is divided by a (linear) border[latex display="true"]\Delta:\lbrace\mathbf{x}\in\mathbb{R}^p:\mathbf{\omega}^T\mathbf{x}+b=0\rbrace[/latex]The distance from point [latex]\mathbf{x}_i[/latex] to [latex]\Delta[/latex] is [latex display="true"]d(\mathbf{x}_i,\Delta)=\frac{\mathbf{\omega}^T\mathbf{x}_i+b}{\|\mathbf{\omega}\|}[/latex]If the space is linearly separable, the problem is ill posed (there is an infinite number of solutions). So consider[latex display="true"]\max_{\mathbf{\omega},b}\left\lbrace\min_{i=1,\cdots,n}\left\lbrace\text{distance}(\mathbf{x}_i,\Delta)\right\rbrace\right\rbrace[/latex]The strategy is to maximize the margin. One can prove that we want to solve [latex display="true"]\max_{\mathbf{\omega},m}\left\lbrace\frac{m}{\|\mathbf{\omega}\|}\right\rbrace[/latex]subject to [latex]y_i\cdot(\mathbf{\omega}^T\mathbf{x}_i)=m[/latex], [latex]\forall i=1,\cdots,n[/latex]. Again, the problem is ill posed (non identifiable), and we can consider [latex]m=1[/latex]: [latex display="true"]\max_{\mathbf{\omega}}\left\lbrace\frac{1}{\|\mathbf{\omega}\|}\right\rbrace[/latex]subject to [latex]y_i\cdot(\mathbf{\omega}^T\mathbf{x}_i)=1[/latex], [latex]\forall i=1,\cdots,n[/latex]. The optimization objective can be written[latex display="true"]\min_{\mathbf{\omega}}\left\lbrace\|\mathbf{\omega}\|^2\right\rbrace[/latex]The primal problemIn the separable case, consider the following primal problem,[latex display="true"]\min_{\mathbf{w}\in\mathbb{R}^d,b\in\mathbb{R}}\left\lbrace\frac{1}{2}\|\mathbf{\omega}\|^2\right\rbrace[/latex]subject to [latex]y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1[/latex], [latex]\forall i=1,\cdots,n[/latex].In the non-separable case, introduce slack (error) variables [latex]\mathbf{\xi}[/latex] : if [latex]y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1[/latex], there is no error [latex]\xi_i=0[/latex].Let [latex]C[/latex] denote the cost of misclassification. The optimization problem becomes[latex display="true"]\min_{\mathbf{w}\in\mathbb{R}^d,b\in\mathbb{R},{\color{red}{\mathbf{\xi}}}\in\mathbb{R}^n}\left\lbrace\frac{1}{2}\|\mathbf{\omega}\|^2 + C\sum_{i=1}^n\xi_i\right\rbrace[/latex]subject to [latex]y_i\cdot (\mathbf{\omega}^T\mathbf{x}_i+b)\geq 1-{\color{red}{\xi_i}}[/latex], with [latex]{\color{red}{\xi_i}}\geq 0[/latex], [latex]\forall i=1,\cdots,n[/latex].Let us try to code this optimization problem. The dataset is heren = length(myocarde[,"PRONO"])myocarde0 = myocardemyocarde0$PRONO = myocarde$PRONO*2-1C = .5and we have to set a value for the cost [latex]C[/latex]. In the (linearly) constrained optimization function in R, we need to provide the objective function [latex]f(\mathbf{\theta})[/latex] and the gradient [latex]\nabla f(\mathbf{\theta})[/latex].f = function(param){  w  = param[1:7]  b  = param[8]  xi = param[8+1:nrow(myocarde)]  .5*sum(w^2) + C*sum(xi)}grad_f = function(param){  w  = param[1:7]  b  = param[8]  xi = param[8+1:nrow(myocarde)]  c(2*w,0,rep(C,length(xi)))}and (linear) constraints are written as [latex]\mathbf{U}\mathbf{\theta}-\mathbf{c}\geq \mathbf{0}[/latex]U = rbind(cbind(myocarde0[,"PRONO"]*as.matrix(myocarde[,1:7]),diag(n),myocarde0[,"PRONO"]),cbind(matrix(0,n,7),diag(n,n),matrix(0,n,1)))C = c(rep(1,n),rep(0,n))Then we useconstrOptim(theta=p_init, f, grad_f, ui = U,ci = C)Observe that something is missing here: we need a starting point for the algorithm, [latex]\mathbf{\theta}_0[/latex]. Unfortunately, I could not think of a simple technique to get a valid starting point (that satisfies those linear constraints).Let us try something else. Because those functions are quite simple: either linear or quadratic. Actually, one can recognize in the separable case, but also in the non-separable case, a classic quadratic program[latex display="true"]\min_{\mathbf{z}\in\mathbb{R}^d}\left\lbrace\frac{1}{2}\mathbf{z}^T\mathbf{D}\mathbf{z}-\mathbf{d}\mathbf{z}\right\rbrace[/latex]subject to [latex]\mathbf{A}\mathbf{z}\geq\mathbf{b}[/latex].library(quadprog)eps = 5e-4y = myocarde[,&#34;PRONO&#34;]*2-1X = as.matrix(cbind(1,myocarde[,1:7]))n = length(y)D = diag(n+7+1)diag(D)[8+0:n] = 0 d = matrix(c(rep(0,7),0,rep(C,n)), nrow=n+7+1)A = Uib = Cisol = solve.QP(D+eps*diag(n+7+1), d, t(A), b, meq=1, factorized=FALSE)qpsol = sol$solution(omega = qpsol[1:7])[1] -0.106642005446 -0.002026198103 -0.022513312261 -0.018958578746 -0.023105767847 -0.018958578746 -1.080638988521(b     = qpsol[n+7+1])[1] 997.6289927Given an observation [latex]\mathbf{x}[/latex], the prediction is[latex display="true"]y=\text{sign}[\mathbf{\omega}^T\mathbf{x}+b][/latex]y_pred = 2*((as.matrix(myocarde0[,1:7])%*%omega+b)&#62;0)-1Observe that here, we do have a classifier, depending if the point lies on the left or on the right (above or below, etc) the separating line (or hyperplane). We do not have a probability, because there is no probabilistic model here. So far.The dual problemThe Lagrangian of the separable problem could be written introducing Lagrange multipliers [latex]\mathbf{\alpha}\in\mathbb{R}^n[/latex], [latex]\mathbf{\alpha}\geq \mathbf{0}[/latex] as[latex display="true"]\mathcal{L}(\mathbf{\omega},b,\mathbf{\alpha})=\frac{1}{2}\|\mathbf{\omega}\|^2-\sum_{i=1}^n \alpha_i\big(y_i(\mathbf{\omega}^T\mathbf{x}_i+b)-1\big)[/latex]Somehow, [latex]\alpha_i[/latex] represents the influence of the observation [latex](y_i,\mathbf{x}_i)[/latex].Consider the Dual Problem, with [latex]\mathbf{G}=[G_{ij}][/latex] and [latex]G_{ij}=y_iy_j\mathbf{x}_j^T\mathbf{x}_i[/latex][latex display="true"]\min_{\mathbf{\alpha}\in\mathbb{R}^n}\left\lbrace\frac{1}{2}\mathbf{\alpha}^T\mathbf{G}\mathbf{\alpha}-\mathbf{1}^T\mathbf{\alpha}\right\rbrace[/latex]subject to [latex]\mathbf{y}^T\mathbf{\alpha}=\mathbf{0}[/latex] and [latex]\mathbf{\alpha}\geq\mathbf{0}[/latex].The Lagrangian of the non-separable problem could be written introducing Lagrange multipliers [latex]\mathbf{\alpha},{\color{red}{\mathbf{\beta}}}\in\mathbb{R}^n[/latex], [latex]\mathbf{\alpha},{\color{red}{\mathbf{\beta}}}\geq \mathbf{0}[/latex], and define the Lagrangian [latex]\mathcal{L}(\mathbf{\omega},b,{\color{red}{\mathbf{\xi}}},\mathbf{\alpha},{\color{red}{\mathbf{\beta}}})[/latex] as[latex display="true"]\frac{1}{2}\|\mathbf{\omega}\|^2+{\color{blue}{C}}\sum_{i=1}^n{\color{red}{\xi_i}}-\sum_{i=1}^n \alpha_i\big(y_i(\mathbf{\omega}^T\mathbf{x}_i+b)-1+{\color{red}{\xi_i}}\big)-\sum_{i=1}^n{\color{red}{\beta_i}}{\color{red}{\xi_i}}[/latex]Somehow, [latex]\alpha_i[/latex] represents the influence of the observation [latex](y_i,\mathbf{x}_i)[/latex].The Dual Problem become with [latex]\mathbf{G}=[G_{ij}][/latex] and [latex]G_{ij}=y_iy_j\mathbf{x}_j^T\mathbf{x}_i[/latex][latex display="true"]\min_{\mathbf{\alpha}\in\mathbb{R}^n}\left\lbrace\frac{1}{2}\mathbf{\alpha}^T\mathbf{G}\mathbf{\alpha}-\mathbf{1}^T\mathbf{\alpha}\right\rbrace[/latex]subject to [latex]\mathbf{y}^T\mathbf{\alpha}=\mathbf{0}[/latex], [latex]\mathbf{\alpha}\geq\mathbf{0}[/latex] and [latex]\mathbf{\alpha}\leq {\color{blue}{C}}[/latex].As previsouly, one can also use quadratic programminglibrary(quadprog)eps = 5e-4y = myocarde[,"PRONO"]*2-1X = as.matrix(cbind(1,myocarde[,1:7]))n = length(y)Q = sapply(1:n, function(i) y[i]*t(X)[,i])D = t(Q)%*%Qd = matrix(1, nrow=n)A = rbind(y,diag(n),-diag(n))C = .5b = c(0,rep(0,n),rep(-C,n))sol = solve.QP(D+eps*diag(n), d, t(A), b, meq=1, factorized=FALSE)qpsol = sol$solutionThe two problems are connected in the sense that for all [latex]\mathbf{x}[/latex][latex display="true"]\mathbf{\omega}^T\mathbf{x}+b = \sum_{i=1}^n \alpha_i y_i (\mathbf{x}^T\mathbf{x}_i)+b[/latex]To recover the solution of the primal problem,[latex display="true"]\mathbf{\omega}=\sum_{i=1}^n \alpha_iy_i \mathbf{x}_i[/latex]thusomega = apply(qpsol*y*X,2,sum)omega                           1                        FRCAR                        INCAR                        INSYS  0.0000000000000002439074265  0.0550138658687635215271960 -0.0920163239049630876653652  0.3609571899422952534486342                        PRDIA                        PAPUL                        PVENT                        REPUL -0.1094017965288692356695677 -0.0485213403643276475207813 -0.0660058643191372279579454  0.0010093656567606212794835while [latex]b=y-\mathbf{\omega}^T\mathbf{x}[/latex] (but actually, one can add the constant vector in the matrix of explanatory variables).More generally, consider the following function (to make sure that [latex]D[/latex] is a definite-positive matrix, we use the nearPD function).svm.fit = function(X, y, C=NULL) { n.samples = nrow(X) n.features = ncol(X) K = matrix(rep(0, n.samples*n.samples), nrow=n.samples) for (i in 1:n.samples){  for (j in 1:n.samples){   K[i,j] = X[i,] %*% X[j,] }} Dmat = outer(y,y) * K Dmat = as.matrix(nearPD(Dmat)$mat)  dvec = rep(1, n.samples) Amat = rbind(y, diag(n.samples), -1*diag(n.samples)) bvec = c(0, rep(0, n.samples), rep(-C, n.samples)) res = solve.QP(Dmat,dvec,t(Amat),bvec=bvec, meq=1) a = res$solution  bomega = apply(a*y*X,2,sum) return(bomega)}On our dataset, we obtainM = as.matrix(myocarde[,1:7])center = function(z) (z-mean(z))/sd(z)for(j in 1:7) M[,j] = center(M[,j])bomega = svm.fit(cbind(1,M),myocarde$PRONO*2-1,C=.5)y_pred = 2*((cbind(1,M)%*%bomega)&#62;0)-1table(obs=myocarde0$PRONO,pred=y_pred)    predobs  -1  1  -1 27  2  1   9 33i.e. 11 misclassification, out of 71 points (which is also what we got with the logistic regression).Kernel Based ApproachIn some cases, it might be difficult to "separate" by a linear separators the two sets of points, like below,It might be difficult, here, because which want to find a straight line in the two dimensional space [latex](x_1,x_2)[/latex]. But maybe, we can distort the space, possible by adding another dimension That's heuristically the idea. Because on the case above, in dimension 3, the set of points is now linearly separable. And the trick to do so is to use a kernel. The difficult task is to find the good one (if any).A positive kernel on [latex]\mathcal{X}[/latex] is a function [latex]K:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}[/latex] symmetric, and such that for any [latex]n[/latex], [latex]\forall\alpha_1,\cdots,\alpha_n[/latex] and [latex]\forall\mathbf{x}_1,\cdots,\mathbf{x}_n[/latex],[latex display="true"]\sum_{i=1}^n\sum_{j=1}^n\alpha_i\alpha_j k(\mathbf{x}_i,\mathbf{x}_j)\geq 0.[/latex]For example, the linear kernel is [latex]k(\mathbf{x}_i,\mathbf{x}_j)=\mathbf{x}_i^T\mathbf{x}_j[/latex]. That's what we've been using here, so far. One can also define the product kernel [latex]k(\mathbf{x}_i,\mathbf{x}_j)=\kappa(\mathbf{x}_i)\cdot\kappa(\mathbf{x}_j)[/latex] where [latex]\kappa[/latex] is some function [latex]\mathcal{X}\rightarrow\mathbb{R}[/latex].Finally, the Gaussian kernel is [latex]k(\mathbf{x}_i,\mathbf{x}_j)=\exp[-\|\mathbf{x}_i-\mathbf{x}_j\|^2][/latex].Since it is a function of [latex]\|\mathbf{x}_i-\mathbf{x}_j\|[/latex], it is also called a radial kernel.linear.kernel = function(x1, x2) { return (x1%*%x2)}svm.fit = function(X, y, FUN=linear.kernel, C=NULL) { n.samples = nrow(X) n.features = ncol(X) K = matrix(rep(0, n.samples*n.samples), nrow=n.samples) for (i in 1:n.samples){  for (j in 1:n.samples){   K[i,j] = FUN(X[i,], X[j,])  } } Dmat = outer(y,y) * K Dmat = as.matrix(nearPD(Dmat)$mat)  dvec = rep(1, n.samples) Amat = rbind(y, diag(n.samples), -1*diag(n.samples)) bvec = c(0, rep(0, n.samples), rep(-C, n.samples)) res = solve.QP(Dmat,dvec,t(Amat),bvec=bvec, meq=1) a = res$solution  bomega = apply(a*y*X,2,sum) return(bomega)}Link to the regressionTo relate this duality optimization problem to OLS, recall that [latex]y=\mathbf{x}^T\mathbf{\omega}+\varepsilon[/latex], so that [latex]\widehat{y}=\mathbf{x}^T\widehat{\mathbf{\omega}}[/latex], where [latex]\widehat{\mathbf{\omega}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}[/latex]But one can also write [latex display="true"]y=\mathbf{x}^T\widehat{\mathbf{\omega}}=\sum_{i=1}^n \widehat{\alpha}_i\cdot \mathbf{x}^T\mathbf{x}_i[/latex]where [latex]\widehat{\mathbf{\alpha}}=\mathbf{X}[\mathbf{X}^T\mathbf{X}]^{-1}\widehat{\mathbf{\omega}}[/latex], or conversely, [latex]\widehat{\mathbf{\omega}}=\mathbf{X}^T\widehat{\mathbf{\alpha}}[/latex].Application (on our small dataset)One can actually use a dedicated R package to run a SVM. To get the linear kernel, uselibrary(kernlab)df0 = dfdf0$y = 2*(df$y=="1")-1SVM1 = ksvm(y ~ x1 + x2, data = df0, C=.5, kernel = "vanilladot" , type="C-svc")Since the dataset is not linearly separable, there will be some mistakes heretable(df0$y,predict(SVM1))         -1 1  -1  2 2  1   1 5The problem with that function is that it cannot be used to get a prediction for other points than those in the sample (and I could neither extract [latex]\omega[/latex] nor [latex]b[/latex] from the 24 slots of that objet). But it's possible by adding a small option in the functionSVM2 = ksvm(y ~ x1 + x2, data = df0, C=.5, kernel = "vanilladot" , prob.model=TRUE, type="C-svc")With that function, we convert the distance as some sort of probability. Someday, I will try to replicate the probabilistic version of SVM, I promise, but today, the goal is just to understand what is done when running the SVM algorithm. To visualize the prediction, usepred_SVM2 = function(x,y){return(predict(SVM2,newdata=data.frame(x1=x,x2=y), type="probabilities")[,2])}plot(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],     cex=1.5,xlab="",     ylab="",xlim=c(0,1),ylim=c(0,1))vu = seq(-.1,1.1,length=251)vv = outer(vu,vu,function(x,y) pred_SVM2(x,y))contour(vu,vu,vv,add=TRUE,lwd=2,nlevels = .5,col="red")Here the cost is [latex]C[/latex]=.5, but of course, we can change itSVM2 = ksvm(y ~ x1 + x2, data = df0, C=2, kernel = "vanilladot" , prob.model=TRUE, type="C-svc")pred_SVM2 = function(x,y){return(predict(SVM2,newdata=data.frame(x1=x,x2=y), type="probabilities")[,2])}plot(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],     cex=1.5,xlab="",     ylab="",xlim=c(0,1),ylim=c(0,1))vu = seq(-.1,1.1,length=251)vv = outer(vu,vu,function(x,y) pred_SVM2(x,y))contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5,col="red")As expected, we have a linear separator. But slightly different. Now, let us consider the "Radial Basis Gaussian kernel"SVM3 = ksvm(y ~ x1 + x2, data = df0, C=2, kernel = "rbfdot" , prob.model=TRUE, type="C-svc")Observe that here, we've been able to separare the white and the black pointstable(df0$y,predict(SVM3))         -1 1  -1  4 0  1   0 6plot(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],     cex=1.5,xlab="",     ylab="",xlim=c(0,1),ylim=c(0,1))vu = seq(-.1,1.1,length=251)vv = outer(vu,vu,function(x,y) pred_SVM3(x,y))contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5,col="red")Now, to be completely honest, if I understand the theory of the algorithm used to compute [latex]\omega[/latex] and [latex]b[/latex] with linear kernel (using quadratic programming), I do not feel confortable with this R function. Especially if you run it several times... you can get (with exactly the same set of parameters)or (to be continued...)