Оптимизация двойной задачи Вольфа в SVM

Question

Я работаю над машиной опорных векторов (SVM), чтобы написать ее с нуля. Как вы знаете, для двойной формулировки SVM требуется оптимизировать двойную задачу Вольфа относительно множителя Лагранжа следующим образом:

Я много искал, чтобы найти простое описание численного решения этой проблемы, но не смог найти подходящий пример. Наконец, после долгих усилий, я решил эту проблему в Python, используя алгоритм Gradient Descent, и решил поделиться своим решением, чтобы узнать отзывы других пользователей об этом. Пожалуйста, запустите мой код и дайте мне знать, если в нем есть ошибки. Я был бы рад, если вы дадите мне свой отзыв.

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(6)
import math
from sklearn import metrics

from sklearn.datasets.samples_generator import make_blobs
X,y =  make_blobs(n_samples=50,n_features=2,centers=2,cluster_std=1.05,random_state=40)
plt.scatter(X[:,0],X[:,1],marker='o',c=y)
plt.axis([-5,10,-12,-1])
plt.show()

y[y==0]=-1

#all the required variables 
w=[] #weights 2 dimensional vector
b=[] #bias

**# find max and min values of predictor variables (here X) to use it to specify initial values of w and b 

max_feature_value=np.amax(X)
min_feature_value=np.amin(X)**

w_optimum = max_feature_value*0.5
b_optimum = min_feature_value*0.5

w = [w_optimum for i in range(X.shape[1])]   # w shoulb be a vector with dimension of the independent features (here:2)
wt_b=w
#wt_b=[0.34,0.34]



# using the Gradient descent method to find maximum value of Wolfe dual Lagrangian function

l_rate=0.01
n_epoch = 100

# Create matrix of xi*xj
k_xi_xj=np.zeros((X.shape[0] , X.shape[0]))

for i in range (X.shape[0]):
    for j in range(X.shape[0]):
        k_xi_xj[i,j]=np.dot(X[i,:],X[j,:])


#alpha_i_t=[0.001 for i in range(X.shape[0])]
alpha_i_t=np.random.uniform(-10,10,X.shape[0])
alpha_i=alpha_i_t

wolfe_dual_list=[]
wolfe_dual_list_total=[]
wolfe_dual_list_total_ave=[]
for epoch in range(n_epoch):
    for i in range (X.shape[0]):
        wolfe_dual = 0
        for i in range (X.shape[0]):
            for j in range(X.shape[0]):
                wolfe_dual +=((np.sum(alpha_i))-0.5*(alpha_i[i]*alpha_i[j]*y[i]*y[j]*k_xi_xj[i,j]))
        wolfe_dual_list.append(wolfe_dual)
        wolfe_dual_list_total.append(wolfe_dual)

        alpha_i[i] = alpha_i[i] + l_rate* np.dot(alpha_i,y)

        print('>epoch=%d, lrate=%.3f,wolfe_dual=%.3f' % (epoch, l_rate,wolfe_dual))
    wolfe_dual_list_total_ave.append(np.mean(wolfe_dual_list_total))
    wolfe_dual_list_total=[]

# plot wolfe dual lagrangian function  values
fig,ax = plt.subplots(1, 1,constrained_layout=True,figsize=(8, 4))
ax.plot(wolfe_dual_list)
ax.set_xlabel('iterate',fontsize=15)
ax.set_ylabel('wolfe dual values',fontsize=15)
#plt.grid(True)
plt.show()

# plot average wolfe dual lagrangian function values in each epoch
fig,ax = plt.subplots(1, 1,constrained_layout=True,figsize=(8, 4))
ax.plot(wolfe_dual_list_total_ave)
ax.set_xlabel('iterate',fontsize=15)
ax.set_ylabel('wolfe dual values',fontsize=15)
#plt.grid(True)
plt.show()

# Compute w and b

# Because Lagrange multipliers for non-support vectors are almost zero, we can also compute w using only support vectors
# data and their multipliers

# find all vectors with positive multipliers (Support vectors have positive multipliers)

support_vectors=X[alpha_i>1e-7]
support_vectors_y=y[alpha_i>1e-7]
sv_alpha_i=alpha_i[alpha_i>1e-7]

#Compute w
wt_b=0
for i in range (support_vectors.shape[0]):
    wt_b=wt_b+sv_alpha_i[i]*support_vectors_y[i]*support_vectors[i,:]

#Compute b
sum=0
for i in range (support_vectors.shape[0]):
    sum=sum+support_vectors_y[i]-np.dot(wt_b,support_vectors[i,:])

b=sum/len(support_vectors)


# find two best support vectors to draw them in plot
positive_instances=[]
negative_instances=[]

for i in range(X.shape[0]):
    num=(np.dot(wt_b,X[i,:]))+b
    if  ((y[i]>=1) and (num>=1)):
        positive_instances.append([num,X[i,:]])
    elif ((y[i]<=-1) and (num<=-1)):
        negative_instances.append([num,X[i,:]]) 

sort_positive=sorted([n for n in positive_instances])
sort_negative=sorted([n for n in negative_instances])

positive_support_vector=sort_positive[0][1]
negative_support_vector=sort_negative[-1][1]

model_support_vectors=np.stack((positive_support_vector,negative_support_vector),axis=-1)

#wt_b=np.asarray([0.32,0.32])
#b=1.49

# visualize the data-set

colors = {1:'r',-1:'b'}
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

plt.scatter(X[:,0],X[:,1],marker='o',c=y)

# plot support vectors
ax.scatter(model_support_vectors[0, :],model_support_vectors[1, :],s=200, linewidth=1,facecolors='none', edgecolors='b')

# hyperplane = x.w+b
# v = x.w+b
# psv = 1
# nsv = -1
# dec = 0
def hyperplane_value(x,w,b,v):
    return (-w[0]*x-b+v) / w[1]

datarange = (min_feature_value*0.9,max_feature_value*1.)
hyp_x_min = datarange[0]
hyp_x_max = datarange[1]

# (w.x+b) = 1
# positive support vector hyperplane
psv1 = hyperplane_value(hyp_x_min, wt_b, b, 1)
psv2 = hyperplane_value(hyp_x_max, wt_b, b, 1)
ax.plot([hyp_x_min,hyp_x_max],[psv1,psv2], 'k')

# (w.x+b) = -1
# negative support vector hyperplane
nsv1 = hyperplane_value(hyp_x_min, wt_b, b, -1)
nsv2 = hyperplane_value(hyp_x_max, wt_b, b, -1)
ax.plot([hyp_x_min,hyp_x_max],[nsv1,nsv2], 'k')

# (w.x+b) = 0
# positive support vector hyperplane
db1 = hyperplane_value(hyp_x_min, wt_b, b, 0)
db2 = hyperplane_value(hyp_x_max, wt_b, b, 0)
ax.plot([hyp_x_min,hyp_x_max],[db1,db2], 'y--')

#plt.axis([-5,10,-12,-1])
plt.show()

# Predict using the optimized values of w and b    
y_predict=[]
for i in range (X.shape[0]):
    y_hat=np.sign(np.dot(wt_b,X[i,:])+b)
    y_predict.append(y_hat)


from sklearn import metrics
metrics.accuracy_score(y_predict, y)