SVM源码实现 李航《统计学习方法》

SVM介绍

  SVM是一种二分类器，当然现在也出现了多分类SVM，针对训练数据数据可分情况将分为以下几种SVM分类算法：

1. 线性可分：采用硬间隔最大化训练出线性可分SVM。
2. 近似线性可分：采用软间隔最大化训练出线性支持SVM。
3. 非线性可分：采用核函数核技巧的方法训练出非线性支持SVM。当然我们知道针对非线性可分数据还存在用感知机方法以及多层感知机（神经网络）方法来训练，但本文仅介绍基于核技巧的SVM。

  支持向量机的学习目的是使得间隔最大化，最后可转化为求解一个凸二次规划问题，本文主要采用QP算法来求解凸二次规划，当数据量特别大时，为了提升计算速度，我们才会考虑使用SMO算法实现求解。

SVM算法理论知识

  在感知机里，我们知道求解分离超平面的策略是误分类最小；而SVM求解分离超平面时利用间隔最大化。

函数间隔

  定义超平面（$\omega ,b$）关于样本点（$x_i,y_i$）的函数间隔为
${\hat{\gamma }}_i=y_i(\omega \cdot x_i+b)$
  超平面关于整个数据集的函数间隔为所有样本点的函数间隔的最小值，即
$\hat{\gamma }=\underset{i=1,\cdots ,N}{\min }{\hat{\gamma }}_i$

间隔最大化等价问题

  1、硬间隔最优化问题：
$\underset{\omega ,b}{\min }\frac{1}{2}||\omega |{|}^2$
$s.t.1-y_i(\omega \cdot x_i+b)⩽0,i=1,2,\cdots ,N$
  与之等价的对偶最优化问题：
$\underset{\alpha }{\min }\frac{1}{2}{\sum }_{i=1}^N$${\sum }_{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j(x_i\cdot x_j)$$-{\sum }_{i=1}^N{\alpha }_i$
$s.t.{\sum }_{i=1}^N{\alpha }_i{y}_i=0$
${\alpha }_i⩾0,i=1,2,\cdots ,N$

  2、软间隔最优化问题：线性不可分意味某些样本点$(x_i,y_i)$不满足函数间隔大于等于1的约束条件，故而引入松弛变量$\xi ⩾0$，使得函数间隔加上松弛变量大于等于1。
  此时最优化问题为：
$\underset{\omega ,b,\xi }{\min }\frac{1}{2}||\omega |{|}^2+C{\sum }_{i=1}^N{\xi }_i$
$s.t.1-{\xi }_i-y_i(\omega \cdot x_i+b)⩽0,i=1,2,\cdots ,N$
${\xi }_i⩾0,i=1,2,\cdots ,N$
  其对偶问题为：
$\underset{\alpha }{\min }\frac{1}{2}{\sum }_{i=1}^N$${\sum }_{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j(x_i\cdot x_j)$$-{\sum }_{i=1}^N{\alpha }_i$
$s.t.{\sum }_{i=1}^N{\alpha }_i{y}_i=0$
$0⩽{\alpha }_i⩽C,i=1,2,\cdots ,N$

  3、非线性SVM最优化问题：将软间隔最优化问题中的内积$(x_i\cdot x_j)$用核函数$K(x_i,x_j)$代替，其对偶问题也是如此。

原问题和对偶问题最优解关系

  以软间隔最优化问题即线性支持向量机为例，找出上面原问题最优解和其对偶问题最优解间的关系。首先我们明确算法总的流程如下：

1. 构造原问题模型，求其最优解为${\omega }^{\ast },b^{\ast }$
2. 得到分离超平面${\omega }^{\ast }x+b^{\ast }=0$，所以最终的分类决策函数：$f(x)=sign({\omega }^{\ast }x+b^{\ast })$

  对于上述流程中第1步，求解原问题最优化问题，应用到拉格朗日对称性，通过求解对偶问题得到原始问题的最优解，本文求解对偶问题主要用到QP算法。
  软间隔最优化问题的拉格朗日函数为：
$L(\omega ,b,\xi ,\alpha ,\mu )=\frac{1}{2}||\omega |{|}^2+C{\sum }_{i=1}^N{\xi }_i+{\sum }_{i=1}^N{\alpha }_i(1-{\xi }_i-y_i(\omega \cdot x_i+b))-{\sum }_{i=1}^N{\mu }_i{\xi }_i$
其中${\alpha }_i⩾0,{\mu }_i⩾0$为拉格朗日乘子。
（1）原问题的对偶问题为极大极小问题：
$\underset{\alpha }{\max }\underset{\omega ,b,\xi }{\min }L(\omega ,b,\xi ,\alpha ,\mu )$
  先求$\underset{\omega ,b,\xi }{\min }L(\omega ,b,\xi ,\alpha ,\mu )$，得到当
$\omega ={\sum }_{i=1}^N{\alpha }_i{y}_i{x}_i$
${\sum }_{i=1}^N{\alpha }_i{y}_i=0$
$C-{\alpha }_i-{\mu }_i=0$
  同时成立时，拉格朗日对$\omega ,b,\xi$取的极小值。
  下面对$\underset{\omega ,b,\xi }{\min }L(\omega ,b,\xi ,\alpha ,\mu )$求$\alpha$的极大，通过转换，最终得到原问题的对偶问题等价形式，相当于求$\underset{\alpha }{\min }\frac{1}{2}{\sum }_{i=1}^N$${\sum }_{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j(x_i\cdot x_j)$$-{\sum }_{i=1}^N{\alpha }_i$
（2）假设已经通过QP算法或者SMO算法求出的对偶问题最优解为${\alpha }^{\ast }=({\alpha }_1^{\ast },{\alpha }_2^{\ast },\cdots ,{\alpha }_N^{\ast }{)}^T$，通过回代计算得到：
${\omega }^{\ast }={\sum }_{i=1}^N{\alpha }_i^{\ast }{y}_i{x}_i$
  选择${\alpha }^{\ast }$满足条件$0<{\alpha }_j^{\ast }<C$的一个分量${\alpha }_j^{\ast }$，由于原问题为一个凸二次规划问题，故而解满足KKT条件，根据KKT条件可计算：
$b^{\ast }=y_j-{\sum }_{i=1}^N{y}_i{\alpha }_i^{\ast }(x_i\cdot x_j)$
  注意：$b^{\ast }$是根据满足$0<{\alpha }^{\ast }<C$这个条件的一个分量确定的，其中的$y_i,x_i$与之对应，由于满足这个条件的${\alpha }_j^{\ast }$不止一个，所以$b^{\ast }$不唯一。
（3）分离超平面和决策函数计算出来后，SVM模型便训练出来了，上面只是理论推导，具体代码实现见后面。

核函数

  当给定的数据完全非线性可分时，我们考虑通过提高其维度使其能够线性可分，也就是存在一个非线性映射$\varphi (x)$，将原空间映射到希尔伯特空间。核函数本质就是一个内积，是两个非线性映射的内积，主要包括“升维”和“求内积”两步。
  核函数:$K(x,z)=\varphi (x)\cdot \varphi (z)$
  核技巧：软间隔只能求解少数线性不可分数据，当数据完全线性不可分，我们要找到一个非线性映射$\varphi (x)$来提升维度使数据线性可分是十分困难的，然而我们在求解该最优化问题时并非需要找到这样的映射$\varphi (x)$，我们目标函数里只有原样本间内积，所以我们用核函数来求映射后样本内积来替换原样本内积就行了（使用核函数就默认存在这样的映射使得数据维度变高，使得样本线性可分）。
  常用核函数：一般自己构造核函数十分困难，所以文献中存在这许多常用的核函数，我们只需要根据自己需求选择使用就行了。常用核函数有：

• 线性核：$K(x,z)=x\cdot zx,z为向量，表示求其内积$
• 高斯核：$K(x,z)=exp(-\frac{||x-z|{|}^2}{2{\sigma }^2})$   径向基核(RBF)的特殊形式.
• 多项式核：$K(x,z)=(a\ast x\cdot z+c{)}^{p}a,c为常数$
• Sigmoid核：$K(x,z)=tanh(\alpha \ast x\cdot z+c)$   常用于深度学习中
• 更多核函数参考核函数详解

  正定核检验：由于我们通过找到映射来构造核函数十分困难，所以我们只用判断我们遇到的一个函数能否作为一个核函数来使用，核函数的检验方法就是判断核函数关于任意样本的Gram矩阵是否为半正定矩阵，如果是，那么这个核函数是正定核。

下面代码是写一个后面要用到的核函数类：

import numpy as np
import numpy.linalg as la  #包含线性代数的函数，求逆矩阵、求特征值、行列式等等。
class kernel(object):
   @staticmethod
   def linear():
       return lambda x, z: np.inner(x, z)

   @staticmethod
   def gaussian(sigma):
       return lambda x, z: np.exp(-la.norm(x-z)**2/(2*sigma**2))  #la.norm表示求2范数

   @staticmethod
   def _polynomial(p, c):
       return lambda x, z: (np.inner(x, z)+c)**p

   @classmethod
   def inhomo_polynomial(cls, p):
       return cls._polynomial(p, 1.0)

   @classmethod
   def homo_polynomial(cls, p):
       return cls._polynomial(p, 0.0)

   @staticmethod
   def sigmoid(gamma, c):
       return lambda x, z: np.tanh(gamma*np.inner(x, z)+c)

   @staticmethod
   def RBF(gamma=10):
       return lambda x, y: np.exp(-gamma*la.norm(x-y))

SVM算法步骤及其Python源码

  我们现在的目的是实现对非线性支持向量机的对偶问题进行求解，首先我们要了解利用python的cvxopt.solvers来求解凸二次规划问题，我们一般要将我们的QP问题转化成标准形式，其重点是将其写成矩阵形式：
  比如针对本文的非线性SVM中的对偶问题，我们得到他们的矩阵形式为：
$P=y{y}^{T}K$  其中K为核函数关于X的Gram矩阵，y为样本标签为列向量
$q=(-1,-1,\cdots ,-1{)}_{1\times N}^T$
$A=y^T$
$b=0$
$G1=-I$  其中$I_{N\times N}$为单位矩阵
$h1=[0{]}_{N\times 1}$  $h1$为0列向量
$G2=I$
$h2=(C,C,\cdots ,C{)}_{1\times N}^T$
  Step1:第一步是先写出根据核函数写出关于样本的Gram矩阵，所以定义以下函数(其中默认选择高斯核作为核函数)：

 def _gram_matrix(self, X):
        n_samples, n_features = X.shape
        K = np.zeros((n_samples, n_samples))
        for i, x_i in enumerate(X):
            for j, x_j in enumerate(X):
                K[i, j] = self._kernel(x_i, x_j)
        return K

  Step2:第二步是利用python的cvxopt.solvers来求解凸二次规划问题（对偶问题）：

    def _compute_multipliers(self, X, y):
        n_samples, n_features = X.shape
        K = self._gram_matrix(X)
        p = cvxopt.matrix(np.outer(y, y)*K)
        q = cvxopt.matrix(-1*np.ones(n_samples))
        A = cvxopt.matrix(y, (1, n_samples))
        b = cvxopt.matrix(0.0)
        G1 = cvxopt.matrix(-1*np.diag(np.ones(n_samples)))
        h1 = cvxopt.matrix(np.zeros(n_samples))
        G2 = cvxopt.matrix(np.diag(np.ones(n_samples)))
        h2 = cvxopt.matrix(self._C*np.ones(n_samples))
        G = cvxopt.matrix(np.vstack((G1, G2)))
        h = cvxopt.matrix(np.vstack((h1, h2)))
        solution = cvxopt.solvers.qp(p, q, G, h, A, b)
        return np.ravel(solution['x'])    

这里返回的结果为对偶问题的一个解，为一个长度为N的行向量。
  Step3:第三步为根据对偶问题的解求出原问题的解，并写出决策函数。

#求决策函数中最优的两个参数,并返回训练好的SVM模型
    def _construct_predictor(self, X, y):
        n_samples, n_features = X.shape
        K = self._gram_matrix(X)
        lagrange_multipliers = self._compute_multipliers(X, y)
        omega_star = np.dot(np.multiply(lagrange_multipliers, np.transpose(y)), X)   #1*N
        alter_lagrange_multipliers = [lagrange_multipliers[i] for i in range(n_samples) if self._C > lagrange_multipliers[i] > 0]
        alter_index = [i for i in range(n_samples) if self._C > lagrange_multipliers[i] > 0]
        alpha_star = alter_lagrange_multipliers[0]
        y_j = y[alter_index[0]]
        x_j = X[alter_index[0], :]
        b_star = y_j - np.dot(np.multiply(lagrange_multipliers, np.transpose(y)), K[:, alter_index[0]])
        return SVMPredictor(kernel = self._kernel, weights = omega_star, bias = b_star)

SVM源码汇总

  根据上面主要部分代码，我们最终将SVM模型（非线性支持向量机）封装成两个类，一个训练类，另一个为一个预测类，完整代码如下：

import numpy as np
import cvxopt.solvers
import kernel

class SVMTrainer(object):
    def __init__(self, kernel, C):
        self._kernel = kernel   #之前要创建一个kernel实例并调用它的方法
        self._C = C

    def _gram_matrix(self, X):
        n_samples, n_features = X.shape
        K = np.zeros((n_samples, n_samples))
        for i, x_i in enumerate(X):
            for j, x_j in enumerate(X):
                K[i, j] = self._kernel(x_i, x_j)
        return K

    def _compute_multipliers(self, X, y):
        n_samples, n_features = X.shape
        K = self._gram_matrix(X)
        p = cvxopt.matrix(np.outer(y, y)*K)
        q = cvxopt.matrix(-1*np.ones(n_samples))
        A = cvxopt.matrix(y, (1, n_samples))
        b = cvxopt.matrix(0.0)
        G1 = cvxopt.matrix(-1*np.diag(np.ones(n_samples)))
        h1 = cvxopt.matrix(np.zeros(n_samples))
        G2 = cvxopt.matrix(np.diag(np.ones(n_samples)))
        h2 = cvxopt.matrix(self._C*np.ones(n_samples))
        G = cvxopt.matrix(np.vstack((G1, G2)))
        h = cvxopt.matrix(np.vstack((h1, h2)))
        solution = cvxopt.solvers.qp(p, q, G, h, A, b)
        return np.ravel(solution['x'])

    #求决策函数中最优的两个参数,并返回训练好的SVM模型
    def _construct_predictor(self, X, y):
        n_samples, n_features = X.shape
        K = self._gram_matrix(X)
        lagrange_multipliers = self._compute_multipliers(X, y)
        omega_star = np.dot(np.multiply(lagrange_multipliers, np.transpose(y)), X)   #1*N
        alter_lagrange_multipliers = [lagrange_multipliers[i] for i in range(n_samples) if self._C > lagrange_multipliers[i] > 0]
        alter_index = [i for i in range(n_samples) if self._C > lagrange_multipliers[i] > 0]
        alpha_star = alter_lagrange_multipliers[0]
        y_j = y[alter_index[0]]
        x_j = X[alter_index[0], :]
        b_star = y_j - np.dot(np.multiply(lagrange_multipliers, np.transpose(y)), K[:, alter_index[0]])
        return SVMPredictor(kernel = self._kernel, weights = omega_star, bias = b_star)

    def train(self, X, y):
        return self._construct_predictor(X, y)

class SVMPredictor(object):
    def __init__(self, kernel, weights, bias):
        self._kernel = kernel
        self._weights = weights
        self._bias = bias

    def predict(self, x):
        result = np.sign(np.dot(self._weights, np.transpose(x))+self._bias)
        return result

该模型主要有两个超参数，一个是kernel（而我们选择的kernel里面还存在其他超参数），另一个为惩罚参数C，我们训练模型的目的即找到最优的参数，其参数要用损失函数关于验证集的大小确定，以高斯核为例，他们参数关系可参考高斯核SVM参数确定.