支持向量机

⽀持向量机（Support Vector Machines，SVM）

• 优点：泛化错误率低，计算开销不⼤，结果易解释。
• 缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适⽤于处理⼆类问题。
  适⽤数据类型：数值型和标称型数据。

间隔和支持向量

假设数据是可分的，该算法就是为了找到一个超平面将两者分开
⽀持向量（support vector）就是离分隔超平⾯最近的那些点。接下来要试着最⼤化⽀持向量到分隔⾯的距离，需要找到此问题的优化求解⽅法

分割超平面的方程为$$w^{T}x+b=0$$
计算点A到分割超平面的距离$$r=\frac{|w^{T}x+b|}{||w||}$$

为了方便计算，类别标签label使用+1和-1：$\mathrm{label}\ast w^{T}x+b$可以同号。
现在的目标就是找出分类器定义中的w和b。为此要找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。对该间隔最大化：$$\arg \underset{w,b}{\max }\left\{\underset{n}{\min }(\mathrm{label}\cdot (w^{\mathrm{T}}x+b))\cdot \frac{1}{\Vert w\Vert }\right\}$$

为了简化问题，令支持向量到超平面距离为1： $r\ge 1$ ，则有 $$\{\begin{array}{ll}{\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_{\mathbf{i}}+b⩾+1,& y_i=+1\\ {\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_{\mathbf{i}}+b⩽-1,& y_i=-1\end{array}$$

两个一类支持向量到超平面的距离之和为 $$\gamma =\frac{2}{||\mathbf{w}||}$$

要使得间隔最大$\max \frac{2}{||w||}$，等价于$\min \frac{1}{2}||w|{|}^2$，即$$\min \frac{1}{2}||w|{|}^{2}s.t.y_i(w^T{x}_i+b)\ge 1$$

对偶问题

问题推导

回顾一下高等数学中约束极值的求法：
$z=f(x,y)$在条件$\psi (x,y)=0$下取得极值，求极值：

1. 构造拉格朗日函数$F(x,y,\lambda )=f(x,y)+\lambda \psi (x,y)$
2. 令$\{\begin{array}{rl}& f_x^{\prime }(x,y)+\lambda {\psi }_x^{\prime }(x,y)=0\\ & f_y^{\prime }(x,y)+\lambda {\psi }_y^{\prime }(x,y)=0\\ & \psi (x,y)=0\end{array}$
3. 将求出的点带入

现问题：求$\frac{1}{2}||w|{|}^2$在$y_i(w^T{x}_i+b)\ge 1$条件下的最小值

解决方法：引入拉格朗日乘子${\alpha }_i⩾0$，构造拉格朗日函数$$L(\mathbf{w},b,\mathbf{\alpha })=\frac{1}{2}\Vert \mathbf{w}{\Vert }^2+\sum _{i=1}^m{\alpha }_i\left(1-y_i({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i+b)\right)$$

其中$\mathbf{\alpha }=({\alpha }_1;{\alpha }_2;\dots ;{\alpha }_m)$ ，对w和b偏导为0可得：$$\begin{array}{rl}\mathbf{w}& =\sum _{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i,\\ 0& =\sum _{i=1}^m{\alpha }_i{y}_i.\end{array}$$

联立后，可解得对偶问题$$\begin{gathered} \underset{\mathbf{\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{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{x}}_j \\ \begin{array}{rl}s.t.& \sum _{i=1}^m{\alpha }_i{y}_i=0,\\ & {\alpha }_i⩾0,i=1,2,\dots ,m.\end{array} \end{gathered}$$

最后模型为$$\begin{gathered} \begin{array}{rl}f(\mathbf{x})& ={\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b\\ & =\sum _{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i^{\mathrm{T}}\mathbf{x}+b.\end{array} \\ s.t.\{\begin{array}{l}{\alpha }_i⩾0;\\ y_{i}f({\mathbf{x}}_i)-1⩾0;\\ {\alpha }_i\left(y_{i}f({\mathbf{x}}_i)-1\right)=0.\end{array} \end{gathered}$$

SMO

推导的式子已经得到了，但是如何求解这些参数？SMO是其中一个高效算法的代表。
SMO的基本思路是先固定${\alpha }_i$之外的所有参数,然后求${\alpha }_i$上的极值.由于存在约束${\sum }_{i=1}^m{\alpha }_i{y}_i=0$,若固定${\alpha }_i$之外的其他变量,则${\alpha }_i$可由其他变量导出.于是,SMO每次选择两个变量${\alpha }_i$和 ${\alpha }_j$,并固定其他参数.这样,在参数初始化后,SMO不断执行如下两个步骤直至收敛:

• 选取一对需更新的变量${\alpha }_i$和 ${\alpha }_j$;
• 固定${\alpha }_i$和 ${\alpha }_j$以外的参数,求解式获得更新后的${\alpha }_i$和${\alpha }_j$.

则约束变为$$\begin{array}{rl}{\alpha }_i{y}_i+{\alpha }_j{y}_j& =c,{\alpha }_i⩾0,{\alpha }_j⩾0\\ & \\ c& =-\sum _{k\ne i,j}{\alpha }_k{y}_k\end{array}$$

确定偏移项b的时候，现实任务中常采用一种更鲁棒的做法：使用所有支持向量求解的平均值（支持向量下标集 S = { i ∣ α i > 0 , i = 1 , 2 , … , m } S=\{i\mid\alpha_{i}>0,i=1,2,\ldots,m\} S={i∣αi​>0,i=1,2,…,m}）：$$b=\frac{1}{|S|}\sum _{s\in S}\left(y_s-\sum _{i\in S}{\alpha }_i{y}_i{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{x}}_s\right)$$

核函数

有时数据并不是线性可分，为解决此问题，可将样本从原始空间映射到一个更高维的特征空间，使得样本在这个特征空间线性可分。

令$\varphi (\mathbf{x})$表示将$\mathbf{x}$映射后的特征向量：$$f(x)={\mathbf{w}}^{\mathrm{T}}\varphi (\mathbf{x})+b$$

对偶问题是$$\begin{gathered} \underset{\mathbf{\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\varphi ({\mathbf{x}}_i{)}^{\mathrm{T}}\varphi ({\mathbf{x}}_j) \\ \begin{array}{rl}s.t.& \sum _{i=1}^m{\alpha }_i{y}_i=0,\\ & {\alpha }_i⩾0,i=1,2,\dots ,m.\end{array} \end{gathered}$$

为了简化$\varphi ({\mathbf{x}}_i{)}^{\mathrm{T}}\varphi ({\mathbf{x}}_j)$，设想一个函数：$$\kappa ({\mathbf{x}}_i,{\mathbf{x}}_j)=⟨\varphi ({\mathbf{x}}_i),\varphi ({\mathbf{x}}_j)⟩=\varphi ({\mathbf{x}}_i{)}^{\mathrm{T}}\varphi ({\mathbf{x}}_j)$$

对偶问题重写为：$$\begin{gathered} \underset{\mathbf{\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\kappa ({\mathbf{x}}_i,{\mathbf{x}}_j) \\ \begin{array}{rl}s.t.& \sum _{i=1}^m{\alpha }_i{y}_i=0,\\ & {\alpha }_i⩾0,i=1,2,\dots ,m.\end{array} \end{gathered}$$

模型为：$$\begin{array}{rl}f(\mathbf{x})& ={\mathbf{w}}^{\mathrm{T}}\varphi (\mathbf{x})+b\\ & =\sum _{i=1}^m{\alpha }_i{y}_i\varphi ({\mathbf{x}}_i{)}^{\mathrm{T}}\varphi (\mathbf{x})+b\\ & =\sum _{i=1}^m{\alpha }_i{y}_i\kappa (\mathbf{x},{\mathbf{x}}_i)+b.\end{array}$$

$\kappa (\ast ,\ast )$ 就是核函数。常用核函数如下。

实验

导入数据

from numpy import *
from time import sleep
import random
def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat
def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj
dataMat,classLabels=loadDataSet('06testSet.txt')

进行数据可视化。

import matplotlib.pyplot as plt
import numpy as np
x = [row[0] for row in dataMat]
y = [row[1] for row in dataMat]
plt.scatter(x, y, c=['blue' if row == 1 else 'red' for row in classLabels])
plt.xlabel('X')
plt.ylabel('Y')

plt.show()

上图可以看出，数据是合法的：可以有一条直线（这条直线就是所需求的超平面）进行划分

实现SMO函数：

创建⼀个alpha向量并将其初始化为0向量
当迭代次数⼩于最⼤迭代次数时（外循环）：
    对数据集中的每个数据向量（内循环）：
        如果该数据向量可以被优化：
            随机选择另外⼀个数据向量
            同时优化这两个向量
        如果两个向量都不能被优化，退出内循环
    如果所有向量都没被优化，增加迭代数⽬，继续下⼀次循环

  Cell In[10], line 1
    创建⼀个alpha向量并将其初始化为0向量
                        ^
SyntaxError: invalid character in identifier

def smoSimple(dataMatIn, classLabels, C, toler, maxIter=50):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    alphas = mat(zeros((m, 1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])  # if checks if an example violates KKT conditions
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i, m)
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L == H:
                    # print("L==H")
                    continue
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - dataMatrix[i, :] * dataMatrix[i, :].T - dataMatrix[
                     j,:] * dataMatrix[j,:].T
                if eta >= 0:
                    # print("eta>=0")
                    continue
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    # print("j not moving enough")
                    continue
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[i, :].T - labelMat[
                    j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].T
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - labelMat[
                    j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].T
                if (0 < alphas[i]) and (C > alphas[i]):
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2) / 2.0
                alphaPairsChanged += 1
                print("iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas


b, alphas = smoSimple(dataMat, classLabels, 0.6, 0.001, maxIter=50)
b=b[0,0]

部分运行结果：

iter: 0 i:1, pairs changed 1
iter: 0 i:3, pairs changed 2
iter: 0 i:4, pairs changed 3
iter: 0 i:8, pairs changed 4
iter: 0 i:10, pairs changed 5
iter: 0 i:14, pairs changed 6
iter: 0 i:29, pairs changed 7
iter: 0 i:55, pairs changed 8
iteration number: 0
iter: 0 i:1, pairs changed 1
iter: 0 i:24, pairs changed 2
iter: 0 i:27, pairs changed 3
iter: 0 i:29, pairs changed 4
iter: 0 i:55, pairs changed 5
iter: 0 i:76, pairs changed 6
iteration number: 0
iter: 0 i:8, pairs changed 1
iter: 0 i:55, pairs changed 2
iter: 0 i:84, pairs changed 3
iteration number: 0
iter: 0 i:4, pairs changed 1
iter: 0 i:11, pairs changed 2
iter: 0 i:46, pairs changed 3
iter: 0 i:52, pairs changed 4
iter: 0 i:54, pairs changed 5
iteration number: 0
iter: 0 i:1, pairs changed 1
iter: 0 i:22, pairs changed 2
iter: 0 i:24, pairs changed 3
iter: 0 i:30, pairs changed 4
iteration number: 0
iter: 0 i:10, pairs changed 1
iteration number: 0
iter: 0 i:18, pairs changed 1
...
iteration number: 38
iteration number: 39
iteration number: 40
iteration number: 41
iteration number: 42
iteration number: 43
iteration number: 44
iteration number: 45
iteration number: 46
iteration number: 47
iteration number: 48
iteration number: 49
iteration number: 50

输出b和alphas

print(b)
print(alphas)

部分运行结果：

-3.7617698556024015
[[0.        ]
	...
 [0.        ]
 [0.12836263]
 [0.        ]
 [0.        ]
 [0.        ]
	 ...
 [0.        ]
 [0.23620881]
 [0.        ]
	...
 [0.        ]]

计算超平面向量$\mathbf{w}={\sum }_{i=1}^m{\alpha }_i{y}_i{\mathbf{x}}_i$

def GetW(alphas,classLabels,dataMat):
    w=np.zeros(2)
    sv=[]
    for i in range(len(dataMat)):
        if alphas[i,0]>0:
            sv.append(i)
            w+= alphas[i,0]*classLabels[i]*np.array(dataMat[i])
    return sv,w
sv,w= GetW(alphas,classLabels,dataMat)
print(sv)
print(w)

[17, 29, 55]
[ 0.80226906 -0.27808456]

由$w_{0}x+w_{1}y+b=0$ 可得$y=\frac{-w_{0}x-b}{w_1}$

x = [row[0] for row in dataMat]
y = [row[1] for row in dataMat]
color=['blue' if classLabel == 1 else 'red' for classLabel in classLabels]
for i in sv:
    color[i]='black'
plt.scatter(x, y, c=color)
plt.xlabel('X')
plt.ylabel('Y')

x = np.linspace(min(x), max(x), 100)
y =(-w[0]*x-b)/w[1]

plt.plot(x, y)

plt.show()

如上图所示，将支持向量标记为黑点，其余平面两边的点分别标记为红点和蓝点。
通过该支持向量机算法，我们可以求得支持向量、求得超平面方程，将数据进行划分。⽀持向量机的泛化错误率较低，也就是说它具有良好的学习能⼒，且学到的结果具有很好的推⼴性。