【复现SVM】使用梯度下降法代替SMO算法求解SVM对偶问题

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
import torch

%matplotlib inline

生成数据

max_r = 0.4
X = []
y = []
for i in range(20):
    r = max_r * np.random.rand()
    rand_num = np.random.rand()
    X.append([0.5 + r * np.cos(2*np.pi * rand_num), 0.5 + r * np.sin(2*np.pi * rand_num)])
    y.append(1)
for i in range(20):
    r = max_r * np.random.rand()
    rand_num = np.random.rand()
    X.append([-0.5 + r * np.cos(2*np.pi * rand_num), -0.5 + r * np.sin(2*np.pi * rand_num)])
    y.append(-1)

plt.subplots(figsize=(10,5))
for i, [x1, x2] in enumerate(X):
    if y[i] == 1:
        plt.scatter( x1, x2, s=30, c='b', marker='o')
    else:
        plt.scatter( x1, x2, s=30, c='r', marker='x')
plt.xlabel("x1") 
plt.ylabel("x2")
plt.show()

高斯核

$$k({\vec{x}}_i,{\vec{x}}_j)=exp(-\frac{||{\vec{x}}_i-{\vec{x}}_j|{|}^2}{2{\sigma }^2})$$

def gaussKernel(vec1, vec2, _sigma = 0.2):
    vec1 = np.array(vec1).reshape(1, -1)
    vec2 = np.array(vec2).reshape(1, -1)
    return torch.FloatTensor([np.exp(- _sigma * np.linalg.norm(vec1- vec2, ord=2)**2)])

SVM

$$\underset{\alpha }{\max }\sum _{i=1}^m{\alpha }_i-\frac{1}{2}\sum _{i=1}^m\sum _{j=1}^m{\alpha }_i{\alpha }_j{y}_i{y}_{j}K({\vec{x}}_i,{\vec{x}}_j)$$
$$s.t.\sum _{i=1}^m{\alpha }_i{y}_i=0$$
$$0\le {\alpha }_i\le C,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}K({\vec{x}}_i,{\vec{x}}_j)-\sum _{i=1}^m{\alpha }_i$$
$${\alpha }_m=-\frac{{\vec{\alpha }}_{1\cdot \cdot m-1}\cdot {\vec{y}}_{1\cdot \cdot m-1}}{y_m}$$

用梯度下降方式求解 $\vec{\alpha }$

M = len(X)
y = torch.Tensor(np.array(y)).reshape(-1, 1)
# vector alpha
alphaLack = torch.zeros((1, M-1), requires_grad = True)
alphaLast = 0

from IPython.display import clear_output 
epochs = 500
lr = 0.1
K = gaussKernel
C = 10.0

for epoch in range(epochs):
    objFunc = 0.0

    for i in range(M):
        for j in range(M):
            if i != M-1 and j != M-1:
                objFunc += 0.5*alphaLack[0,i] * alphaLack[0,j] * y[i,0] * y[j,0] * K(X[i],X[j])
            elif  i == M-1 and j != M-1: 
                objFunc += 0.5*alphaLast * alphaLack[0,j] * y[i,0] * y[j,0] * K(X[i],X[j])
            elif  i != M-1 and j == M-1: 
                objFunc += 0.5*alphaLack[0,i] * alphaLast * y[i,0] * y[j,0] * K(X[i],X[j])
            else:
                objFunc += 0.5*alphaLast * alphaLast * y[i,0] * y[j,0] * K(X[i],X[j])
    objFunc -= (torch.sum(alphaLack)+ alphaLast)
    objFunc.backward()
#     if epoch%300 == 0:
#         lr = lr/2
    alphaLack.data -= lr*alphaLack.grad
    alphaLack.grad.zero_()
    
    # 软间隔
    for k in range(alphaLack.shape[1]):
        if alphaLack[0,k] > C:
            alphaLack[0,k] = C
        elif alphaLack[0,k] < 0:
            alphaLack[0,k] = 0
    alphaLack = torch.FloatTensor(alphaLack.data).requires_grad_(True)
    
    alphaLast = - (torch.mm(alphaLack, y[0:-1,0].reshape(-1, 1))/y[-1])[0]
    if alphaLast < 0:
        alphaLast = torch.tensor([0])
    elif alphaLast > 10:
        alphaLast = torch.tensor([10])
    
    # 迭代信息
    clear_output(wait=True)
    print("epoch: ", epoch + 1)
    print("objFunc = ", objFunc)
    print("alphaLack = ", alphaLack)
    print("alphaLast = ", alphaLast)

epoch:  500
objFunc =  tensor([-6.3391], grad_fn=<SubBackward0>)
alphaLack =  tensor([[0.0000, 0.0000, 0.0000, 5.8167, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
         0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
         0.2983, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 6.1422,
         0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
         0.0000, 0.0000, 0.0000]], requires_grad=True)
alphaLast =  tensor([0])

合成$\vec{\alpha }$&找出所有的支持向量对应下标

alpha = alphaLack.detach().numpy().reshape(1,-1)
alpha = np.append(alpha[0], alphaLast.detach().numpy().reshape(1,-1))
points = []
for i in range(M):
    if alpha[i] > 0.0001:
        points.append(i)
print("alpha = ",alpha)

alpha =  [0.         0.         0.         5.8166852  0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.29833502 0.         0.         0.         0.         0.
 0.         0.         6.14217091 0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.        ]

求解$b$

$$b=\frac{1}{S}\sum _s(y_{s\in S}-\sum _{s\in S}{\alpha }_i{y}_i{\vec{x}}_i^T{\vec{x}}_s)$$
S = { i ∣ α i > 0 , i = 1 , 2... m } S = \{i|\alpha_i \gt 0, i = 1,2...m\} S={i∣αi​>0,i=1,2...m}

ans = []
b = 0
for p in points:
    ys = y[p, 0]
    sum_record = 0.0
    K = gaussKernel
    for i in points:
        b += alphaLack[0,i] * y[i,0] * K(X[i], X[p])
    b = (1/ys) - b
    ans.append(b.detach().numpy())
b = sum(ans)/len(ans)
print("b = ", b)

b =  [0.04126883]

求解$\vec{w}$

$$\vec{w}=\sum _{i=1}^m{\alpha }_i{y}_i{\vec{x}}_i$$

W = torch.FloatTensor(np.zeros((2,1)))
for i in range(M):
    W = W + alpha[i] * y[i,0] * torch.FloatTensor(np.array(X[i]).reshape(-1, 1))
W = W.detach().numpy()
print("W = ", W)

W =  [[4.8299384]
 [2.977161 ]]

绘制划分方程&标记出支持向量

$$\left[\begin{array}{cc}x1& x2\end{array}\right]\left[\begin{array}{c}w1\\ w2\end{array}\right]+b=0$$

def func(x):
    return - (W[0]*x + b)/W[1]

nums = np.arange(-0.8,0.8,step=.01)
fig, ax = plt.subplots(figsize=(8,7))
ax.plot(nums, func(nums))

for i, [x1, x2] in enumerate(X):
    if y[i] == 1:
        c, marker = 'b', 'o'
    else:
        c, marker = 'r', 'x'
    if alpha[i] > 0.0001:
        c = 'g'
    plt.scatter( x1, x2, s=30, c=c, marker=marker)

plt.xlabel("x1") 
plt.ylabel("x2")
plt.show()