这次一定要弄懂-SVM-4-总结Hard Margin SVM的求解过程，并用python实现

4-1 总结Hard Margin SVM的数学推导过程

我们从原问题：
$$\underset{\omega }{\min }\frac{||\omega |{|}^2}{2}$$$$s.t.y^{(i)}({\omega }^T{\mathbf{x}}^{(\mathbf{i})}+b)⩾1\text{i=1,2,...m}$$

转化为对偶问题：
$$\underset{\alpha }{\min }\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}_i{y}_j{x}^{(i)}\cdot x^{(j)}-\sum _{i=1}^m{\alpha }_i$$

$s.t.$
$${\alpha }_i⩾0i=1,2,...,m$$$$\sum _{i=1}^m{\alpha }_i{y}_i=0$$
同时因为原问题满足Slater条件，所以原问题与对偶问题为强对偶关系，最优解是一致的。所以接下来全力以赴解决对偶问题就可以啦。

解决对偶问题的SMO算法：

1. 为${\alpha }_i$设置初始值，让${\alpha }_i=0，i=1,2,...m$

2. 外层循环：
   （1）根据KKT条件选择两个优化变量
   变量一：根据KKT条件选择出第一个违反KKT条件的变量，作为${\alpha }_i$
   $$\{\begin{array}{ll}{\alpha }_i>0& y_i{u}_i=1\\ {\alpha }_i=0& y_i{u}_i⩾1\end{array}$$
   其中 $u_i=\sum _{j=1}^m{y}_j{\alpha }_j{X}_j\cdot X_i+b$

   第二个变量${\alpha }_j$的寻找
   通过在所有不违反KKT条件的乘子中，寻找使$|E_i-E_j|$最大的${\alpha }_j$
   其中$E_i=u_i-y_i$

   （2）求解子问题
   ${\alpha }_j$的值为
   $${\alpha }_j^{new}=\{\begin{array}{rlr}& L& \text{if }{\alpha }_j^{best}<L\\ & {\alpha }_j^{best}& \text{if }L⩽{\alpha }_j^{best}⩽H\\ & H& \text{if }{\alpha }_j^{best}>H\end{array}$$

   其中${\alpha }_j^{best}={\alpha }_j^{\ast }+\frac{y_j(E_i-E_j)}{\eta }$
   $E_i=u_i-y_i$，$\eta =K_{ii}+K_{jj}-2K_{ij}$

   { L = 0 , H = α i + α j y i y j = 1 L = max ⁡ { 0 , α i − α j } , H = + ∞ y i y j = − 1 \begin{cases} L=0,&amp;H=\alpha_i+\alpha_j&amp;\quad y_iy_j=1\\ L=\max\{0,\alpha_i - \alpha_j\},&amp;H=+\infin&amp;\quad y_iy_j=-1 \end{cases} {L=0,L=max{0,αi​−αj​},​H=αi​+αj​H=+∞​yi​yj​=1yi​yj​=−1​
   ${\alpha }_i$的值
   ${\alpha }_i^{new}={\alpha }_i^{\ast }+s({\alpha }_j^{\ast }-{\alpha }_j^{new})$

   更新b的值
   $b_1^{new}=b^{\ast }-E_1+y_1{K}_{11}({\alpha }_1^{\ast }-{\alpha }_1^{new})+y_2{K}_{21}({\alpha }_2^{\ast }-{\alpha }_2^{new})$
   $b_2^{new}=b^{\ast }-E_2+y_1{K}_{12}({\alpha }_1^{\ast }-{\alpha }_1^{new})+y_2{K}_{22}({\alpha }_2^{\ast }-{\alpha }_2^{new})$
   最终$b^{new}$的取值为
   $$b^{new}=\frac{b_1^{new}+b_2^{new}}{2}$$
   对b的更新还不是十分确定，先暂时按这样的方式实现下代码

   更新$E_k$的值
   $$E_k^{new}=\sum _{i=1}^m{y}_i{\alpha }_i{K}_{ik}+b^{new}-y_k$$

3. 判断是否收敛，如果已经收敛，则退出结束循环，否则继续循环

思路大体是清晰了，接下来就试着用代码实现下简单的Hard Margin SVM

4-2 python实现Hard Margin SVM

4-2-1 建立数据集

为了方便验证计算的结果，手动建立数据集

import numpy as np
import matplotlib.pyplot as plt

X=np.array([
            [-1,6],
            [1,5],
            [1,7],
            [3,3],
            [5,4],
            [2,0]])
y=np.array([1,1,1,-1,-1,-1])

plt.scatter(X[y==1,0],X[y==1,1],label='+',color='r')
plt.scatter(X[y==-1,0],X[y==-1,1],label='-',color='b')
plt.legend()

可以很明显看出支撑向量分别为(1,5),(3,3)

4-2-2 初始化

我们需要初始化的量：

• $\alpha$=0
• b=0
• E的值
• m的值

def calE(i):
    return u(i)-y[i]

def u(i):
    return np.sum(X.dot(X[i])*alpha*y)+b

4-2-3 设置循环，寻找第一个不满足kkt条件的$\alpha$

$$\{\begin{array}{ll}{\alpha }_i>0& y_i{u}_i=1\\ {\alpha }_i=0& y_i{u}_i⩾1\end{array}$$

def KKT(i):
    yu=y[i]*u(i)
    if alpha[i]==0:
        return yu >= 1
    else:
        return yu == 1

def search_alpha():
    index_list = np.where(alpha>0)
    non_satisfy = np.where(alpha==0)
    index_list= np.append([index_list],[non_satisfy]) #为了优先搜索alpha>0的项
    for i in index_list:
        if KKT(i):
            continue
        
        E_temp=np.abs(E-E[i]) 
        j=np.argsort(E_temp)[-1]#取-1是因为argsort是升序排列
        return i,j

4-2-4 求解子问题

确定${\alpha }_j$的边界
{ L = 0 , H = α i + α j y i y j = 1 L = max ⁡ { 0 , α i − α j } , H = + inf ⁡ y i y j = − 1 \begin{cases} L=0,&amp;H=\alpha_i+\alpha_j&amp;\quad y_iy_j=1\\ L=\max\{0,\alpha_i - \alpha_j\},&amp;H=+\inf &amp;\quad y_iy_j=-1 \end{cases} {L=0,L=max{0,αi​−αj​},​H=αi​+αj​H=+inf​yi​yj​=1yi​yj​=−1​

def HL(i,j):
    if y[i]*y[j]==1:
        L=0
        H=alpha[i]+alpha[j]
    else:
        L=np.max([0,alpha[i]-alpha[j]])
        H=2147483647 #numpy为int32
    return L,H

确定${\alpha }_j$在无约束情况下的最值点
$${\alpha }_j^{best}={\alpha }_j^{\ast }+\frac{y_j(E_i-E_j)}{\eta }$$
其中$E_i=u_i-y_i$，$\eta =K_{ii}+K_{jj}-2K_{ij}$

def K(i,j):
    return X[i].dot(X[j])

def alpha_unclip(i,j):
    eta = K(i,i)+K(j,j)-2*K(i,j)
    if eta<0:
        print(eta<0)
    return alpha[j]+y[j]*(E[i]-E[j])/eta

确定${\alpha }_j^{new}$和${\alpha }_i$的值
$${\alpha }_j^{new}=\{\begin{array}{rlr}& L& \text{if }{\alpha }_j^{best}<L\\ & {\alpha }_j^{best}& \text{if }L⩽{\alpha }_j^{best}⩽H\\ & H& \text{if }{\alpha }_j^{best}>H\end{array}$$
$${\alpha }_i^{new}={\alpha }_i^{\ast }+s({\alpha }_j^{\ast }-{\alpha }_j^{new})$$

def alpha_new(i,j):
    alpha_new_j=0
    alpha_new_i=0
    
    L,H=HL(i,j)
    unclip=alpha_unclip(i,j)
    
    if unclip<L:
        alpha_new_j=L
    elif L<=unclip<=H:
        alpha_new_j=unclip
    else:
        alpha_new_j=H
    
    alpha_new_i=alpha[i]+y[i]*y[j]*(alpha[j]-alpha_new_j)
    return alpha_new_i,alpha_new_j

更新b
$$b_1^{new}=b^{\ast }-E_1+y_1{K}_{11}({\alpha }_1^{\ast }-{\alpha }_1^{new})+y_2{K}_{21}({\alpha }_2^{\ast }-{\alpha }_2^{new})$$
$$b_2^{new}=b^{\ast }-E_2+y_1{K}_{12}({\alpha }_1^{\ast }-{\alpha }_1^{new})+y_2{K}_{22}({\alpha }_2^{\ast }-{\alpha }_2^{new})$$
最终$b^{new}$的取值为
$$b^{new}=\frac{b_1^{new}+b_2^{new}}{2}$$
对b的更新还不是十分确定，先暂时按这样的方式实现下代码

def find_b_new(i,j,alpha_new_i,alpha_new_j):
    b1_new = b - E[i] + y[i]*K(i,i)*(alpha[i]-alpha_new_i)+y[j]*K(j,i)*(alpha[j]-alpha_new_j)
    b2_new = b - E[j] + y[i]*K(i,j)*(alpha[i]-alpha_new_i)+y[j]*K(j,j)*(alpha[j]-alpha_new_j)
    if alpha_new_i>0:
        return b1_new
    elif alpha_new_j>0:
        return b2_new
    else:
        return (b1_new+b2_new)/2

更新$E_k$的值
$$E_k^{new}=\sum _{i=1}^m{y}_i{\alpha }_i{K}_{ik}+b^{new}-y_k$$

def E_update(i,j):
    E[i]=calE(i)
    E[j]=calE(j)

4-2-5 设置完整的计算函数

1. 初始化
2. 循环不满足KKT条件的alpha，并更新

def ifstop():
    for i in range(len(X)):
        if ~KKT(i):
            return False    
    return True

m = len(X)
global b
b = 0
alpha = np.zeros(m)
E = np.array([calE(i) for i in range(m)])

def fit(X,y):
    
    for k in range(10000):
        i,j=search_alpha()
        alpha_new_i,alpha_new_j = alpha_new(i,j)
        b_new = find_b_new(i,j,alpha_new_i,alpha_new_j)
        
        
        
        alpha[i]=alpha_new_i
        alpha[j]=alpha_new_j
        E_update(i,j)
        b=b_new
       
        
    w=np.sum(X*np.tile(y.reshape(-1,1),(1,X.shape[1]))*np.tile(alpha.reshape(-1,1),(1,X.shape[1])),axis=0)
    return w,b

w,b = fit(X,y)

w

array([-0.48169935,  0.73137255])

b

-2.175163398692811

x_plot=np.linspace(-1,5,1000)
y_plot=(-b-w[0]*x_plot)/w[1]

plt.scatter(X[y==1,0],X[y==1,1])
plt.scatter(X[y==-1,0],X[y==-1,1])
plt.plot(x_plot,y_plot)

看起来好还可以，但实际上和sklearn中的svm结果做对比就会发现。。。

from sklearn.svm import LinearSVC
clf=LinearSVC()
clf.fit(X,y)

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=None, tol=0.0001,
     verbose=0)

x_plot=np.linspace(-1,5,1000)
y_plot=(-b-w[0]*x_plot)/w[1]
y_plot2 = (-clf.intercept_-clf.coef_[0,0]*x_plot)/clf.coef_[0,1]
plt.scatter(X[y==1,0],X[y==1,1])
plt.scatter(X[y==-1,0],X[y==-1,1])

plt.plot(x_plot,y_plot)
plt.plot(x_plot,y_plot2)

肯定有问题的

4-3 总结

虽然表面上依照原问题->对偶问题->求解对偶问题，完成了最终的编码，但肯定哪里是有问题。决定把含有松弛变量和惩罚因子的内容再推导一遍后，也许就清晰了。