支持向量机原理简介
Kotori Y 27 Posts
 2022-08-26 08:27 2022-08-30 19:12      

SVM有三宝：间隔、对偶、核技巧

在机器学习中，支持向量机（Support Vector Machine, SVM）是在分类与回归分析中分析数据的监督式学习模型与相关的学习算法。给定一组训练实例，每个训练实例被标记为属于两个类别中的一个或另一个，SVM训练算法创建一个将新的实例分配给两个类别之一的模型，使其成为非概率二元线性分类器。SVM模型是将实例表示为空间中的点，这样映射就使得单独类别的实例被尽可能宽的明显的间隔分开。然后，将新的实例映射到同一空间，并基于它们落在间隔的哪一侧来预测所属类别。

1. 预备

1.1 点到超平面的距离

在二维空间中的一点到直线的距离为：

将式扩展到样本空间中的任意点到超平面的距离：

其中为法向量，决定了超平面的法向量；为位移项，决超平面与原点之间的距离。

1.2 不等式约束与KKT条件

对于给定的目标函数与约束条件，求函数的极值问题可以构建一个拉格朗日函数:

分别对式中的, , 求偏导，并令结果等于0，可解出此方程。

而对于带有不等式约束，同样可以构建拉格朗日函数，并为两种情况讨论

  1. 原函数的最值在约束范围内，此时约束条件并没有起到作用，满足条件：

  2. 原函数的最值在约束范围外，此时的最值一定在约束条件的边界，满足条件：

合并上述情况的条件，即KKT（Karush-Kuhn-Tucker）条件：

1.3 凸二次规划问题

凸函数的定义为:

仿射函数为最高次数为1的多项式函数，仿射函数是一种特殊的凸函数:

凸优化问题是一种特殊的约束最优化问题，其中要求目标函数和不等式约束是可微的凸函数，等式约束是可微的仿射函数。当目标函数是二次型函数:

且不等式约束函数是仿射函数时，就变成一个凸二次规划问题。

1.4 矩阵（向量）求导的几个公式

2. 最大间隔与支持向量

对于给定的训练样本集，。假设超平面能将训练样本进行正确的分类，那么对于样本将满足以下条件：

式中的两个向量与被称为“支持向量”，两个支持向量到超平面的距离之和为：

它被称为“间隔”（margin）。

SVM的主要原理就是寻找具有“最大间隔”（maximum margin）的划分超平面，也就是要找到能满足中的约束参数和，使得最大：

式等效为：

这就是SVM的基本型。

3. 最大间隔的对偶问题

可以观察到式是一个凸二次规划问题，虽然可以使用一些软件包求解，但存在更高效的方法。首先，构造式的拉格朗日函数：

因为的拉格朗日函数一定有：

当拉格朗日函数越接近目标函数，最小值越有可能出现，那么原始问题等效为：

为了方便求解，将转换为其对偶问题（dual problem）:

首先求，式对求偏导：

令式的结果为0，可得出：

式对求偏导：

令式的结果为0：

将和带入式：

对偶问题为：

若能解出，即可得到超平面方程：

另外上述条件将符合KKT条件：

可以观察到，训练样本中必有或。符合前者的训练样本将不会体现在中，最终的模型不必保留这些点。而符合后者的样本将位于最大间隔边界上。

4. 软间隔支持向量机

上述中的间隔被称为“硬间隔”（hard margin），而实际中的训练数据不可避免的存在噪声，对于这种情况需要引入“软间隔”（soft margin）的概念，即允许一定的错误，即不满足：

同时，不满足约束的样本应尽可能的少，在原优化目标基础上加上损失函数：

其中是一个常数，是“0/1损失函数”：

另外还有hinge损失，式变为：

引入”松弛变量“（slack variables），将重写为：

构建的拉格朗日函数:

原问题为：

对偶问题为：

式对，，求偏导并令结果为零，可得：

将对-代入对偶问题：

相较于唯一的不同是约束条件的差异，软间隔支持向量机的KKT条件为：

5. 非线性支持向量机与核函数

在之前的章节中，我们假设样本在它们原来的空间中是线性可分的，但对于真实世界的样本往往无法找到一个这样的超平面将两类样本正确的分开。幸运的是，如果原始空间是有限维，及属性数有限，那么一定存在一个高维特征空间使得样本线性可分。

令表示将映射后的特征向量，超平面方程可表示为：

类似

对偶问题为：

假设从存在这样一个函数：

通过上面式可以避免直接计算高维空间中的内积，函数就是核函数（kernel function）。

式重写为：

求解后可得：

常用的核函数有：

名称 表达式 参数
线性核
多项式核 为多项式的次数
高斯核（RBF核） 为高斯核的带宽
拉普拉斯核
Sigmoid核 为双曲正切函数，，

6. SMO求解拉格朗日参数

6.1 权重的更新

SMO的思想类似坐标上升算法，我们需要优化一系列的的值，我们每次选择尽量少的来优化，不断迭代直到函数收敛到最优值。
式加入核技巧后：

SMO的核心思想是，对于一系列待优化的参数，先选取两个参数，固定剩下的参数进行优化。

选择式中的与，将其他视为常数项：

优化目标及约束：

为了叙述简洁，记：

将式代入的目标函数：

由式的约束条件可得：

将式代入:

这样就变成了只含有的函数，式对求导并令结果等于0：

将-代入得：

令：

将代入:

的取值应该符合的约束条件。

  • 如果，有:

  • 如果：

因此：

并且：

6.2 偏置的更新

  • 当：

    因此：

    有式得：

    将代入:

  • 类似地，当：

  • 其他情况：

7. 代码简易实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
import numpy as np
from numpy.typing import ArrayLike


def linear_kernel(xi: np.matrix, xj: np.matrix):
    assert xi.shape[0] == xj.shape[0]
    assert xi.shape[1] == xj.shape[1] == 1

    return xi.T * xj


class SVMClassifier:

    def __init__(self, gamma: float, kernel: str):
        self._C = np.mat([[gamma]])
        self._kernel = kernel
        self._alphas = None
        self._b = None
        self._w = None

    @property
    def C(self):
        return self._C

    @property
    def kernel(self):
        return self._kernel

    @property
    def alphas(self):
        return self._alphas

    @property
    def b(self):
        return self._b

    @property
    def w(self):
        return self._w

    def random_choose_index_j(self, i, m):
        j = np.random.randint(0, m)
        if j == i:
            return self.random_choose_index_j(i, m)
        return j

    def calculate_low_and_high(self, alpha_i_old, alpha_j_old, y_i, y_j):

        if y_i != y_j:
            low = max(np.mat([[0]]), alpha_j_old - alpha_i_old)
            high = min(self.C, self.C + alpha_j_old - alpha_i_old)

            return low, high

        low = max(np.mat([[0]]), alpha_j_old + alpha_i_old - self.C)
        high = min(self.C, alpha_j_old + alpha_i_old)

        return low, high

    def kernel_func(self, *args):
        if self.kernel == "linear":
            return linear_kernel(*args)

    def _gx(self, X_mat, y_mat, i):
        return np.multiply(self.alphas, y_mat).T * (X_mat * X_mat[i, :].T) + self.b

    def calculate_eta(self, X_i, X_j):

        k_ii = self.kernel_func(X_i, X_i)
        k_jj = linear_kernel(X_j, X_j)
        k_ij = linear_kernel(X_i, X_j)

        return k_ii + k_jj - 2 * k_ij

    def clip_alpha(self, alpha, low, high):
        if low <= alpha <= high:
            return alpha

        if alpha < low:
            return low

        return high

    def fit(self, X_train: ArrayLike, y_train: ArrayLike, epochs):
        m, n = X_train.shape

        self._alphas = np.mat(np.zeros((m, 1)))
        self._b = np.mat([[0]])

        X_mat = np.mat(X_train)
        y_mat = np.mat(y_train.reshape(-1, 1))

        for epoch in range(epochs):
            for i in range(m):
                gx_i = self._gx(X_mat, y_mat, i)
                ex_i = gx_i - y_mat[i, :]

                j = self.random_choose_index_j(i, m)
                gx_j = self._gx(X_mat, y_mat, j)
                ex_j = gx_j - y_mat[j, :]

                y_i = y_mat[i, :]
                y_j = y_mat[j, :]
                X_i = X_mat[i, :].T
                X_j = X_mat[j, :].T
                alpha_i_old = self.alphas[i, :]
                alpha_j_old = self.alphas[j, :]

                eta = self.calculate_eta(X_i, X_j)
                if eta <= 0:
                    continue

                _alpha_j = alpha_j_old + (y_j * (ex_i - ex_j)) / eta

                low, high = self.calculate_low_and_high(alpha_i_old, alpha_j_old, y_i, y_j)

                alpha_j_new = self.clip_alpha(_alpha_j, low=low, high=high)
                alpha_i_new = alpha_i_old + y_i * y_j * (alpha_j_old - alpha_j_new)

                b1 = - ex_i - y_i * linear_kernel(X_i, X_i) * (alpha_i_new - alpha_i_old) - \
                     y_j * linear_kernel(X_j, X_i) * (alpha_j_new - alpha_j_old) + self.b

                b2 = - ex_j - y_i * linear_kernel(X_i, X_j) * (alpha_i_new - alpha_i_old) - \
                     y_j * linear_kernel(X_j, X_j) * (alpha_j_new - alpha_j_old) + self.b

                if np.mat([[0]]) < alpha_i_new < self.C:
                    self._b = b1
                elif np.mat([[0]]) < alpha_j_new < self.C:
                    self._b = b2
                else:
                    self._b = (b1 + b2) / 2

                self._alphas[i, 0] = alpha_i_new
                self._alphas[j, 0] = alpha_j_new

        self._w = np.multiply(self.alphas, y_mat).T * X_mat

在鸢尾花数据集上进行测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

iris = load_iris()
_X_train, _y_train = iris.data[:, :2], iris.target
_X_train = _X_train[_y_train != 2]
_y_train = _y_train[_y_train != 2]
_y_train = np.where(_y_train == 0, -1, 1)

model = SVMClassifier(gamma=10, kernel="linear")
model.fit(_X_train, _y_train, 500)

fig, ax = plt.subplots()
ax.scatter(*_X_train[_y_train == 1].T, s=25)
ax.scatter(*_X_train[_y_train == -1].T, s=25)

w_0, w_1 = np.array(model.w).flatten()
# w_0 * X_0 + w_1 * X_1 + b = 0
# => w_1 * X_1 = - (b +  w_0 * X_0)
# => X_1 = - (b +  w_0 * X_0) / w_1
x_0 = np.arange(4, 7, 0.01)
x_1 = - (model.b[0, 0] + w_0 * x_0) / w_1

# support vec
x_1_p = (- (model.b[0, 0] + w_0 * x_0) + 1) / w_1
x_1_n = (- (model.b[0, 0] + w_0 * x_0) - 1) / w_1

ax.plot(x_0, x_1, color='k', lw=1)
ax.plot(x_0, x_1_p, color='k', lw=1, ls='--')
ax.plot(x_0, x_1_n, color='k', lw=1, ls='--')

support = _X_train[np.array(model.alphas != 0).flatten()]
ax.scatter(*support.T, s=100, color='m', zorder=-1, label='support vector')

ax.legend()
plt.show()

8. 参考资料

  1. 机器学习 - 周志华
  2. 统计学习方法 - 李航
  • Post title：支持向量机原理简介
  • Post author：Kotori Y
  • Create time：2022-08-26 08:27
  • Update time：2022-08-30 19:12
  • Post link：https://blog.iamkotori.com/2022/08/26/支持向量机原理简介/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments