Suppport Vector Machines

Maszyna wektorów nośnych / maszyna wektorów podpierających

Podstawowe pojęcia z geometrii analitycznej

Odległość punktu \(M=(x_0, y_0)\) od prostej \(l\) danej równaniem w postaci ogólnej:

\(l: Ax + By + C = 0\)

wyraża się wzorem:

\(d(M,l) = \frac{|Ax_0 + By_0 + C|}{\sqrt{A^2+B^2}}\)

Moduł wektora (długość wektora, norma):

\(\mathbf{u}=[u_1, u_2]\) \(u=\sqrt{u_1^2 + u_2^2}\)

ogólnie:

\(||\mathbf{u}||=\sqrt{u_1^2 + u_2^2 +\cdots + u_p^2}\)

Hierpłaszczyzna rozdzielająca

Rozważmy zadanie wyznaczenia “marginesów” ograniczonych dwiema równoległymi hiperpłaszczyznami, we wnętrzu których nie leży ani jeden element próby uczącej:

\(x_1, \cdots, x_n \in R^p\)

z przypisanymi etykietami klas:

\(y_1, \cdots, y_n \in \{-1,1\}\)

Zdefiniujmy parametry hiperpłaszczyzny: \(\mathbf{w}=[w_1, w_2, \cdots, w_p]\)

\(\mathbf{w}^T \mathbf{x} +b \geq 1\) gdy \(y_i=1\)

oraz

\(\mathbf{w}^T\mathbf{x} +b \leq -1\) gdy \(y_i=-1\)

\(H_1\) oraz \(H_2\) - hiperpłaszczyzny:

\(H_1: \mathbf{w}^T\mathbf{x} +b = 1\)

\(H_2: \mathbf{w}^T\mathbf{x} +b = -1\)

Punkty na hiperpłaszczyznach \(H_1\) oraz \(H_2\) są końcami wektorów podpierających.

Płaszczyzna rozdzielająca \(H_0\) leżąca pomiędzy:

\(H_0: \mathbf{w}^T\mathbf{x} +b = 0\)

Odległości:

\(d_+\) najmniejsza odległość od hiperpłaszczyzny rozdzielającej \(H_0\) do punktu Pozytywnego (klasa y=1)

\(d_-\) najmniejsza odległość od hiperpłaszczyzny rozdzielającej \(H_0\) do punktu Negatywnego (klasa y=-1)

Szerokość marginesu \(d= d_+ + d_-\)

Odległość hiperpłaszczyzn \(H_0\) oraz \(H_1\):

\(\frac{|\mathbf{w}^T\mathbf{x} +b|}{||\mathbf{w}||}=\frac{1}{||\mathbf{w}||}\)

a całkowita odległość:

\(d= \frac{2}{||\mathbf{w}||}\)

Zadanie optymalizacji

Aby zmaksymalizować szerokość marginesu, należy zminimalizować \(||\mathbf{w}|| \to min\)

przy ograniczeniach:

\(\mathbf{w}^T \mathbf{x} +b \geq 1\) gdy \(y_i=1\)

\(\mathbf{w}^T\mathbf{x} +b \leq -1\) gdy \(y_i=-1\),

które można zapisać :

\(y_i ( \mathbf{w}^T {x_i} +b) \geq 1\) dla \(i=1 \cdots n\)

Hiperpłaszczyzna rozdzielająca: które punkty należy uwzględnić ? - wszystkie (regresja liniowa, sieci neuronowe) - tylko “kłopotliwe”, leżące blisko granicy decyzyjnej : SVM

W rzeczywistości minimalizujemy: \(\frac{||\mathbf{w}||^2}{2}\)

Rozwiązania zadania optymalizacji z wykorzystaniem mnożników Lagrangea’

Rozwiązanie zadania z wykorzystaniem metody mnożników Lagrange’a:

\(L=L(\mathbf{w}, w_0; \boldsymbol{\alpha})=\frac{1}{2}\mathbf{w}'\mathbf{w}-\sum_{j-1}^{n}\alpha_{j}[y_j(\mathbf{w}'\mathbf{x_{j}}+w_{0})-1]\)

gdzie

\(\boldsymbol{\alpha}=(\alpha_{1}, \alpha_{1}, \dots , \alpha_{n})' \ge 0\)

Różniczkując funkcję Lagrange’a kolejno względem \(w\) oraz \(w_{0}\) oraz przyrównując pochodne cząstkowe do zera, otrzymujemy odpowiednio równości:

\(\mathbf{w}=\sum_{j=1}^{n} \alpha_{j}y_{j}\mathbf{x_{j}}\)

\(\sum_{j=1}^{n}\alpha_{j}y_{j} = 0\) , \(\boldsymbol{\alpha} \ge 0\)

Podstawiając do równania powyżej mamy:

\(L_{D}(\boldsymbol{\alpha})=\sum_{j=1}^{n}\alpha_{j}-\frac{1}{2}\sum_{j=1}^{n}\sum_{k=1}^{n} \alpha_{j} \alpha_{k} y_{j} y_{k} \mathbf{x_{j}}'\mathbf{x_{k}}\)

Formułując problem w postaci dualnej zwiększamy wymiar przestrzeni parametrów z \(p+1\) (wektor wag oraz stała \(w_{0}\) ) do \(n\) (mnożniki Lagrange’a ), tj. liczebności próby. Zaletą formy dualnej jest to, że wyraża kryterium optymalizacji w terminach iloczynów skalarnych obserwacji \(x_{j}\).

Forma dualna w notacji macierzowej:

\(L_{D}(\boldsymbol{\alpha})=\boldsymbol{\alpha}'\mathbf{1}-\frac{1}{2}\boldsymbol{\alpha}'\mathbf{D}\boldsymbol{\alpha}\)

gdzie \(\mathbf{1}\) jest n-wymiarowym wektorem jednostkowym,

\(\mathbf{D}\) macierzą symetryczną o elementach:

\(d_{j,k}=y_j y_k \mathbf{x_j}'\mathbf{x_k}\) , \(j,k=1,2,\dots, n\)

Szukamy punktu \(\boldsymbol{\alpha}\) dla którego funkcja \(L_D\) osiąga maksimum przy warunku pobocznym:

\(\boldsymbol{\alpha}'\mathbf{y} = 0\)

gdzie

\(\mathbf{y} = (y_1, y_2, \dots , y_n)'\)

Zgodnie z warunkami Kuhna-Tuckera rozwiązanie powyższego problemu związane jest z dodatkową zależnością postaci:

\(\alpha_{j}[y_j(\mathbf{w}'\mathbf{x_{j}}+w_{0})-1]=0\) dla \(j=1,2,\dots, n\)

Wynika z niej następująca implikacja:

\(\alpha_{j} > 0 \Rightarrow y_j(\mathbf{w}'\mathbf{x_j}+w_{0}) = 1\)

Zatem wszystkie obserwacje, którym odpowiadają niezerowe mnożniki Lagrange’a leżą na hiperpłaszczyznach kanonicznych. Obserwacje te nazywamy wektorami nośnymi.

\(\mathbf{w} = \sum_{ \mathbf{x_{j}} \in SV} \alpha_{j} y_j \mathbf{x_{j}}\)

gdzie SV oznacza zbiór wszystkich wektorów nośnych.,

O ostatecznej postaci funkcji dyksryminacyjnej g decydują wyłącznie wektory nośne. Im większa wartość \(\alpha_{j}\) wektora nośnego \(\mathbf{x_{j}}\), tym większy jego wpływ na kształt granic decyzyjnych. Usunięcie z próby którejkolwiek (a nawet wszystkich) z pozostałych obserwacji nie wpłynie na postać hiperpłaszczyzny .

library(ggplot2)
library(ggrepel)


x1<-runif(50)*20
x2<-runif(50)*20




y<- ifelse( (x1 + 2*x2-30 > 0 ), 1, -1)
beta<-c(1,2,-30)
y<- as.factor(ifelse( (x1 + 2*x2-30 > 0 ), 1, -1))
df<-data.frame(x1,x2,y)
df$id<- as.numeric(row.names(df))
p<- ggplot(df, aes(x=x1))+geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  +scale_shape_manual(values=c(1,2)) + geom_abline(intercept=15, slope=-0.5) + geom_text_repel(aes(y=x2,label=id,color=as.factor(y)),size=3)
p

Klasy separowalne liniowo

library(e1071)


m <- svm(y~x1+x2, df, kernel="linear",cost=100000, type="nu-classification", nu=0.04)
summary(m)
## 
## Call:
## svm(formula = y ~ x1 + x2, data = df, kernel = "linear", cost = 1e+05, 
##     type = "nu-classification", nu = 0.04)
## 
## 
## Parameters:
##    SVM-Type:  nu-classification 
##  SVM-Kernel:  linear 
##       gamma:  0.5 
##          nu:  0.04 
## 
## Number of Support Vectors:  4
## 
##  ( 2 2 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1
#indeksy punktow (wektory podpierajace)
#index of resulting support vectors in the data matrix 
m$index
## [1]  1 42 23 24
#negative intercept 
m$rho
## [1] 0.4764308
#the corresponding coefficients times the training labels
m$coefs
##               [,1]
## [1,]  4.094287e+01
## [2,]  2.500162e+01
## [3,] -6.594449e+01
## [4,] -1.464262e-14
svm<-df[m$index,]
svm_pos<-svm[svm$y==1,]
svm_pos
##           x1        x2 y id
## 1  14.632142  7.989934 1  1
## 42  3.568406 13.640343 1 42
svm_neg<-svm[svm$y==-1,]
head(svm_neg)
##          x1       x2  y id
## 23 10.08174 9.579436 -1 23
## 24 10.41537 8.517665 -1 24
## rownanie pierwszej prostej 

if (nrow(svm_pos) > nrow(svm_neg)) {
  idxs<-as.vector(svm_pos[,4])
  idx<-svm_neg[1,4]
} else {
  idxs<-as.vector(svm_neg[,4])
  idx<-svm_pos[1,4]
}

regr<-lm(x2~x1, df[idxs,])
coef<-coefficients(regr)
beta0 <-coef[1]
beta1 <-coef[2] 


ggplot(df, aes(x=x1)) + 
   geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  +
   scale_shape_manual(values=c(1,2)) +
   geom_abline(intercept=beta0, slope=beta1)  +
   geom_text_repel(aes(y=x2,label=id,color=as.factor(y)),size=3)

#rownanie drugiej prostej 

df[idx,]
##         x1       x2 y id
## 1 14.63214 7.989934 1  1
beta02 <- df[idx,2] - beta1*df[idx, 1]


beta0
## (Intercept) 
##    41.66519
beta02
##      x1 
## 54.5576
beta1
##        x1 
## -3.182559
ggplot(df, aes(x=x1))+geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  + scale_shape_manual(values=c(1,2)) +  geom_abline(intercept=beta02, slope=beta1, color="#00BFC4") +  geom_text_repel(aes(y=x2,label=id,color=as.factor(y)),size=3) +  geom_abline(intercept=beta0, slope=beta1, color="#F8766D")

Klasy liniowo nieseparowalne

Jeżeli nie ma płaszczyzny rozdzielającej, przyjmujemy, że część punktów może się znaleźć “po drugiej stronie”, jednak wiąże się to z pewną karą \(\xi\) przypisaną dla takiego punktu. Do funkcji kryterium wprowadzona jest wówczas funkcja kary:

\(\frac{||\mathbf{w}||^2}{2} + C*\sum_{i=1}^{n}\xi_i\)

a ograniczenia przyjmują postać:

\(y_i ( \mathbf{w}^T {x_i}+b) \geq 1 -\xi_i\) dla \(i=1 \cdots n\)

\(\xi_i \geq 0\) dla \(i=1 \cdots n\)

Nieliniowe granice klas

Klasyfikator z liniowym rozgraniczeniem klas może być przedstawiony w postaci iloczynu skalarnego:

\(f(x)=\beta_{0} + \sum_{i=1}^{n}\alpha_{i} \langle x, x_{i} \rangle\) , gdzie

\(\langle x_{i}, x_{k} \rangle\) oznacza iloczyn skalarny wektorów:

\[\langle x_{i}, x_{k} \rangle = \sum_{j=1}^{p} x_{i,j}x_{k,i}\] Współczynniki \(\alpha_{i}\) przyjmują niezerowe wartości tylko dla wektorów nośnych , a więc jeżeli przez \(\mathbf{SV}\) oznaczymy zbiór indeksów odpowiadających wektorom nośnym to:

\[f(x)=\beta_{0} + \sum_{i \in \mathbf{SV} }\alpha_{i} \langle x, x_{i} \rangle\] Jeżeli iloczyn skalarny dwóch wektorów uogólnimy do funkcji K : \(K(x_{i}, x_{k})\) , zwanej jądrem (kernel), która określa podobieństwo dwóch obserwacji:

\[f(x)=\beta_{0} + \sum_{i \in \mathbf{SV} } \alpha_{i} K( x, x_{i})\]

Kernels

Polynomial kernel

\[K(x_{i},x_{k}) = (1 + \sum_{j=1}^{p} x_{i,j} x_{k,j})^d \]

Radial kernel

\[K(x_{i},x_{k}) = exp(-\gamma \sum_{j=1}^{p}(x_{i,j} – x_{k,j})^2 )\]

require(e1071)
require(ElemStatLearn) 
## Loading required package: ElemStatLearn
library(ggplot2)
#Loading the data
attach(mixture.example) 
## The following object is masked _by_ .GlobalEnv:
## 
##     y
dataset<-mixture.example

df<-data.frame("x1"=dataset$x[,1], "x2"=dataset$x[,2], "y"=dataset$y)

head(df)
##            x1         x2 y
## 1  2.52609297  0.3210504 0
## 2  0.36695447  0.0314621 0
## 3  0.76821908  0.7174862 0
## 4  0.69343568  0.7771940 0
## 5 -0.01983662  0.8672537 0
## 6  2.19654493 -1.0230141 0
p<- ggplot(df, aes(x=x1))+geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  +scale_shape_manual(values=c(1,2)) + geom_abline(intercept=15, slope=-0.5) 
p

#
radial.svm=svm(factor(y) ~ .,data=df,kernel="radial",cost=100,scale=F)

support.vectors<-data.frame(radial.svm$SV)
#dodajemy identyfikator klasy 
support.vectors<-merge(support.vectors, df, by=c("x1", "x2"))
head(support.vectors)
##            x1        x2 y
## 1 -0.01983662 0.8672537 0
## 2 -0.07407910 0.9197022 0
## 3 -0.07711446 0.5051217 0
## 4 -0.19127236 2.4936857 0
## 5 -0.19624633 0.5514036 1
## 6 -0.20068316 2.2815850 1
#wygenerujemy siatkę punktów  
xgrid=expand.grid(x1=px1,x2=px2) #generating grid points
ygrid=predict(radial.svm,newdata = xgrid) #ygird consisting of predicted Response values


#Wykres z zaznaczonymi wektorami podpierajacymi: 
df1<-data.frame("x1"=xgrid$x1, "x2"=xgrid$x2, "y"=ygrid)

p<- ggplot(df, aes(x=x1))+geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  

p<- p +   geom_point( data=df1, aes(x=x1, y=x2,  color=as.factor(y)),  shape=19  , alpha=0.1)

p<- p +   geom_point( data=support.vectors, aes(x=x1, y=x2 ,color=as.factor(y) ),  shape=3  ,  size = 3, alpha=1)

p

Wybór najlepszych parametrów dla modelu

require(e1071)
svm_tune <- tune("svm",  data=df, train.x=as.formula("factor(y) ~ ."), kernel="radial", ranges=list(cost=10^(-2:2), gamma=2^(-2:2)))

svm_tune
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost gamma
##     1     4
## 
## - best performance: 0.15
ygrid=predict(svm_tune$best.model,newdata = xgrid) #ygird consisting of predicted Response values


#Wykres z zaznaczonymi wektorami podpierajacymi: 
df1<-data.frame("x1"=xgrid$x1, "x2"=xgrid$x2, "y"=ygrid)
p<- ggplot(df, aes(x=x1))+geom_point(aes( y=x2,shape=as.factor(y) ,color=as.factor(y)))  
p<- p +   geom_point( data=df1, aes(x=x1, y=x2,  color=as.factor(y)),  shape=19  , alpha=0.1)
p<- p +   geom_point( data=support.vectors, aes(x=x1, y=x2 ,color=as.factor(y) ),  shape=3  ,  size = 3, alpha=1)
p