支持向量机模型原理推导及python实现

SVM

支持向量机（support vector machines, SVM）,是一种二分类模型，属于判别模型。

基本原理是在特征空间中找到一个间隔最大的分离超平面，将正负类分开。

1. 基本模型定义：

（硬间隔）支持向量机：

假设给定线性可分训练集 $T = \{(x_1,y_1), (x_2. y_2), \cdots ,(x_N, y_N)\}$

通过间隔最大化或者等价地求解相应的凸二次规划问题，从而学习到的分隔超平面为： $w^* · x + b^* = 0$

以及相应的分类决策函数为 $f(x) = sign(w^* · x + b^* )$

以上模型即（硬间隔）支持向量机

其中

$w^*, b^*$ 均为模型参数

$w^*$ 表示分离超平面的法向量

$b^*$ 为分离超平面的截距

支持向量机的目标是最大化分超平面与最近的二分类点的距离

样本空间中任意一点到分离超平面的距离推导：

设一点 $x_0$ , 求它到超平面 $S：w · x + b = 0$ 的距离，
设 $x_0$ 在 $S$ 上的投影点为 $x_1$ , 则有 $w ·x_1 + b = 0$ ，
同时向量 $\overrightarrow{x_0x_1}$ 与超平面的法向量 $w$ 平行，故有，
$|w ·\overrightarrow{x_0x_1} | = |w ·| |\overrightarrow{x_0x_1} | = ||w|| d$ 其中 d 是我们要求的距离，而
$w ·\overrightarrow{x_0x_1} = w · x_0 - w · x_1$
代入 $w · x_1 + b = 0$ 得到 $w ·\overrightarrow{x_0x_1} = w · x_0 + b$
即 $||w|| d = |w · x_0 + b| \Rightarrow d = \frac{1}{||w||} |w· x + b|$

因此支持向量机的目标就是通过特征空间，找到 $（w, b）$ 来最大化这个距离d

2. 函数间隔和几何间隔

函数间隔定义为

对于给定的训练数据集 $T$ 和超平面 $(w·x+b = 0)$ , 定义样本点 $(x_i, y_i)$ 到超平面的函数间隔为：

$\hat \gamma = y_i(w · x_i + b)$

PS：定义训练数据集 $T$ 关于超平面 $(w,b)$ 的函数间隔为所有样本点 $(x_i,y_i)$ 的函数间隔中的最小值即 $\hat \gamma = min (\hat{\gamma_i} ), i =1,2,\cdots, N$

几何间隔定义为：

对于给定的训练数据集 $T$ 和超平面 $(w·x+b = 0)$ , 定义样本点 $(x_i, y_i)$ 到超平面的几何间隔为：

$\gamma = y_i(\frac{w}{||w||} · x_i + \frac{b}{||w||})$

定义训练数据集 $T$ 关于超平面 $(w,b)$ 的几何间隔为所有样本点 $(x_i,y_i)$ 的几何间隔中的最小值即 $\gamma = min ({\gamma_i} ), i =1,2,\cdots, N$

所以当||w|| = 1 时函数间隔与几何间隔等价。

3. 几何间隔最大化

支持向量机的基本目标是找到能够正确划分数据集并且最大化几何间隔的超平面。

和感知机原理不同，感知机中能够正确划分数据集的分离超平面可能会有很多。

SVM中只有一个符合目标的超平面

数学表示为约束最优化问题：
$max_{(w,b)} \gamma$

$s.t. y_i(\frac{w}{||w||} · x_i + \frac{b}{||w||}) >= \gamma, i = 1, 2, \cdots, N$

考虑到几何间隔和函数间隔的关系, 有：

$max_{(w,b)} \frac{\hat\gamma}{||w||}$

$s.t. y_i(w · x+b) > \hat\gamma, i = 1,2,\cdots, N$

假设将 $w$ 和 $b$ 按比例改变为 $\lambda w$ 和 $\lambda b$ ,此时函数间隔将变为 $\lambda \hat\gamma$ ,因此函数间隔 $\hat \gamma$ 的取值对于上述最优化问题没有影响。

不妨取 $\hat\gamma =1$

意识到最大化 $\frac{1}{||w||}$ 和最小化 $\frac{1}{2} ||w||^2$ 是等价的

则上述最优化问题变为

$min_{(w,b)} \frac{1}{2} ||w|| ^2$

$s.t. y_i(w·x+b)-1 \geq 0, i=1,2,\cdots,N$

这是个凸二次规划问题。

可以采用拉格朗日乘子法，将其变换成对偶问题。

4. 拉格朗日法求解凸二次规划问题

定义拉格朗日函数为：

$L(w, b, \alpha) = \frac{1}{2} ||w||^2 - \sum_{i =1}^{N} \alpha_i(y_i(w · x_i + b) -1)$

其中

$\alpha = (\alpha_1 ,\alpha_2, \cdots, \alpha_N)^T$ 为拉格朗日乘子向量

分别对 $w,b$ 并令其 = 0得：

$\nabla _w L(w, b, \alpha) = w - \sum_{i=1}^{N}\alpha_iy_ix_i = 0$

$\nabla _b L(w, b, \alpha) = -\sum_{i =1}^{N}a_iy_i= 0$

得 $w = \sum_{i=1}^{N}\alpha_iy_ix_i$

$\sum_{i =1}^{N}a_iy_i = 0$

代回拉格朗日函数

$L(w,b, \alpha) =\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_iy_ix_i \alpha_j y_j x_j - \sum_{i =1}^{N} \alpha_i(y_i(\sum_{j=1}^{N} \alpha_j y_j x_j) · x_i+ b-1)$

$= -\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_iy_ix_i \alpha_j y_j x_j + \sum_{i =1}^{N} \alpha_i$

即转换后得对偶问题为

$max -\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_iy_ix_i \alpha_j y_j x_j + \sum_{i =1}^{N} \alpha_i$

$s.t. \sum_{i=1}^{N} \alpha_{i} y_{i} = 0$

$\alpha_i \geq 0, i=1,2,\cdots,N$

将目标函数由求极大转为求极小，得到与之等价的对偶问题为

$min \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_iy_ix_i \alpha_j y_j x_j -\sum_{i =1}^{N} \alpha_i$

$s.t. \sum_{i=1}^{N} \alpha_{i} y_{i} = 0$

$\alpha_i \geq 0, i=1,2,\cdots,N$

相应的KKT条件为：

$\alpha_i \geq 0, i=1,2,\cdots,N ;$

$y_i(w · x_i +b -1) \geq 0;$

$\alpha_i(y_i (w · x_i+b) -1) = 0.$

5. 软间隔支持向量机

如果训练样本含有噪声，这个时候支持向量机可能就不存在这个超平面可以将两个样本分开或者达到比较好的最大间隔，为了解决这个问题引入一个新概念：软间隔。

即： $y_i(w· x_i +b) \geq 1 - \xi_ i$

同时我们的约束变为

$min_{(w,b)} \frac{1}{2} ||w|| ^2 + C \sum _{i=1}^{N} \xi_i$

$s.t. y_i(w·x+b) \geq 1-\xi_i , i=1,2,\cdots,N$

$\xi_i \geq 0, i = 1,2,\cdots, N$

同之前硬间隔一样，变换为对偶问题，最后变成：

$min_{\alpha} \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_i \alpha_j y_i y_j(x_i ·x_j）-\sum_{i =1}^{N} \alpha_i$

$s.t. \sum_{i=1}^{N} \alpha_{i} y_{i} = 0$

$0 \leq \alpha_i \leq C , i=1,2,\cdots,N$

所以我们只要求解满足约束的最优解的 $\alpha$ 出来就计算 $w,b$ 就可以找到超平面了。

6. 核函数

实际上到此我们已经掌握了SVM，但是如果我们就这样直接去寻找最优解，速度会非常慢, 这里的核技巧和后面序列最优化算法（SMO）都是为了改善SVM的速度而提出的。

$min_{\alpha} \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_i \alpha_j y_i y_j(x_i ·x_j）-\sum_{i =1}^{N} \alpha_i$
$s.t. \sum_{i=1}^{N} \alpha_{i} y_{i} = 0$
$0 \leq \alpha_i \leq C , i=1,2,\cdots,N$

上式中 $x_i,x_j 均为向量$ ，如果直接计算，对于整个训练集而言，计算起来会有巨大的时间开销，核函数可以减少这里时间开销。

$K(x , z) = \phi(x) · \phi(z) = (x, z)$

即向量 $x,z$ 的点积可以用 $K(x,z)$ 来表示

比如高斯核：

$K(x , z) =exp (- \frac{||x-z||^2}{2 \sigma ^2})$

因此使用的时候我们把 $(x_i, x_j)$ 替换成 $K(x_i,x_j)$ 就好

7. 序列最优化算法

将上一部分的核函数替换掉点积，我们得到:

$min_{\alpha} \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_i \alpha_j y_i y_j K(x_i ·x_j）-\sum_{i =1}^{N} \alpha_i$

$s.t. \sum_{i=1}^{N} \alpha_{i} y_{i} = 0$

$0 \leq \alpha_i \leq C , i=1,2,\cdots,N$

我们需要做的就是计算 $\alpha$

由于 $\alpha = \{\alpha_1, \alpha_2, \cdots, \alpha_N\}$ ,直接求解N个未知量难度很大

序列最优化算法就是为了解决这个问题。

其核心思想就是，先将 $\alpha_1，\alpha_2$ 视为变量, 将其他的 $\alpha_i$ 全部当成参数，然后迭代使得 $\alpha_1$ 满足KKT条件。

由第一个约束有： $\alpha_1 = - \frac{1}{y_1} \sum_{i = 2}^{N} \alpha_i y_i$

所以将原式中的变量 $\alpha_1, \alpha_2$ 显式写出，最终整理得如下：

$min_{(\alpha_1,\alpha_2)} \frac{1}{2} K_{11} \alpha_1^2 +\frac{1}{2} K_{22} \alpha_2^2 + y_1y_2K_{22} \alpha_1\alpha_2 - (\alpha_1+\alpha_2) + y_1\alpha_1 \sum_{i=3}^{N} y_i \alpha_i K_{i1} + y_2\alpha_2 \sum_{i=3}^{N} y_i \alpha_i K_{i2}$

$s.t. \alpha_1 y_1+ \alpha_2 y_2 = \sum_{i=3}^{N} y_i \alpha_i = \zeta$

$0 \leq \alpha_i \leq C, i =1, 2$

根据书中表述，假设上述问题得初始可行解为 $(\alpha_1^{old}, \alpha_2^{old})$ ，最优解为 $(\alpha_1^{new}, \alpha_2^{new})$ ,则有：

$L \leq \alpha_2^{new} \leq H$

其中 $L$ 和 $H$ 是 $\alpha_2^{new}$ 所在对角线端点的界。

如果 $y_1 \neq y_2 \Rightarrow \alpha_1 -\alpha_2 = k$ 则 $L = max(0, \alpha_2^{old} - \alpha_1^{old}), H=min(C,C+\alpha_2^{old} - \alpha_1^{old})$

如果 $y_1 = y_2 \Rightarrow \alpha_1 +\alpha_2 = k$ ,则 $L = max(0, \alpha_2^{old} + \alpha_1^{old}-C), H=min(C,\alpha_2^{old} + \alpha_1^{old})$

假设在沿着约束方向未剪辑时 $\alpha_2$ 的最优解为 $\alpha^{new,unc}$

为了表示 $\alpha^{new}$ ,记 $g(x) = \sum_{i=1}^{N} \alpha_i y_i K(x_i, x) + b， E_i = g(x_i) - y_i, i=1,2$

8.变量的选择方法

第一个变量的选择

从训练样本中选择违反KKT条件最严重的点即检查KKT条件

第二个变量的选择

选择最大化 $|E_1- E_2|$ 的点。

9. python代码实现

实现遵循文中表述

import time
import numpy as np
import math
import random

'''
数据集: MNIST
训练集大小: 60000（实际使用1000）
测试集大小: 10000（实际使用100）
-----
运行结果：

the accuracy is : 0.92
time span: 62.679516553878784 s
'''


def loadData(fileName):
    '''
    加载数据集
    :param fileName: 文件路径
    :return: dataList，labelList分别为特征集X和标记Y. 均为list

    '''
    dataList = []  #
    labelList = []

    f = open(fileName, 'r')

    for line in f.readlines():
        curline = line.strip().split(',')
        '''
        这里考虑到我用的文件是csv格式，所以用split(',')

        Mnsit有0-9十个标记
        文件每行开头第一个数字为该行的label标记
        这里为了简化，限定二分类任务，所以将标记0的作为1(正例)，其余为0(反例)
        '''
        if (int(curline[0]) == 0):
            labelList.append(1)
        else:
            labelList.append(0)

        '''
        加入特征集X
        由于开头第一个数字为标记，故从下标 1 开始
        这里转为int类型
        /255 是为了归一化，有效减少数字爆炸。
        '''
        dataList.append([int(num) / 255 for num in curline[1:]])

    # 读取完毕关闭文件
    f.close()

    # 返回数据集和标记
    return dataList, labelList

class SVM_CLF:
    def __init__(self, X_train,y_train, sigma = 10, C = 200, toler = 0.001):
        """
        相关参数初始化
        :param X_train: 训练数据集
        :param y_train: 训练数据集标记
        :param sigma: 高斯核中的\sigma
        :param C: 惩罚项系数
        :param toler: 松弛变量
        """
        self.trainData = np.mat(X_train)
        self.trainLabel = np.mat(y_train).T
        self.m, self.n = np.shape(self.trainData) #训练及大小，循环中会用到
        self.sigma = sigma                        #高斯核参数
        self.C = C
        self.toler = toler
        self.k = self.calcKernel()                 #初始化时提前计算高斯核
        self.b = 0                                #初始化偏置项为0
        self.alpha = [0] * self.m
        self.E = [0 * self.trainLabel[i,0] for i in range(self.m)]
        self.supportVectorIndex = []


    def calcKernel(self):
        """
        计算核函数，本例使用的是高斯核
        :return: 高斯核矩阵
        """
        #初始化高斯核矩阵 大小=[m * m] ，m为训练集样本数量
        k_matrix = np.zeros((self.m, self.m))
        for i in range(self.m):
            X = self.trainData[i,:]
            for j in range(i, self.m):
                Z = self.trainData[j,:]
                XZ = (X - Z) * (X - Z).T
                result = np.exp(-1 * XZ / (2 * self.sigma **2))
                k_matrix[i][j] = result
                k_matrix[j][i] = result

        return k_matrix

    def isSatisfyKKT(self, i):
        """
        判断alpha i 是否符合KKT条件
        i为下标
        :return: True (满足)或者 False (不满足)
        """
        gxi = self.calc_gxi(i)
        yi = self.trainLabel[i]
        if (math.fabs(self.alpha[i] < self.toler)) and (yi * gxi >= 1):
            return True
        elif(math.fabs(self.alpha[i] - self.C) < self.toler) and (yi * gxi )<=1:
            return True
        elif(math.fabs(self.alpha[i]) > -self.toler) and (self.alpha[i] < self.C)\
            and (math.fabs(yi * gxi) - 1) < self.toler:
            return True
        return False

    def calc_gxi(self, i):
        gxi = 0
        index = [i for i,alpha in enumerate(self.alpha) if alpha != 0]
        for j in index:
            gxi += self.alpha[j] * self.trainLabel[j] * self.k[j][i]

        gxi += self.b

        return gxi

    def calc_EI(self, i):
        '''
        计算Ei
        根据书中章节7.4.1 两个变量二次规划的求解方法 中 式7.105
        :param i: E的下标
        :return: E2,和下标
        '''
        gxi = self.calc_gxi(i)
        return gxi - self.trainLabel[i]

    def get_Alpha2(self, E1, i):
        '''
        选择第二个alpha
        :param E1:
        :param i:
        :return:
        '''
        E2 = 0
        maxE1_E2 = -1
        E2_index = -1

        noZero_E = [i for i ,Ei in enumerate(self.E) if Ei != 0]
        for j in noZero_E:
            E2_tmp = self.calc_EI(j)

            if (math.fabs(E1 - E2_tmp) > maxE1_E2):
                maxE1_E2 = math.fabs(E1 - E2_tmp)
                E2 = E2_tmp
                E2_index = j

            if E2_index == -1:
                E2_index = i
                while E2_index == i:

                    E2_index = random.randint(0, self.m)
                E2 = self.calc_EI(E2_index)

            return E2, E2_index


    def train(self, Maxiter = 100):
        '''
        用训练数据集学习模型
        :param Maxiter: 最大迭代次数
        :return: w， b

        '''
        iterstep = 0
        is_Alpha_changed = 1
        while(iterstep < Maxiter) and (is_Alpha_changed > 0):

            print('iterstep:%d:%d'%(iterstep, Maxiter))

            iterstep += 1
            is_Alpha_changed = 0
            for i in range(self.m):
                if self.isSatisfyKKT(i) == False:
                    E1 = self.calc_EI(i)
                    E2, j = self.get_Alpha2(E1, i )

                    y1 = self.trainData[i]
                    y2 = self.trainLabel[j]

                    alpha_old_1 = self.alpha[i]
                    alpha_old_2 = self.alpha[j]

                    if y1 != y2:
                        L = max(0, alpha_old_2 - alpha_old_1)
                        H = min(self.C, self.C+ alpha_old_2 - alpha_old_1)

                    else:
                        L = max(0, alpha_old_2 + alpha_old_1 - self.C)
                        H = min(self.C, alpha_old_2 + alpha_old_1)
                    #如果L 和 H相等 ，说明变量无法优化， 直接跳过
                    if L == H: continue

                    #根据书中式7.106
                    # 计算alpha的新值
                    K11 = self.k[i][i]
                    K22 = self.k[j][j]
                    K12 = self.k[i][j]

                    alpha_new_2 = alpha_old_2 + y2 * (E1 - E2) / (k11 + k22 - 2 * k12)
                    #剪切alpha2
                    if alpha_new_2 < L : alpha_new_2 = L
                    elif alpha_new_2 > H: alpha_new_2 = H

                    #根据 7.109式 更新alpha1
                    alpha_new_1 = alpha_old_1 + y1 * y2 * (alpha_old_2 - alpha_new_2)

                    #根据书中“7.4.2 变量的选择方法” 第三步式 7.115 和 7.116
                    b1_new = -E1 - y1 * K11 * (alpha_new_1 - alpha_old_1) \
                             - y2* K21 * (alpha_new_2 - alpha_old_2) + self.b
                    b2_new = -E2 - y1 * K12 * (alpha_new_1 - alpha_old_1) \
                             - y2* K22 * (alpha_new_2 - alpha_old_2) + self.b
                    if(alpha_new_1 > 0) and (alpha_new_1 < self.C):
                        b_new = b1_new
                    elif(alpha_new_2 > 0) and (alpha_new_2 < self.C):
                        b_new = b2_new
                    else:
                        b_new = (b1_new + b2_new) / 2

                    #将各参数更新
                    self.alpha[i] = alpha_new_1
                    self.alpha[j] = alpha_new_2
                    self.b = b_new

                    self.E[i] = self.calc_EI(i)
                    self.E[j] = self.calc_EI(j)

                    if math.fabs(alpha_new_2 - alpha_old_2) > 0.0001:
                        is_Alpha_changed += 1

                    #print("iter: %d i: %d, pairs changed %d" %(iterstep, i, is_Alpha_changed))

                    #全部计算结束后遍历一遍alpha ,找出支持向量。
                    for i in range(self.m):

                        if self.alpha[i] > 0:
                            self.supportVectorIndex.append(i)


    def calc_Single_Kernel(self, x1, x2):
        '''
        单独计算核函数
        在predict的时候用到

        :param x1:
        :param x2:
        :return:
        '''
        x1_x2 = (x1 - x2) * (x1 - x2).T
        result = np.exp(- x1_x2 / (2 * self.sigma ** 2))
        return result

    def predict(self, x):
        '''
        预测样本
        :param x: 预测的样本
        :return: 预测结果
        '''

        result = 0
        for i in self.supportVectorIndex:
            tmp = self.calc_Single_Kernel(np.mat(x), self.trainData[i,:])
            result += self.alpha[i] * self.trainLabel[i] * tmp

        result += self.b

        return np.sign(result)

    def test(self, testData, testLabel):
        '''
        测试测试集
        :param testData: 测试数据集
        :param testLabel: 测试标记集
        :return:  正确率

        '''
        errorCount = 0
        for i in range(len(testData)):
            print('testing : %d: %d'%(i, len(testData)))

            result = self.predict(testData[i])
            if result != testLabel[i]:
                errorCount += 1

        return 1 - errorCount/len(testLabel)


if __name__ == '__main__':
    start = time.time()

    print('read trainSet')
    trainData, trainLabel = loadData('D:/PythonLearn/MLA/mnist_train.csv')

    print('read testSet')
    testData, testLabel = loadData('D:/PythonLearn/MLA/mnist_test.csv')

    print('init SVM ')
    svm = SVM_CLF(trainData[:1000], trainLabel[:1000], sigma=10, C = 200,toler = 0.001)

    print('training')
    svm.train()

    print('testing')
    accuracy = svm.test(testData[:100], testLabel[:100])

    end = time.time()

    print('the accuracy is :', accuracy)
    print('time span:', end - start, 's')