统计学习方法(1) 梯度下降法和SMO算法实现SVM

SVM

SVM是深度学习之前的一种最常用的监督学习方法, 可以用来做分类也可以做回归. 它的本质和感知机类似, 但是额外增加了大间隔的优化目标. 结合前面提到的核技巧, 我们可以把SVM推广到非线性. 这样实现的分类器将有着非常好的性质, 让它一度成为"最优解".

LibSVM

在线性二分类SVM中,我们不止设置一个决策平面$w^{T}x+b=0$,还会有两个支持平面$w^{T}x+b=-1$和$w^{T}x+b=1$. 如果我们假设数据是线性可分的, 那么决策平面一定位于两类样本之间. 而现在我们想从两类样本之间确定一个最好的超平面,这个超平面满足的性质是到两类样本的最小距离最大.很自然的, 我们可以看出这个超平面到两类样本的最小距离应该相等. 这时上面的支持超平面就派上用场了. 设到超平面最近的两个样本是$x_{+}$和$x_{-}$, 样本距离$d=\frac{w^{T}x+b}{||w||}$, 也就是$w^T{x}_{+}+b=-(w^T{x}_{-}+b)=C$. 又因为w的缩放对超平面的性质没有影响, 我们不妨让C = 1, 即两个样本$x_{+}$和$x_{-}$会分别落在$w^{T}x+b=1$和$w^{T}x+b=-1$超平面上. 其他样本都在这两个超平面以外.
如果用上面的一系列假设, 我们的损失函数就应该是: 如果正样本越过了$w^{T}x+b=1$超平面, 即$w^{T}x+b<1$那么它受到随距离线性增长的惩罚. 同样如果正样本越过了$w^{T}x+b=-1$超平面, 即$w^{T}x+b>-1$也受到惩罚. 如果我们设y是样本标签,正样本为+1,负样本为-1. 则我们有保证两类样本到决策超平面的最小距离相等的新型感知机
$$L(w,b)=\frac{ReLU(-y_i(w^T{x}_i+b)+1)}{||w||}$$
同样的道理, 等价于
$$L(w,b)=ReLU(-y_i(w^T{x}_i+b)+1)$$
到此还不算结束, 因为我们还想要这个最小距离最大化. 我们知道$d=\frac{w^{T}x+b}{||w||}$, 而$|w^T x_++b| = |w^T x_-+b| = 1 $,即$d=\frac{1}{||w||}$, 也就是最大化距离, 只需要最小化w即可. 最小化一个参数常用的做法是正则化, 也就是对w施加L2的惩罚(L1的惩罚性质较差, 容易让w为0).这样我们就得到了最大化最小距离的损失.
$$minimize\frac{1}{2}||w|{|}_2^2$$
我们把上面两个损失函数加起来, 并给感知机正确分类的一项乘上一个常数系数, 用来调控(正确分类) vs (最大化间隔)两个目标的重要程度.
$$minimizeL(w,b)=\frac{1}{2}||w|{|}_2^2+C\ast ReLU(-y_i(w^T{x}_i+b)+1)$$
优化这个无约束的优化问题就能得到线性的二分类SVM, 我们用Pytorch的自动求导来实现梯度下降, 优化这个目标函数看看效果.

import torch
import numpy as np
import matplotlib.pyplot as plt
import torch.nn.functional as F
import torch.nn as nn
import matplotlib
import random

torch.manual_seed(2020)
n = 50

Xp = torch.randn(n,2)
Xp[:,0] = (Xp[:,0]+3)
Xp[:,1] = (Xp[:,1]+4)
A = torch.tensor([[-1.,0.8],[2,1]])
Xp = Xp.mm(A)

Xn = torch.randn(n,2)
Xn[:,0] = (Xn[:,0]-1)
Xn[:,1] = (Xn[:,1]-2)
A = torch.tensor([[-0.5,1],[1,0.5]])
Xn = Xn.mm(A)

X = torch.cat((Xp,Xn),axis = 0)
y = torch.cat((torch.ones(n),-1*torch.ones(n)))
y = y.reshape(-1,1)


class LibSVM:
    def __init__(self, C = 1, lr = 0.01):
        '''
        超参数包括松弛变量C和学习率lr
        要学习的参数是一个线性超平面的权重和偏置
        '''
        self.C = C
        self.lr = lr
        self.weights = None
        self.bias = None
    
    def train(self):
        self.weights.requires_grad = True
        self.bias.requires_grad = True
        
    def eval(self):
        self.weights.requires_grad = False
        self.bias.requires_grad = False
    
    def fit(self, X, y, max_iters = 1000):
        '''
        X是数据张量, size(n,m)
        数据维度m, 数据数目n
        y是二分类标签, 只能是1或-1
        '''
        n,m = X.shape
        y = y.reshape(-1,1)
        self.weights = torch.randn(m,1)
        self.bias = torch.randn(1)
        self.train()
        
        for step in range(max_iters):
            out = X.mm(self.weights)+self.bias # 前向计算
            # 损失计算
            loss = 0.5*self.weights.T.mm(self.weights)+\
            self.C*torch.sum(F.relu(-y*out+1))
            # 自动求导
            loss.backward()
            # 梯度下降
            self.weights.data -= self.lr*self.weights.grad.data
            self.bias.data -= self.lr*self.bias.grad.data
            self.weights.grad.data.zero_()
            self.bias.grad.data.zero_()
            
        return loss
    
    def predict(self, x, raw = False):
        self.eval()
        out = x.mm(self.weights)+self.bias
        if raw: return out
        else: return torch.sign(out)
        
    def show_boundary(self, low1, upp1, low2, upp2):
        # meshgrid制造平面上的点
        x_axis = np.linspace(low1,upp1,100)
        y_axis = np.linspace(low2,upp2,100)
        xx, yy = np.meshgrid(x_axis, y_axis)
        data = np.concatenate([xx.reshape(-1,1),yy.reshape(-1,1)],axis = 1)
        data = torch.tensor(data).float()
        # 用模型预测
        out = self.predict(data, raw = True).reshape(100,100)
        out = out.numpy()
        Z = np.zeros((100,100))
        Z[np.where(out>1)] = 1.
        Z[np.where(out<-1)] = -1.
        # 绘制支持边界
        plt.figure(figsize=(8,6))
        colors = ['aquamarine','seashell','palegoldenrod']
        plt.contourf(xx, yy, Z, cmap=matplotlib.colors.ListedColormap(colors))
        # 绘制决策边界
        slope = (-self.weights[0]/self.weights[1]).data.numpy()
        b = (-self.bias/self.weights[1]).data.numpy()
        xx = np.linspace(low1, upp1)
        yy = slope*xx+b
        plt.plot(xx,yy, c = 'black')

model = LibSVM()
model.fit(X,y)
model.show_boundary(-5,12,-5,10)
plt.scatter(X[:n,0].numpy(),X[:n,1].numpy() , marker = '+', c = 'r')
plt.scatter(X[n:,0].numpy(),X[n:,1].numpy() , marker = '_', c = 'r')

Kernel SVM

有时, 我们需要的是超越超平面的决策边界, 它应该能容忍非线性可分的两类样本, 并给出非线性的解. 这一点仅靠上面的方法是无法实现的, 为此我们可以把核技巧带到LibSVM中, 得到二分类的KerbSVM. 复习一下核方法, 我们使用核方法的核心思想是在向量内积得到的标量上进行非线性变换, 从而间接实现高维映射.
那么怎么把内积带到上面的线性模型里呢? 观察到权重w和数据点x的维度相同, 那么如果我们把w写成x的线性组合, 就能得到x和x的内积. 即
$$w=\sum _{i=1}^N{\alpha }_i{x}_i^T$$
$$y=\sum _{i=1}^N{\alpha }_{i}x{x}_i^T$$
在内积上套用核函数就有非线性的广义线性模型, 常用的非线性核函数有高斯核, 多项式核等等
$$k(x,y)=(x{y}^T{)}^d$$
$$k(x,y)=exp(-\sigma ||x-y|{|}^2)$$
y = ∑ i = 1 N α i ϕ ( x ) ϕ ( x i ) T + b = ∑ i = 1 N α i k ( x , x i ) + b y = \sum_{i=1} ^{N} \alpha_i \phi(x) \phi(x_i)^T+b = \sum_{i=1} ^{N} \alpha_i k(x,x_i)+b y=