Maszyna wektorów nośnych / maszyna wektorów podpierających
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}\)
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}||}\)
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ą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
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")
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\)
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})\]
\[K(x_{i},x_{k}) = (1 + \sum_{j=1}^{p} x_{i,j} x_{k,j})^d \]
\[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
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