Máquina de Vectores de Soporte (Support Vector Machines, SVM)

Mariano Rivera

25 Sept 2017

Referencias

[1] T. Hastie et al., The elements of Statistical Learning, 2nd Ed. Springer, 2009.

[2] C. Bishop, Pattern Recogniton and Machine Larning, Springer, 2006.

Unas notas previas:
∇x y⊤x=y∇x Ax=A∇x x⊤Ax=A⊤x+Ax∇x x⊤Ax=2Ax,    si  A⊤=A∇x ∥x−y∥2=∇x (x−y)⊤(x−y)=2x⊤(x−y)∇xx2 x⊤Ax=2A,    si  A⊤=A \nabla_x \, y^\top x = y \\ \nabla_x \, A x = A \\ \nabla_x \, x^\top A x = A^\top x + A x \\ \nabla_x \, x^\top A x = 2 Ax, \;\; si \; A^\top=A \\ \nabla_x \, \|x-y\|^2 = \nabla_x \, (x-y)^\top (x-y)= 2 x^\top (x-y) \\ \nabla^2_{xx} \, x^\top A x = 2A, \;\; si \; A^\top=A

Datos para Clasificación Binaria (dos clases).

Sea xi∈Rn\mathbf{x_i} \in \mathbb{R}^n un vector que corresponde al ii-ésimo dato tal que la jj-ésima entrada xij\mathbf{x}_{ij} corresponde a la jj-ésima característica. Luego arreglamos los datos de la forma

(1)
X=[x1⊤x2⊤⋮xm⊤], X = \begin{bmatrix} \mathbf{x}^\top_1 \\ \mathbf{x}^\top_2 \\ \vdots \\ \mathbf{x}^\top_m \\ \end{bmatrix}, \\
de tal que cada dato corresponde a un renglón de la matrix X∈Rm×nX \in \mathbb{R}^{m \times n} y tenemos una etiqueta yi∈{−1,1}\mathbf{y}_i \in \{-1, 1\} asociada a cada dato, que indica la clase a que pertenece. Luego denotamos por

(2)
y=[y1y2⋮ym], y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \\ \end{bmatrix},\\
el vector de etiquetas

SVMs, caso linealmente separable

Considere el conjunto de datos como el mostrado en la figura siguiente, donde el color indica la clase a que pertenece.

Note que podemos trazar una linea recta que defina la frontera entre los dos conjuntos (las clases). Por ejemplo, considere las fronteras propueetas en la figura siguiente. Las tres propuestas son correctas en el sentido de que no se equivocan al dejar datos mal clasificados. Pero aparentemente la linea roja es mejor, o nos parece mejor, la pregunta es ¿porque?

La razón estriba en que sabemos que los datos de entrenamiento son sólo un subconjunto del universo de datos posibles, y además generalmente los datos tienen ruido: si los observamos múltiples veces su posición cambiará ligeramente. Por ello, si dejamos fronteras muy cerca de los datos de entrenamiento podriamos tener problemas con datos que no has sido presentados a la red. Por ello, debemos mantener la frontera lo mas posible alejada de los datos de entrenamiento. Debemos maximizar el margen δ\delta de acuerdo con la siguente figura.

Para ilustrar como calcular el separador lineal de máximo margen, consideremos los datos de la iulstración pasada. Luego les añadimos como entrada adicional su respectiva etiqueta, −1-1 o 11. La gráfica de los puntos (ahora tri-dimensionales) queda como en el panel de la izquierda de la siguiente figura.

Para simplificar aún mas la ilustración, asumamos datos originalmente unidimensinales, como en el panel de la derecha de la figura anterior.

La estrategia es encontrar el plano con menor pendiente tal que los puntos con etiqueta 11 estén por debajo y los puntos con etiqueta −1-1 esten por encima del plano. Es como si tomaramos una varilla sobre nuestra figura anterior que asemeja clavos hacia arriba y hacia abajo, luego la desplazamos y acostamos hasta que se apoye en un clavo de la clase −1-1 y en otro de la clase 11. Vea la siguiente figura.

Matemáticamente esto se escribe como

(3)
arg⁡min⁡w,b12w⊤ws.a                                                                        w⊤xi+b≥ 1,      yi=1                            w⊤xi+b≤−1,    yi=−1 \underset{w,b}{\arg \min} \frac{1}{2} w^\top w \\ \text{s.a} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ \;\;\;\;\;\;\;\;\;\;\;\; w^\top x_i + b \ge \,1, \;\;\; y_i=1 \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; w^\top x_i + b \le -1, \;\; y_i=-1
El problema de Programación Quadrática (Programa Cuadrático) dado por (3) es nuestra primera versión de SVM, ahora sobre las restricciones notamos que si multiplicamos
las restricciones por yiy_i, tenemos:

(4)
yi(w⊤xi+b≥ 1)→yi(w⊤xi+b)−1≥0 y_i(w^\top x_i + b \ge \,1) \rightarrow y_i(w^\top x_i + b) -1 \ge 0

(5)
yi(w⊤xi+b≤ −1)→yi(w⊤xi+b)−1≥0 y_i(w^\top x_i + b \le \,-1) \rightarrow y_i(w^\top x_i + b) -1 \ge 0

Note que, (4) y (5) son iguales por lo que, el problema (3) se puede escribir como

(6)
arg⁡min⁡w,b  12w⊤w                                        s.a            yi(w⊤xi+b−1≥0,      ∀i \underset{w,b}{\arg \min} \;\frac{1}{2} w^\top w \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ \text{s.a} \;\;\;\;\;\; y_i (w^\top x_i + b -1 \ge 0, \;\;\; \forall i
Que es la forma primal de las SVMs.

Una vez que se resuelve el Programa Cuadrático (problema de opmimización con función objetivo cuadrática y restricciónes lineales) obtenemos los parámetros del hiperplano ww y bb separador. Luego para un dato xkx_k que no se haya usado para entrenar (no visto por el entrenamiento) su clase yky_k se calcula con

yk=signo(w⊤xk+b) y_k = signo(w^\top x_k +b)
donde la función signo se define

signo(z)={1,z≥00,otro signo(z) = \left\{ \begin{matrix} 1, & z \ge 0 \\ 0, & \text{otro} \end{matrix} \right.

Forma Dual de las SVMs

Todo problema de optimización con restricciones tienen una forma dual. Es decir, las condiciones KKT son solución para dos problemas: El problema Primal y el problema Dual.

El problema original se conoce como problema primal y su solución requiere resolver las KKT para dos grupos variables: las variables del problema original [ww y bb en el caso (6)] y las variables duales o multiplicadores de Lagrange (aqui las denotaremos por el vector α\alpha).

El problema dual es un problema escrito en términos de las variables duales y cuyos multiplicadores de Lagrange son las variables del problema primal.

Las denominaciones de problemas primal y dual son arbitrarias: generalmente llamamos primal al problema inicialmente planteado, y el dual del dual es el primal.

Para encontrar la forma dual de (6) escribrimos el Lagrangiano:

(7)
min⁡w,bmax⁡α  L(w,b,α)=12w⊤w−∑iαi[yi(w⊤xi+b)−1] \underset{w,b}{\min} \underset{\alpha}{\max} \;\mathcal{L}(w,b, \alpha) = \frac{1}{2} w^\top w - \sum_i \alpha_i[y_i (w^\top x_i + b) -1]

Luego las KKT correspondientes a L\mathcal{L} estan dadas por

(8)
∇wL=0→w−∑iαiyixi=0                             w=∑iαiyixi \nabla_w \mathcal{L}=0 \rightarrow w - \sum_i \alpha_i y_i x_i=0 \\ \;\;\;\;\;\;\;\;\;\;\;\; \;\;\ w = \sum_i \alpha_i y_i x_i\\

(9)
∇bL=0→∑iαiyi=0 \nabla_b \mathcal{L}=0 \rightarrow \sum_i \alpha_i y_i =0

(10)
αi[yi(w⊤xi+b−1]=0 \alpha_i[y_i (w^\top x_i + b -1] =0

(11)
αi≥0 \alpha_i \ge 0

Ahora rescibimos el Lagrangiano (7):

(12)
L(w,b,α)=12w⊤w−∑iαiyiw⊤xi−b∑iαiyi+∑iαi \mathcal{L}(w,b, \alpha) = \frac{1}{2} w^\top w - \sum_i \alpha_i y_i w^\top x_i - b \sum_i \alpha_i y_i + \sum_i \alpha_i

y sustituimos (8) y (9) en (12):

(13)
L(w,b,α)=12∑i∑iαiyixi⊤xjyjαj−∑i∑jαiyixi⊤xjyjαj+∑iαi                  =−12∑i∑iαiyixi⊤xjyjαj+∑iαi \mathcal{L}(w,b, \alpha) = \frac{1}{2} \sum_i \sum_i \alpha_i y_i x_i^\top x_j y_j \alpha_j - \sum_i \sum_j \alpha_i y_i x_i^\top x_j y_j \alpha_j + \sum_i \alpha_i \\ \;\;\;\;\;\;\;\;\; = - \frac{1}{2} \sum_i \sum_i \alpha_i y_i x_i^\top x_j y_j \alpha_j + \sum_i \alpha_i

Ahora cambiamos max⁡αL=min⁡α−L\max_{\alpha} \mathcal{L} = \min_{\alpha} -\mathcal{L} para obtener la forma dual de (6)

(15)
min⁡α12α⊤Qα−1⊤αs.a.                                          y⊤α≥0                        α⪰0 \min_{\alpha} \frac{1}{2} \alpha^\top Q \alpha - {\bf 1}^\top \alpha \\ \text{s.a.} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ \;\;\;\;\;\; y^\top \alpha \ge 0 \\ \;\;\;\;\;\;\;\;\;\;\;\; \alpha \succeq 0 \\
y la matriz Q se calcula como

(16)
Qij=yixi⊤xjyi Q_{ij} = y_i x_i^\top x_j y_i

El problema dual (15) es preferido sobre el primal (6) debido a

Una vez que se resuleve el problema dual para α\alpha, ww se calcula con (8)

(8)
w=∑iαiyixi w = \sum_i \alpha_i y_i x_i
y b se calcula resolviendo las restricciones activas primales :

yi(w⊤xi+b)−1=0,    para    αi>0 y_i(w^\top x_i + b) - 1 = 0, \;\;\text{para} \;\; \alpha_i >0

Basta con una, pero promediar reduce el error, luego

(17)
b=∑i:αi>0(yi−w⊤xi) b = \sum_{i: \alpha_i>0} (y_i - w^\top x_i)

SVMs Robustas

** Clases linealmente separables con dátos atípicos**

Un problema al que nos podemos enfrentar es que los datos de las clases tengan unos cuantos datos que no permitan separar las clases mediante un hiperplano; vea la siguiente figura.

Si pudieramos de alguna manera eliminar esos datos atípico (indicados con marcas rojas en la figura anterior), mediante la SVM podriamos calcular el plano cuya intersección con el espacio de los datos es la línea negra. Luego, visto en una dimensión, el gráfico se veria como sigue.

En el panel de la izquierda se muestran los datos con los datos atípicos (outliers). En el panel central nos hemos concentrado en uno de ellos en particular. El problema que tenemos es que, dada la solución del problema sin outliers, el punto en cuestión (con yi=−1y_i=-1) no satisfacen la resticción de que los datos con etiqueta −1-1 deben estar por arriba del hiperplano (recta en el caso ilustrado) w⊤xi+bw^\top x_i +b.

w⊤xi+b≤−1,    para    yi=−1 w^\top x_i +b \le -1, \;\; \text{para} \;\; y_i =-1

Tenemos dos opciones para hecer que la restricción se cumpla, ambas implican modificarla y conservan el sentido del problema.

  1. Detectar los outliers y multiplicarlos por cero (eliminar esos datos que nos causan problemas). Esto se puede lograr con algo como
    βi[yi(w⊤xi+b)−1]≥0,    ∀yi \beta_i [y_i(w^\top x_i +b) -1] \ge 0, \;\; \forall y_i
    con βi∈{0,1}\beta_i \in \{0,1\} donde esperamos que los datos sean pesados por 11 y sólo por cero en caso de ser outliers. Aun si relajamos βi∈[0,1]\beta_i \in [0,1] (con mayor razón), esta estrategia de solución implica que ww y β\beta se multipliquen lo que nos lleva a restricciones no-lineales.

  2. La otra estrategia es cambiar el resultado evaluar el plano, como si usaramos variables de olgura y las denotamos por rr. Esto es:

(18)
yi(w⊤xi+b)−1+ri≥0,    ∀yi y_i(w^\top x_i +b) -1 + r_i \ge 0, \;\; \forall y_i
por supuesto que implica añadir el costo de usra esta “olguras”, en otro para cualquier plano, siempre podemos encontrar un vector rr que satisfaga (10).

El problema primal de SVM robustas queda entonces

(19)
arg⁡min⁡w,b  12w⊤w+C 1⊤r                                        s.a            yi(w⊤xi+b)+ri−1≥0,      ∀i                                                                r⪰0 \underset{w,b}{\arg \min} \;\frac{1}{2} w^\top w + C\, {\bf 1}^\top r \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ \text{s.a} \;\;\;\;\;\; y_i (w^\top x_i + b) + r_i-1 \ge 0, \;\;\; \forall i \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; r \succeq 0

Para encontrar la forma dual, procedemos como en el caso anterior, determinamos el Lagrangiano:

(20)
min⁡w,b,rmax⁡α,μ  L(w,b,r,α,μ)=12w⊤w+C 1⊤r−∑iαi[yi(w⊤xi+b)+ri−1]−μ⊤r \underset{w,b,r}{\min} \underset{\alpha, \mu}{\max} \;\mathcal{L}(w,b,r, \alpha, \mu) = \frac{1}{2} w^\top w + C\, {\bf 1}^\top r - \sum_i \alpha_i[y_i (w^\top x_i + b) + r_i -1] - \mu^\top r
donde C≥0C \ge 0 es un parámetro de sensibilidad a outliers. Entonces, las KKTs son:

(21)
∇wL=0→w=∑iαiyixi \nabla_w \mathcal{L} =0 \rightarrow w = \sum_i \alpha_i y_i x_i

(22)
∇bL=0→α⊤y=0 \nabla_b \mathcal{L} =0 \rightarrow \alpha^\top y = 0

(23)
∇rL=0→C 1−α−μ=0 \nabla_r \mathcal{L} =0 \rightarrow C \,{\bf 1} - \alpha - \mu =0

(24)
αi[yi(w⊤xi+b)+ri−1]=0 \alpha_i[y_i (w^\top x_i + b) + r_i -1] =0

(25)
αi,μi≥0 \alpha_i,\mu_i \ge 0

Ahora sustituimos (8) en el Lagrangiano (20):

(26)
L(w,b,r,α,μ)=12α⊤Qαj+C 1⊤r−α⊤Qαj−b α⊤y⎵=0 por (22)−α⊤r+1⊤αi−μ⊤r                  =−12α⊤Qαj+1⊤αi+C 1⊤r−α⊤r−μ⊤r=−12α⊤Qαj+1⊤αi+[C 1−α−μ]⊤r⎵=0 por (23) \mathcal{L}(w,b,r, \alpha, \mu) = \frac{1}{2} \alpha ^\top Q \alpha_j + C \, {\bf 1}^\top r - \alpha ^\top Q \alpha_j - b \underbrace{\,\alpha^\top y}_{=0\text{ por (22)}} - \alpha^\top r + {\bf 1}^\top \alpha_i - \mu^\top r \\ \;\;\;\;\;\;\;\;\; = -\frac{1}{2} \alpha ^\top Q \alpha_j + {\bf 1}^\top \alpha_i + C \, {\bf 1}^\top r - \alpha^\top r - \mu^\top r \\ = -\frac{1}{2} \alpha ^\top Q \alpha_j + {\bf 1}^\top \alpha_i + \underbrace{\left[ C \, {\bf 1} - \alpha - \mu \right]^\top r}_{=0 \text{ por (23)}}

Finalmente el problema Dual de (19) queda

(27)
min⁡α12α⊤Qα−1⊤αs.a.                                          y⊤α≥0                        0⪯α⪯C \min_{\alpha} \frac{1}{2} \alpha^\top Q \alpha - {\bf 1}^\top \alpha \\ \text{s.a.} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\ \;\;\;\;\;\; y^\top \alpha \ge 0 \\ \;\;\;\;\;\;\;\;\;\;\;\; 0 \preceq \alpha \preceq C

donde la última restricción vienen de (23): dado que μ⪰0\mu \succeq 0 entonces C−α⪰0C -\alpha \succeq 0, o equivalentemente $\alpha \preceq C $.

</video width=“320” height=“240” controls>

Kernel-SVMs

# Ejemplo de SVM para clasificar Números Escritos

Mariano Rivera

25 Sept 2017

Función para leer la base de datos de Digitos manuscritos MNIST

import os
import struct
import numpy as np

#--------------------------------------------------------------------
def load_mnist (path, kind='train'):
    """
    Load MNIST data from `path`
    """
    labels_path = os.path.join(path,'%s-labels-idx1-ubyte' % kind) 
    images_path = os.path.join(path,'%s-images-idx3-ubyte' % kind)
    with open(labels_path, 'rb') as lbpath:
        magic, n = struct.unpack('>II', lbpath.read(8))
        labels = np.fromfile(lbpath, dtype=np.uint8)
    with open(images_path, 'rb') as imgpath:
        magic, num, rows, cols = struct.unpack(">IIII",imgpath.read(16))
        images = np.fromfile(imgpath, dtype=np.uint8).reshape(len(labels), 784)
    return images, labels 
#--------------------------------------------------------------------

Directorio donde estan los datos (pueded ser relativo al actual)

pathmnist = r'./../../Z Appendix/Datos/mnist/' 
print(os.listdir(pathmnist))

['.DS_Store', 't10k-images-idx3-ubyte', 't10k-images-idx3-ubyte.gz', 't10k-labels-idx1-ubyte', 't10k-labels-idx1-ubyte.gz', 'train-images-idx3-ubyte', 'train-images-idx3-ubyte.gz', 'train-labels-idx1-ubyte', 'train-labels-idx1-ubyte.gz']

Lectura de datos de entrenamiento

X_train, y_train = load_mnist(pathmnist, kind='train')
print('Rows: %d, columns: %d' % (X_train.shape[0], X_train.shape[1]))

Rows: 60000, columns: 784

Lectura de datos de prueba

X_test, y_test = load_mnist(pathmnist, kind='t10k')
print('Rows: %d, columns: %d' % (X_test.shape[0], X_test.shape[1]))

Rows: 10000, columns: 784

Los datos estan en forma de cada digito es un reglon de dimension 756 (28*26 pixeles)

Aqui depliegamos los primerso 700 digitos

import matplotlib.pyplot as plt
import numpy as np
plt.imshow(X_train[:700,:], cmap='gray')
plt.title('Datos')
plt.show()

png

Deplegamos undato, enforma de imagen, pasamos de un reblos a una matriz de 28x28 para
que parezca imagen

Leemos la imagen y su etiqueta

idx = 1601
I = np.reshape((X_train[idx,:]),[28,28], order='C')
label = y_train[idx]
plt.imshow(I, cmap='gray')
plt.title(label)
plt.show()

png

Seleciconamos dos clases

idClass0 =3

idClass1 =8

y se hara la clasificacion binaria

# se seleccionan de entre los primeros elemntos para la clseificacion binaria
# y se asigana las clases 0 y 1

idClass0 =3
idClass1 =8

Idx0 = (y_train==idClass0).astype(int)
Idx1 = (y_train==idClass1).astype(int)

X0 =  X_train[Idx0[:60000]==True,:]
X1 =  X_train[Idx1[:60000]==True,:]

Y0 = np.zeros(X0.shape[0])
Y1 = np.ones (X1.shape[0]).astype(int)

Xtrain = np.concatenate((X0,X1), axis=0)
Ytrain = np.concatenate((Y0,Y1), axis=0).astype(int)

print(Xtrain.shape, Ytrain.shape)
print('número de ejemplos totales = ',  Xtrain.shape[0])
print('número de ejemplos en Clase 0 = ', X0.shape[0])
print('número de ejemplos en Clase 1 = ', X1.shape[0])

(11982, 784) (11982,)
número de ejemplos totales =  11982
número de ejemplos en Clase 0 =  6131
número de ejemplos en Clase 1 =  5851

Paquete con la SVM

probamos solo el caso Lineal, robusto (C=1)

from sklearn import svm 

# crear la SVM y le pasamos los datso de entrenamiento
'''
Cambiamos en la SVM: 
    shrinking= False 
'''

# Clasificador lineal: SVM
clf = svm.LinearSVC(random_state=0)

# Clasificador no-lienal Kernel-SVM
#clf = svm.SVC(kernel='rbf', shrinking= False)

clf.fit(Xtrain, Ytrain)

LinearSVC(C=1.0, class_weight=None, dual=True, fit_intercept=True,
     intercept_scaling=1, loss='squared_hinge', max_iter=1000,
     multi_class='ovr', penalty='l2', random_state=0, tol=0.0001,
     verbose=0)

Generamos las etiquetas predichas para los datos de entrenamiento

Ypred=clf.predict(Xtrain) 
# errores
Err= np.abs(Ytrain-Ypred)
# numero de errores
print('número de errores =', np.sum(Err))

número de errores = 356

Desplegamos una imagens para ver el nivel de ruido

idx = 100
label = int(Ytrain[idx])

plt.subplot(121)
I = np.reshape((Xtrain[idx,:]),[28,28], order='C')

plt.imshow(I, cmap='gray')
label = 'Clase: ' + np.str(label) 
plt.title(label)

plt.subplot(122)
In = np.reshape((Xtest[idx,:]),[28,28], order='C')
plt.imshow(In, cmap='gray')
plt.title(label)

plt.show()

png

Prueba

La validación anterior fue con los datos de Entrenamiento (test) debe probarse con lo de prueba (test)

Idx0test = (y_test==idClass0).astype(int)
Idx1test = (y_test==idClass1).astype(int)

X0test =  X_test[Idx0test[:10000]==True,:]
X1test =  X_test[Idx1test[:10000]==True,:]

Y0test = np.zeros(X0test.shape[0])
Y1test = np.ones (X1test.shape[0]).astype(int)

Xtest = np.concatenate((X0test,X1test), axis=0)
Ytest = np.concatenate((Y0test,Y1test), axis=0).astype(int)

print(Xtest.shape, Ytest.shape)
print('número de Pruebas totales = ',  Xtest.shape[0])
print('número de Pruebas en Clase 0 = ', X0test.shape[0])
print('número de Pruebas en Clase 1 = ', X1test.shape[0])

(1984, 784) (1984,)
número de Pruebas totales =  1984
número de Pruebas en Clase 0 =  1010
número de Pruebas en Clase 1 =  974

Ypred=clf.predict(Xtest) 
# errores
Err= np.abs(Ytest-Ypred)
# numero de errores
print('número de errores =', np.sum(Err))


número de errores = 84

Ejercicios:

  1. Probar clasificación multiclase con SVM
  2. Comparar lineal-SVM vs kernel-SVM

Para ver una implementación de algoritmo de optimización para SVM ver la siguiente liga