so as the title suggest I am trying to implement myself SVM method in python, using a polynomial kernel and soft-margin. I'll first discuss my understanding of the algorithmic process to check if there is any misunderstanding, and then i'll discuss about my code and where I think it might fail.
I first implemented the quadratic maximisation problem to solve in the dual formulation:
$$\begin{array}{ll} \text{maximize} & f(c_1,...,c_n) = \displaystyle\sum_{i=1}^n c_i - \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n y_i c_i k(x_i,x_j) y_j c_j\\ \text{subject to} & \displaystyle\sum_{i=1}^n c_i y_i = 0\\ & 0 \leq c_i \leq \frac{1}{2 n \lambda}, \quad \forall i\end{array}$$
For this I used CVXOPT, which need the following formulation :
$$\begin{array}{ll} \text{minimize} & 1/2 x^T P x + q^T x\\ \text{subject to} & Gx \leq h\\ & Ax=b\end{array}$$
I identified everything and got the following matrices, with some help of Math Stack Exchange:
$P_{ij}=y_ik(x_i, x_j)y_j$, $G=\begin{bmatrix} I_n \\ -I_n \end{bmatrix}$, b = [0,...,0] of size 1*n, A is the diagonal matrix of size n*n with the $y_i$ on the diagonal and $h=[\frac{1}{2 \pi n \lambda},\frac{1}{2 \pi n \lambda},...,\frac{1}{2 \pi n \lambda},0,0,....,0]$ of size 1*2n
I believe this part is correct. Thus we obtain the coefficients $c_1,...,c_n$. I'll put all my code in the end.
From my understanding of the wikipedia article, we then have to find some i such that $c_i<\frac{1}{2 \pi n \lambda}$.
After finding this particular value of i, we compute the intercept b using the following formula :
$$ b =[ \sum_{i=1}^n c_j y_j k(x_j,x_i) ] - y_i$$
this conclude the training of the classifier and we can then classify some new input vectors z using the following formula :
$$ z_{class} = sgn([ \sum_{i=1}^n c_i y_i k(x_i,z) ] - b) $$
implementation
For the kernel function, I chose the very simple polynomial kernel and I've trained my classifier on a very simple toy dataset, using the sklearn method make circle. Here is the full code :
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 17 02:53:27 2020
@author: nathan
"""
import numpy as np
import sklearn.datasets
import matplotlib.pyplot as plt
import csv
import math
from cvxopt import matrix
from cvxopt import solvers
import sklearn.svm
#lam is the parameters of the L2 regularization
#sim_mat is the similarity matrix
#y are the labels
#output : returns the vector c optimising the quadratic problem
#for the polynomial kernel with L2 regularisation
def solve_max_quadratic(y,sim_mat,lam):
n=len(sim_mat)
P=np.zeros((n,n))
A=y*np.eye(n)
for i in range(len(P)):
for j in range(i+1):
P[i,j]=A[i,i]*sim_mat[i,j]*A[j,j]
P[j,i]=P[i,j]
A=matrix(A,tc='d')
P=matrix(P,tc='d')
q=matrix(-1*np.ones(n),tc='d')
G=matrix(np.vstack((np.eye(n),-np.eye(n))),tc='d')
h=matrix(np.hstack(((1/(2*n*lam))*np.ones(n),np.zeros(n))),tc='d')
b=matrix(np.zeros(n),tc='d')
sol = solvers.qp(P,q,G,h,A,b)
c=sol['x']
return np.matrix(c)
#return the first indice i such that ci <= 1/(2*n*lam)
def find_in_boundary(c,lam):
n=len(c)
for i in range(n):
if c[i] <= 1/(2*n*lam) :
return i
#if it's not found
return 100000000000
#c is the vector of solution of the quadratic maximisation problem
#y is the vector of labels
#index is the indice computed by find_in_boundary
#sigma is the parameter for the gaussian kernel
#A is the matrix of vectors
def compute_intercept_classifier(A,c,y,index,sigma) :
n=len(c)
opt_var=c
b=0
for i in range(n) :
b+=opt_var[i]*y[i]*poly_kernel(A[i,:],A[index,:],sigma)
return b-y[index]
#predict z with the trained classifier
def predict(z,A,c,y,b,sigma) :
n=len(c)
pred=0
for i in range(n) :
pred+=c[i]*y[i]*poly_kernel(A[i,:],z,sigma)
pred=pred-b
return np.sign(pred)
def euclidean_distance(x_1, x_2):
sum_squared_distance = 0
for i in range(len(x_1)):
sum_squared_distance += math.pow(x_1[i] - x_2[i], 2)
return math.sqrt(sum_squared_distance)
def poly_kernel(xi,xj,sigma) :
return (np.sum(xi*xj)+sigma)**2
def similarity_matrix(A,sigma) :
sim=np.zeros((len(A),len(A)))
for i in range(len(A)):
if(i%10 ==0) :
print("similarity matrix : "+str(i)+" lines computed")
for j in range(i+1):
sim[i,j]=poly_kernel(A[i,:],A[j,:],sigma)
sim[j,i]=sim[i,j]
return sim
np.random.seed(0)
d=2
n_sample=100
data_train, label_train = sklearn.datasets.make_circles(n_samples=n_sample, factor=.3, noise=0.08)
label_train[label_train==0]= -1
plt.figure()
plt.title("Original space")
reds = label_train == -1
blues = label_train == 1
plt.scatter(data_train[reds, 0], data_train[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(data_train[blues, 0], data_train[blues, 1], c="blue",
s=20, edgecolor='k')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
sigma=0.1
sim=similarity_matrix(data_train,sigma)
lam=0.01
c=solve_max_quadratic(label_train,sim,lam)
index=find_in_boundary(c,lam)
b=compute_intercept_classifier(data_train,c,label_train,index,sigma)
y_pred=np.zeros(n_sample)
i=0
for data in data_train :
y_pred[i]=predict(data,data_train,c,label_train,b,sigma)
i=i+1
plt.figure()
plt.title("predicted value")
reds = y_pred == -1
blues = y_pred == 1
plt.scatter(data_train[reds, 0], data_train[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(data_train[blues, 0], data_train[blues, 1], c="blue",
s=20, edgecolor='k')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
The code output two figures, the first one is the plot of the 2D datasets with the train vectors and their labels represented with different colors, the second one are the test vectors (here the same as the train vectors) with the predicted labels. For now, it doesn't work at all.
After some debugging, I believe that the solution returned by the solver of CVXOPT is weird. All coefficients are between $1e^{-25}$ and $1e^{-28}$, hence the intercept will be either very close to 1 or to -1, and every prediction will either be labelled 1 or labelled -1, depending of the particular i chosen.
This tells me that the $c_i$ are the problems somehow and they should be bigger. At this point I don't know what to do, I believe everything else is correct. I also tried on another dataset and I had very similar results for the $c_i$, which reinforce my conviction that they are the problem and therefore that there is a problem with the way I use CVXOPT.
