SVM浅显理解

0. 前言

  最近在学习一些机器学习的知识，看的书主要是李航的《统计学习方法》。看到了支持向量机这一章，花了不少时间，把一些基础知识掌握了。写一篇博客记录一下我的理解。
  首先支持向量机有三种分类，分别是线性可分支持向量机、线性支持向量机以及非线性支持向量机。后面的类型比前面的复杂，用到的技术也是逐渐增加。线性可分SVM用到了硬间隔最大化，线性SVM则用到条件松一些的软间隔最大化，非线性SVM则用到核技巧以及软间隔最大化。

1. 线性可分支持向量机

  我们先看最简单的线性可分SVM。首先我们需要知道了解两个定义：函数间隔以及几何间隔。
  函数间隔的定义：

其中y代表分类是否正确，wx+b代表分类预测的置信度。函数间隔越小，说明样本跟超平面越接近。
  几何间隔的定义：

  几何间隔就是我们平时理解的点到直线的距离（二维情况下）。
  需要注意的两个点（后面都会用到）是：

1. 函数间隔可以随着w跟b的倍数增加而增加，取了两倍、三倍的w跟b，超平面是不变的，但是函数间隔会变成原来的两倍、三倍。
2. 函数间隔跟几何间隔的关系是差了一个$\frac{1}{\Vert w\Vert }$。

  接下来进行建模，优化问题写成如下的样子：

  这个公式的意思可以这样理解：所有样本点中，会有一个离超平面最近的点，这个最近的距离就是$\gamma$，其他所有距离都会大于这个数，这就是约束条件。同时我们想要把这个最近的距离尽可能地放大（这就是间隔最大化）。
  然后我们用函数间隔跟几何间隔的关系（上面第二点），把$\gamma$换成$\hat{\gamma }$

  然后我们取$\hat{\gamma }$为1，方便后面计算。为什么可以这样做呢，就是上面的第一点，函数间隔是可以变化的，你可以通过调整w，b的倍数把函数间隔调整到某一个数，因此$\hat{\gamma }$可以是任意一个数。
  其次把问题转化一下，最大化$\frac{1}{\Vert w\Vert }$转化成最小化$\frac{1}{2}{\Vert w\Vert }^2$，这个证明我没去看过，但是也很容易理解，我们可以把两个函数画出来，就可以看出两个问题是一样的。于是优化问题变成这个样子。
  这是一个凸二次优化问题，化成这个样子可以方便优化，前人已经总结了专门一套方法来对付这个问题，只要你能化成这个样子，总是可以求出一个最好的解。不过这里面涉及到最优化的知识，我对这方面了解也不是很深入，所以书里提到通过对偶算法去优化这部分我就不写了。

2. 线性不可分支持向量机

  接下来就看一下线性SVM，用到的知识是软间隔最大化。
  刚才的线性可分SVM并不能处理特异点。特意点就是下图这种，你无法找出一条线把所有样本点准确分出来，就算找到一个很好的超平面，也总有一些奇怪的点在对面。

  所以我们需要对约束条件进行调整，使得它的要求不那么严格。
  刚才我们把函数间隔设为1，理论上所有点的函数间隔都应该大于等于1。但是特异点的函数间隔就是小于1的，所以也很好办，引进一个松弛变量$\xi \ge 0$使得函数间隔加上它会大于1。约束条件就长这个样子。

  引入了松弛变量，目标函数也是稍作修改，变成这个样子.
其中c是可调的惩罚参数。
  所以问题求解变成这个样子。

  这个跟前面的优化过程也是一样的，通过对偶算法去求解。同样的这部分我也不懂，就不写了。
  值得一提的是，跟线性可分SVM可以求出唯一的w跟b不同，线性不可分SVM只能求出唯一的w，但b是不唯一的。

3. 非线性支持向量机

  不是所有的数据都是线性可分的，但是某些时候可以把他们映射到另外一个空间，在这个空间他们是线性可分的。下面就是一个例子。

  但是很多时候这种映射是找不到的，你不知道数据长什么样子，有可能是很高维的，而且杂乱无章没有规律可循。所以非线性SVM用到了核技巧，这样不用显式地定义映射函数。核函数跟映射函数的关系略微复杂，有兴趣的小伙伴可以去看一下书。我在这里只写几个常用的核函数。
  如果我们使用多项式核函数

这里的x跟z为数据集里任意两个样本点
  那么最后得到的决策函数就为

  如果我们使用高斯核函数

  那么最后的决策函数就为

  SVM的理论就说到这里了，可能有些地方没讲清楚，或者有错误，请大家多多见谅，因为我也刚学不久。然后我还跑了网上别人给的一些代码，如果大家有想实践也可以跑一下试试。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.model_selection import  train_test_split
import matplotlib.pyplot as plt
# %matplotlib inline
# data
class SVM:
    def __init__(self, max_iter=100, kernel='linear'):
        self.max_iter = max_iter
        self._kernel = kernel

    def init_args(self, features, labels):
        self.m, self.n = features.shape
        self.X = features
        self.Y = labels
        self.b = 0.0

        # 将Ei保存在一个列表里
        self.alpha = np.ones(self.m)
        self.E = [self._E(i) for i in range(self.m)]
        # 松弛变量
        self.C = 1.0

    def _KKT(self, i):
        y_g = self._g(i) * self.Y[i]
        if self.alpha[i] == 0:
            return y_g >= 1
        elif 0 < self.alpha[i] < self.C:
            return y_g == 1
        else:
            return y_g <= 1

    # g(x)预测值，输入xi（X[i]）
    def _g(self, i):
        r = self.b
        for j in range(self.m):
            r += self.alpha[j] * self.Y[j] * self.kernel(self.X[i], self.X[j])
        return r

    # 核函数
    def kernel(self, x1, x2):
        if self._kernel == 'linear':
            return sum([x1[k] * x2[k] for k in range(self.n)])
        elif self._kernel == 'poly':
            return (sum([x1[k] * x2[k] for k in range(self.n)]) + 1)**2

        return 0

    # E（x）为g(x)对输入x的预测值和y的差
    def _E(self, i):
        return self._g(i) - self.Y[i]

    def _init_alpha(self):
        # 外层循环首先遍历所有满足0<a<C的样本点，检验是否满足KKT
        index_list = [i for i in range(self.m) if 0 < self.alpha[i] < self.C]
        # 否则遍历整个训练集
        non_satisfy_list = [i for i in range(self.m) if i not in index_list]
        index_list.extend(non_satisfy_list)

        for i in index_list:
            if self._KKT(i):
                continue

            E1 = self.E[i]
            # 如果E2是+，选择最小的；如果E2是负的，选择最大的
            if E1 >= 0:
                j = min(range(self.m), key=lambda x: self.E[x])
            else:
                j = max(range(self.m), key=lambda x: self.E[x])
            return i, j

    def _compare(self, _alpha, L, H):
        if _alpha > H:
            return H
        elif _alpha < L:
            return L
        else:
            return _alpha

    def fit(self, features, labels):
        self.init_args(features, labels)

        for t in range(self.max_iter):
            # train
            i1, i2 = self._init_alpha()

            # 边界
            if self.Y[i1] == self.Y[i2]:
                L = max(0, self.alpha[i1] + self.alpha[i2] - self.C)
                H = min(self.C, self.alpha[i1] + self.alpha[i2])
            else:
                L = max(0, self.alpha[i2] - self.alpha[i1])
                H = min(self.C, self.C + self.alpha[i2] - self.alpha[i1])

            E1 = self.E[i1]
            E2 = self.E[i2]
            # eta=K11+K22-2K12
            eta = self.kernel(self.X[i1], self.X[i1]) + self.kernel(
                self.X[i2],
                self.X[i2]) - 2 * self.kernel(self.X[i1], self.X[i2])
            if eta <= 0:
                # print('eta <= 0')
                continue

            alpha2_new_unc = self.alpha[i2] + self.Y[i2] * (
                E1 - E2) / eta  #此处有修改，根据书上应该是E1 - E2，书上130-131页
            alpha2_new = self._compare(alpha2_new_unc, L, H)

            alpha1_new = self.alpha[i1] + self.Y[i1] * self.Y[i2] * (
                self.alpha[i2] - alpha2_new)

            b1_new = -E1 - self.Y[i1] * self.kernel(self.X[i1], self.X[i1]) * (
                alpha1_new - self.alpha[i1]) - self.Y[i2] * self.kernel(
                    self.X[i2],
                    self.X[i1]) * (alpha2_new - self.alpha[i2]) + self.b
            b2_new = -E2 - self.Y[i1] * self.kernel(self.X[i1], self.X[i2]) * (
                alpha1_new - self.alpha[i1]) - self.Y[i2] * self.kernel(
                    self.X[i2],
                    self.X[i2]) * (alpha2_new - self.alpha[i2]) + self.b

            if 0 < alpha1_new < self.C:
                b_new = b1_new
            elif 0 < alpha2_new < self.C:
                b_new = b2_new
            else:
                # 选择中点
                b_new = (b1_new + b2_new) / 2

            # 更新参数
            self.alpha[i1] = alpha1_new
            self.alpha[i2] = alpha2_new
            self.b = b_new

            self.E[i1] = self._E(i1)
            self.E[i2] = self._E(i2)
        return 'train done!'

    def predict(self, data):
        r = self.b
        for i in range(self.m):
            r += self.alpha[i] * self.Y[i] * self.kernel(data, self.X[i])

        return 1 if r > 0 else -1

    def score(self, X_test, y_test):
        right_count = 0
        for i in range(len(X_test)):
            result = self.predict(X_test[i])
            if result == y_test[i]:
                right_count += 1
        return right_count / len(X_test)

    def _weight(self):
        # linear model
        yx = self.Y.reshape(-1, 1) * self.X
        self.w = np.dot(yx.T, self.alpha)
        return self.w

def create_data():
    iris = load_iris()
    df = pd.DataFrame(iris.data, columns=iris.feature_names)
    df['label'] = iris.target
    df.columns = [
        'sepal length', 'sepal width', 'petal length', 'petal width', 'label'
    ]
    data = np.array(df.iloc[:100, [0, 1, -1]])
    for i in range(len(data)):
        if data[i, -1] == 0:
            data[i, -1] = -1
    # print(data)
    return data[:, :2], data[:, -1]
X, y = create_data()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)
svm = SVM(max_iter=200)
svm.fit(X_train, y_train)
result1 = svm.score(X_test, y_test)
print("own svm result:", result1)

from sklearn.svm import SVC
clf = SVC()
clf.fit(X_train, y_train)
result2 = clf.score(X_test, y_test)
print("sklearn svm result:", result2)