来理解Support Vector Machine吧！

---

本文内容参考机器学习与自然语言算法公众号，知乎，周志华Machine Learning 等整理而来，原版资料已在文末以Ref标明，如有疏忽敬请谅解。

---

第一部分：感知机的介绍

感知机是个相当简单的模型，但它既可以发展成支持向量机（通过简单地修改一下损失函数）、又可以发展成神经网络（通过简单地堆叠），所以它也拥有一定的地位。
在此，我们先讨论二分类问题，样本类别为正样本 和负样本。

1. 感知机能干嘛

• 感知机能且一定能将线性可分的数据集分开，但是不能处理线性不可分的数据集

什么叫线性可分？在二维平面上、线性可分意味着能用一条线将正负样本分开，在三维空间中、线性可分意味着能用一个平面将正负样本分开。
直观的显示可分与不可分就是这样：

• // gif图片不可显示 请于附加的word中查看

2. 搭建一个感知机模型

为了搭建感知机，我们先要知道数学上的的线性可分是指什么。

• 定义超平面
  $\qquad\Pi:\qquad$$w * x + b = 0$
  其中，w和x都是n维向量，n=2时，有$w_1x_1 + w_2x_2 + b = 0$，这就是大学以前学的直线方程。

• 感知机的线性可分
  对于数据集D，样本标签$y_i$为1则为正样本$:$$w*x+b$$>$0，反之为负样本。感知机只有w和b两个参数，其要做的就是根据样本的信息来逐步更新w和b,从而使得超平面$\Pi$慢慢将数据集D的样本类别分开。

3. 训练感知机，更新参数

为了保证收敛性，我们需要把w初始化为n维0向量，将b初始化为标量0。

import numpy as np
class Perceptron:
    def __init__(self):
        self.w = self.b = None
    def fit(self, x, y, lr=0.01, epoch=1000):
    # 将输入的样本x和标签y转化为NumpyN维数组对象
    x, y = np.asarray(x, np.float32), np.asarray(y, np.float32)
    self.w = np.zeros(x.shape[1])
    self.b = 0.

fit函数中的lr和epoch分别代表梯度下降法中的学习速率和迭代上限。
梯度下降法的步骤：
- 对损失函数求导得到梯度值。
- 梯度是函数值增长最快的方向，我们想让损失函数最小化，就要让参数沿着梯度的反方向更新。
那么感知机的损失函数是什么呢？？
我们知道：
$若w*x+b>0,则L(x,1)=0, L(x,-1)=w*x+b$
$若w*x+b<0,则L(x,1)=-(w*x+b), L(x,-1)=0$
综上：
$\color{#F00}{损失函数可写为L(x,y) = max(-y(w*x+b),0)}$
$\color{#F00}{样本(x,y)被正确分类\Longleftrightarrow y(w*x+b)>0}$
因为只有错分类的点才贡献梯度，所以在训练时，我们挑选一个使$L(x,y)$最大的样本$(x_i,y_i)$,如果$(x_i,y_i)$被正确分类的话，说明所有样本均已正确分类，则停止模型的训练。
$L(x,y) = -y_i(w*x_i+b) = -y_ix_i*w - y_ib$

$$\left. \frac{{\rm d}L}{{\rm d}w} \right| _{(x_i,y_i)} =-y_ix_i, \left. \frac{{\rm d}L}{{\rm d}b} \right| _{(x_i,y_i)} = -y_i$$
代码体现损失函数减小和参数更新，实现线性可分

import numpy as np
for _ in range(epoch):
    # 先计算w*x+b
    y_pred = x.dot(self.w) + self.b
    # 选出使损失函数最大的样本
    idx = np.argmax(np.maximun(0,-y*y_pred))    # SVM这不同,还有间隔的最大化
    # 若该样本被正确分类，则停止模型训练
    if y[idx] * y_pred[idx] > 0:             # SVM这不同
        break
    # 否则，让参数沿着负梯度方向走一步
    delta = lr * y[idx]
    self.w += delta * x[idx]
    self.b += delta

---

第二部分：介绍LinearSVM

在这里，我们先讲SVM的原始版本，其本身只是一个线性模型，只有在应用了核方法后，SVM才会升级成为非线性模型。第三部分我们再讲解带核技巧的SVM。

1. SVM与感知机

第一部分中我们提过，LinearSVM可以说是改了感知机的损失函数而已。感知机中样本被正确分类的标志是$y(w*x+b)>0$，而SVM中样本被正确分类的标志是$y(w*x+b)\geq1$。为什么是大于1，这个1是怎么来的？我们先来看一下感知机决策面的劣势。
如(word图1)中，感知机要求样本被正确分类，而不要求样本$\color{#F00}{“被很好的分类”}$。那么如何从这些很多的决策面中选择一个最好的，使该划分的超平面对训练样本的局部扰动“容忍”性最好。$\color{#F00}{直观的解决方法就是让超平面离正负样本点的间隔都尽可能大}$。

2. SVM模型搭建

• 首先也是定义超平面
  $\qquad\Pi:\qquad$$w * x + b = 0$
  此处先给出直线外一点到已知直线的距离：

$P(x_0,y_0)到Ax+By+C = 0的距离为：$

$$D = \frac{Ax_0+By_0+C}{\sqrt{A^2+B^2}}$$

其中A,B决定了直线的法线方向。
所以有，样本空间中任意点$x_i$到超平面$\Pi$的距离可写为

$$r = \frac{|w * x_i + b|}{||w||}$$
硬间隔最大化
假设存在超平面$\Pi$能将训练样本正确分类，即有：对于$(x_i,y_i)$$\in$D，若$y_i = +1$，则$w*x_i+b>0$；若$y_i = -1$，则$w*x_i+b<0$；

定义一个很好的分类：令

$$\begin{cases} w*x_i+b\geq1, & \text {$y_i = +1$} \\ w*x_i+b\leq1, & \text{$y_i = -1$} \end{cases}$$
如下图，$\color{#F00}{距离超平面最近的这几个训练样本使上式等号成立，它们称为“支持向量(Support Vector)”}$。

所以两个异类支持向量到超平面的距离之和为 $$\gamma = \frac{2}{||w||}$$ ，称$\gamma$为“间隔”。
欲使划分的超平面对训练样本的局部扰动“容忍”性最好，需要使间隔最大。故得到如下的约束优化模型。（见word图2）
软间隔最大化
但是，对于上面的问题，$\color{#F00}{如果数据集线性不可分时，模型必定找不到合适的超平面而一直震荡}$，也就是说$w*x_i+b\geq1$的这个约束太硬了（硬间隔）。缓解该问题的一个方法是允许支持向量在一些样本上出错(见下图和word图3)，为此对于线性不可分的情况引入“软间隔”的概念(word图4)。

返回来，我们再看一下感知机的损失函数是写作： $$\color{#F00}{L(x_i,y_i) =\sum_{i=1}^N max(-y_i(w*x_i+b),0)}$$
而SVM的损失函数呢？
刚才看到了：模型的优化目标为 $$\begin{cases}min\left[ \frac{||w||^2}{2}+C\sum_{i=1}^N\xi_i\right]\\s.t.\quad y_i(w*x_i+b)\geq1-\xi_i; \quad\xi_i\geq0, i=1,2,\cdots,m\end{cases}$$
现在对上面的模型进行转化：
$$y_i(w*x_i+b)\geq1-\xi_i\Longleftrightarrow\xi_i\geq1-y_i(w*x_i+b)$$
为使$\left[\frac{||w||^2}{2}+C\sum_{i=1}^N\xi_i\right]$最小，必有：
$$\begin{cases}1-y_i(w*x_i+b)\geq0时，\xi_i=1-y_i(w*x_i+b)\\ 1-y_i(w*x_i+b)<0时，\xi_i=0\end{cases}$$
也就是$\xi_i = max(1-y_i(w*x_i+b),0)。$
所以,SVM的损失函数就是： $$\color{#F00}{L(x_i,y_i)= \frac{||w||^2}{2}+C\sum_{i=1}^N\max(1-y_i(w*x_i+b),0)}$$
于是综上可以看出，LinearSVM在形式上和感知机的损失函数确实长得很像。

另外：$max(1-y_i(w*x_i+b),0)$就是hinge损失，也可以写作：$[1-y_i(w*x_i+b)]_+$

3. LinearSVM的训练

在这里，由于对偶问题的复杂性，我们“偷懒”一下，使用朴素的随机梯度下降（SGD）来求解参数$w,b$。
首先计算损失函数对参数的导数，得到了下降的梯度之后，再进行参数更新。

代码实现线性不可分情况

import numpy as np
class LinearSVM:
    def __init__(self):
        self.w = self.b = None
    def fit(self, x, y, lr=0.01, epoch=10000):
        # 将输入的样本x和标签y转化为Numpy数组
        x, y = np.asarray(x, np.float32), np.asarray(y, np.float32)
        # 随机选取batch_size个样本
        batch = np.random.choice(len(x), batch_size)
        x_batch, y_batch = x[batch], y[batch]
        self.w = np.zeros(x_batch.shape[1])
        self.b = 0.
        for _ in range(epoch):
            self.w *= 1 + lr
            err = 1-y_batch * self.predict(x_batch, true)
            if np.max(err) <= 0  
            # 全部样本均已分类正确,就截断当前的梯度下降，但是不能终止训练，还需要使超平面间隔最大
                #self.w += (-lr * w)
                continue
            # 因为分类正确的样本,惩罚部分梯度为0,所以我们这里只分类误分类的样本
            mask = err > 0  
            delta = lr * c * y_batch[mask]
            # 取各梯度平均并做进一步梯度下降
            self.w += np.mean(delta * x_batch[mask], axis=0)
            self.b += np.mean(delta)
    def predict(self, x, raw=False)  # 补充一句，当要用到类中定义的属性时，类方法中写上self,相当于this。
        x = np.asarray(x, np.float32)
        y_pred = x.dot(self.w) + self.b
        if raw:
            return y_pred
        return np.sign(y_pred).astype(np.float32)

当然，朴素的梯度下降法是有问题存在的。当数据的scale很大时，模型对参数极为敏感，从而会导致超平面持续的震荡，原因估计是因为数据太大梯度很大导致跑得太快了，解决的方法可以进行数据的归一化$X \leftarrow \frac{X-X_均}{\sqrt{Var(X)} \quad }$。
还有一个做法是：用更好的梯度下降法来代替朴素的SGD，如Adam。

最后，请允许我不加证明地给出两个结论（因为结论直观且证明太长……）：
- 若数据集线性可分，则最大硬间隔分离超平面存在且唯一
- 若数据集线性不可分，则最大软间隔分离超平面的解存在但不唯一，其中：
法向量（$w$）唯一
偏置量（$b$）可能不唯一

---

这里，我们给出LinearSVM的官方代码，由于复杂，只挑选了部分

```
from .base import _fit_liblinear, BaseSVC, BaseLibSVM
from ..base import BaseEstimator, RegressorMixin
from ..linear_model.base import LinearClassifierMixin, SparseCoefMixin, LinearModel
from ..utils import check_X_y
from ..utils.validation import _num_samples
from ..utils.multiclass import check_classification_targets
class LinearSVC(BaseEstimator, LinearClassifierMixin,
SparseCoefMixin):
    def __init__(self, penalty='l2', loss='squared_hinge', dual=True, tol=1e-4, C=1.0, multi_class='ovr', fit_intercept=True, intercept_scaling=1, class_weight=None, verbose=0,
random_state=None, max_iter=1000):
    def fit(self, X, y, sample_weight=None):
        # FIXME Remove l1/l2 support in 1.0 
        msg = ("loss='%s' has been deprecated in favor of "
               "loss='%s' as of 0.16. Backward compatibility"
               " for the loss='%s' will be removed in %s")
        if self.loss in ('l1', 'l2'):
            old_loss = self.loss
            self.loss = {'l1': 'hinge', 'l2': 'squared_hinge'}.get(self.loss)
            warnings.warn(msg % (old_loss, self.loss, old_loss, '1.0'),
                          DeprecationWarning
        if self.C < 0:
            raise ValueError("Penalty term must be positive; got (C=%r)"
                             % self.C)
        X, y = check_X_y(X, y, accept_sparse='csr',
                         dtype=np.float64, order="C")
        check_classification_targets(y)
        self.classes_ = np.unique(y)
        self.coef_, self.intercept_, self.n_iter_ = _fit_liblinear(X, y, self.C, self.fit_intercept, self.intercept_scaling, self.class_weight, self.penalty, self.dual, self.verbose, self.max_iter, self.tol, self.random_state, self.multi_class, self.loss, sample_weight=sample_weight)
        if self.multi_class == "crammer_singer" and len(self.classes_) == 2:
            self.coef_ = (self.coef_[1] - self.coef_[0]).reshape(1, -1)
            if self.fit_intercept:
                intercept = self.intercept_[1] - self.intercept_[0]
                self.intercept_ = np.array([intercept])
        return self
        # ----------------------------------------------------------------------------------------------
from .base import _fit_liblinear
    def _fit_liblinear(X, y, C, fit_intercept, intercept_scaling, class_weight, penalty, dual, verbose, max_iter, tol, random_state=None, multi_class='ovr',loss='logistic_regression', epsilon=0.1, sample_weight=None):
         solver_type = _get_liblinear_solver_type(multi_class, penalty, loss, dual)
    raw_coef_, n_iter_ = liblinear.train_wrap(X, y_ind, sp.isspmatrix(X), solver_type, tol, bias, C, class_weight_, max_iter, rnd.randint(np.iinfo('i').max), epsilon, sample_weight)
    return coef_, intercept_, n_iter_
     # ----------------------------------------------------------------------------------------------
     cimport liblinear
     def train_wrap(X, np.ndarray[np.float64_t, ndim=1, mode='c'] Y, bint is_sparse, int solver_type, double eps, double bias, double C, np.ndarray[np.float64_t, ndim=1] class_weight, int max_iter, unsigned random_seed, double epsilon, np.ndarray[np.float64_t, ndim=1, mode='c'] sample_weight):
        model = train(problem, param)
    # ----------------------------------------------------------------------------------------------
model* train(const problem *prob, const parameter *param)
return model_;

源码如此复杂，还是放弃吧。。。
那现在，再来看一下怎么用：

Examples_1  回归
--------
>>> from sklearn.svm import LinearSVR
>>> from sklearn.datasets import make_regression
>>> X, y = make_regression(n_features=4, random_state=0)  #n_features表示一个样本有4个数据 4*1
>>> regr = LinearSVR(random_state=0)
>>> regr.fit(X, y)
LinearSVR(C=1.0, dual=True, epsilon=0.0, fit_intercept=True, intercept_scaling=1.0, loss='epsilon_insensitive', max_iter=1000, random_state=0, tol=0.0001, verbose=0)
>>> print(regr.coef_) #1*4
[ 16.35750999  26.91499923  42.30652207  60.47843124]
>>> print(regr.intercept_)
[-4.29756543]
>>> print(regr.predict([[0, 0, 0, 0]]))
[-4.29756543]
Examples_2  分两类
--------
>>> from sklearn.svm import LinearSVC
>>> from sklearn.datasets import make_classification
>>> X, y = make_classification(n_features=4, random_state=0)
>>> clf = LinearSVC(random_state=0)
>>> clf.fit(X, y)
LinearSVC(C=1.0, class_weight=None, dual=True, fit_intercept=True, intercept_scaling=1， loss='squared_hinge', max_iter=1000, multi_class='ovr', penalty='l2', random_state=0, tol=0.0001, verbose=0)
>>> print(clf.coef_)
[[ 0.08551385  0.39414796  0.49847831  0.37513797]]
>>> print(clf.intercept_)
[ 0.28418066]
>>> print(clf.predict([[0, 0, 0, 0]]))
[1]   
Examples_3  分多类
--------
# 对于一对一的分类，如果要分n类，则需要构造n*(n-1)/2个分类器
>>> X = [[0], [1], [2], [3]]
>>> Y = [0, 1, 2, 3]
>>> clf = svm.SVC(decision_function_shape='ovo')
>>> clf.fit(X, Y) 
SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0, decision_function_shape='ovo', degree=3, gamma='auto', kernel='rbf', max_iter=-1, probability=False, random_state=None, shrinking=True, tol=0.001, verbose=False)
>>> dec = clf.decision_function([[1]]) # 得到的是矩阵
>>> dec.shape[1] # 4 classes: 4*3/2 = 6 # 矩阵才可以.shape[],得到的是分类器的个数
6
>>> clf.decision_function_shape = "ovr"
>>> dec = clf.decision_function([[1]])
>>> dec.shape[1] # 4 classes
4
>>> print(clf.predict([[5]]))
[3]

行文至此，可以做个小结，不准确的说，SVM它本质上即是一个分类方法，用$w^T+b$定义分类函数，于是求w、b，为寻最大间隔，引出$1/2||w||^2$，继而引入拉格朗日乘子，化为对拉格朗日乘子a的求解（求解过程中会涉及到一系列最优化或凸二次规划等问题），如此，求w，b与求a等价，而a的求解可以用一种快速学习算法SMO。
至于核函数，是为处理非线性情况，若直接映射到高维计算唯恐维度爆炸，故在低维计算，等效高维表现。

---

第三部分：核方法

前面我们说的都是建立在线性分类问题上（线性可分与线性不可分），对于非线性的分类问题，SVM引入了核方法来解决问题。

1. 什么是核方法

往简单里说，核方法是将一个低维的线性不可分的数据映射到一个“高维”的空间、并期望映射后的数据在“高维空间”里是线性可分的。如下图很好的表示了这样的一个过程。

Cover定理从理论上说明了：当空间维数d越大时，其中的点线性可分的概率就越大。所以，在解决线性不可分的问题时，我们要做的就是寻找合适的映射$\phi$，使得数据集在被它映射到高维空间后变得线性可分。
令$\phi(x)$表示将$x$映射后的特征向量，于是在特征向量中划分超平面所对应的模型可表示为$f(x) = w*\phi(x)+b$。类似于上面的硬间隔最大化，优化模型为：

$$\begin{cases}min\frac{||w||^2}{2}\\s.t.\quad y_i(w*\phi(x_i)+b)\geq1; \quad i\leq0, i=1,2,\cdots,m\end{cases}$$
其对偶问题是
$$\begin{cases}max_\alpha\sum_{i=1}^m\alpha_i - \frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_j\color{#F00}{\phi(x_i)^T\phi(x_j)} \\s.t.\quad \sum_{i=1}^m\alpha_iy_i = 0 ; \quad \alpha_i\geq0, i=1,2,\cdots,m \end{cases}$$
可以看到，映射后的特征空间中涉及到计算$\phi(x_i)^T\phi(x_j)$。
但是，现实任务中的数据要复杂得多，直接构造一个恰当的$\phi$难度可能要高于原问题本身，而且有时候特征空间的维数是无穷维，直接计算$\phi(x_i)^T\phi(x_j)$通常是很困难的。为了避开这个难题，设想了这样一个函数：称为核函数。
$$k(x_i,x_j) = <\phi(x_i),\phi(x_j)> = \phi(x_i)^T\phi(x_j)$$
即$x_i与x_j$在特征空间的内积等于它们在原始空间通过核函数$k(x_i,x_j)$计算的结果。

2. 核函数的思想

• 将算法表述成样本点内积的组合：$f(x) = w*\phi(x)+b=\sum_{i=1}^m\alpha_iy_i\phi(x_i)^T\phi(x)+b = \sum_{i=1}^m\alpha_iy_ik(x,x_i)+b$，进而表示为核函数的线性组合。

• 设法找到核函数，就不必直接去计算高维甚至无穷维特征空间中的内积。

• 核函数的引入完成了从线性算法到非线性算法的转换

但是，每一个核函数都能找到一个映射与之对应吗？显而易见这是不可能的，并不是所有的$k(x_i,y_i)$都能拆成$\phi(x_i)^T\phi(x_j)$。1909年Mercer定理给出了解决方案：想让这个核函数总能找到内积映射，也就是说这个核函数总能对应一个特征空间的话，我们就要找到一个对称函数且对应的核矩阵半正定，那么这个函数就可以作为核函数使用，定义的空间称为：再生核希尔伯特空间(RKHS)。
核函数选择的是否合适决定了能否能将样本映射到合适的希尔伯特空间，因此核函数的选择是做非线性分类的重要问题。
如下图，是应用核方法的分类器的分类结果：

3. 常用核函数

• 线性核 $\quad k(x_i,y_i) = x_i^Tx_j$
• 多项式核 $\quad k(x_i,y_i) = (x_i^Tx_j)^d$
• 高斯核（RBF核） $\quad k(x_i,y_i) = exp(-\frac{||x_i-x_j||^2}{2\sigma^2})$
• 拉普拉斯核 $\quad k(x_i,y_i) = exp(-\frac{||x_i-x_j||^2}{\sigma})$
• sigmord核 $\quad k(x_i,y_i) = tanh(\beta x_i^Tx_j+\theta)$

再看一张图：

这个例子从侧面简单说明了SVM使用非线性分类器的优势，而logstic回归（感知机线性分类器）以及决策树模式都是使用了直线方法。
看了这么多，说白了$\color{#F00}{核函数=映射=内积=相似度}$。
其实，核函数的应用场景非常广泛。表示定理中说了：

第四部分：完整代码

from sklearn.datasets import load_digits
from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC
from sklearn.metrics import classification_report
digits = load_digits()
print digits.data.shape #load datas
x_train, x_test, y_train, y_test = train_test_split(digits.data, digits.target, test_size = 0.25, random_state = 33)
print x_train.shape,x_test.shape,y_test.shape,y_train.shape #split train datas and test datas
ss = StandardScaler()  # 标准化数据，使每个维度上的特征数据方差为1，均值为0，使得预测结果不会被某些维度过大的特征值而主导。
x_train = ss.fit_transform(x_train)
x_test = ss.transform(x_test)  #standard datas
lsvc = LinearSVC()  # 线性分类，没有核函数；非线性分类使用SVC(kernel='rbf')或者SVC(kernel='poly')等。
lsvc.fit(x_train,y_train)
y_predict = lsvc.predict(x_test)  #use svc(support vectors classification) to train and predict
print 'The Accurary of linear SVC is',lsvc.score(x_test,y_test)
print classification_report(y_test, y_predict, target_names=digits.target_names.astype(str))
# from sklearn import svm
#
# x = [[0,0],[1,1],[1,0],[13,113]]
# y = [0,1,1,1]
# clf = svm.SVC()
# clf.fit(x,y)
# print clf.predict([[2,2]])   #result
# print clf.support_vectors_   #get support vectors
# print clf.support_           #get indices of support vectors
# print clf.n_support_         #get number of support vectors for each class

---

Ref:
周志华《机器学习》
范筱《python机器学习及其实战》
微信公众号：射命丸咲 机器学习算法与自然语言处理
https://github.com/scikit-learn/scikit-learn/tree/ef5cb84a805efbe4bb06516670a9b8c690992bd7/sklearn/svm 机器学习官方代码
还有很多从博客或者是网上看到后做的笔记，已经找不到是在哪里了，如果侵犯，敬请谅解。