浏览本文须要的背景知识点：硬距离反对向量机、松弛变量、一丢丢编程常识

一、引言

  后面一节咱们介绍了一种最根底的反对向量机模型——硬距离反对向量机，该模型能对线性可分的数据集进行分类，但事实中的数据集往往是线性不可分的，那么这一节咱们来介绍反对向量机中的第二种——软距离反对向量机 ¹（Soft-margin Support Vector Machine），来解决下面说的数据集线性不可分的问题。

二、模型介绍

原始模型

  咱们先来看下硬距离反对向量机的原始模型如下：

$$
\begin{array}{}
\underset{w, b}{\max} \frac{1}{2} w^Tw \\
\text {s.t.} \quad y_{i}\left(w^{T} x_{i}+b\right) \geq 1 \quad i=1,2, \cdots, N
\end{array}
$$

<center> 式 2 -1</center>

  约束条件 y(wx + b) 大于等于 1 是为了使得所有样本点都在正确的分类下，这也是为什么称为硬距离的起因。而当初数据集无奈用一个超平面齐全离开，这时就须要容许局部数据集不满足上述约束条件。

  批改下面的代价函数，加上分类谬误的样本点的惩办项，其中 1(x) 在后面章节中提到过，即批示函数 ²（indicator function），当满足上面不等式时函数返回 1，不满足时函数返回 0。同时通过常数 C 来管制惩办力度，其中 C > 0，能够看到当 C 取无限大时，会迫使每一个样本点分类正确，等同于硬距离反对向量机。

$$
\begin{array}{}
\underset{w, b}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} 1_{x_{i}}\left(y_{i}\left(w^{T} x_{i}+b\right) \lt 1\right) \\
\text {s.t.} \quad C>0
\end{array}
$$

<center> 式 2 -2</center>

  式 2 - 2 中的批示函数既不是凸函数又不是连续函数，解决起来比拟麻烦，这时能够用 max 函数来替换，模式如下：

$$
\begin{array}{}
\underset{w, b}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \max\left(0, 1 – y_{i}\left(w^{T} x_{i}+b\right)\right) \\
\text {s.t.} \quad C>0
\end{array}
$$

<center> 式 2 -3</center>

  这时再引入松弛变量 ³（slack variable）ξ，同时加上如下的条件，进一步将 max 函数去掉，式 2 - 3 就是软距离反对向量机的原始模型：

$$
\begin{array}{}
\underset{w, b, \xi}{\min} \frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i} \\
\text {s.t.} \quad C>0 \quad \xi_{i} \geq 0 \quad \xi_{i} \geq 1-y_{i}\left(w^{T} x_{i}+b\right) \quad i=1,2, \cdots, N
\end{array}
$$

<center> 式 2 -4</center>

对偶模型

  与硬距离反对向量机的解决形式一样，同样对原始模型用拉格朗日乘子法进行转换：

（1）应用拉格朗日乘子法，引入两类拉格朗日参数 λ、μ，失去拉格朗日函数

（2）拉格朗日函数对 w 求偏导数并令其为零

（3）同硬距离反对向量机一样失去 w 的解析解

（4）拉格朗日函数对 b 求偏导数并令其为零

（5）拉格朗日函数对 ξ 求偏导数并令其为零

（6）失去 C = λ + μ

$$
\begin{aligned}
L(w, b, \xi, \lambda, \mu) &=\frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}\left(1-\xi_{i}-y_{i}\left(w^{T} x_{i}+b\right)\right)-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (1)\\
\frac{\partial L}{\partial w} &=w-\sum_{i=1}^{N} \lambda_{i} y_{i} x_{i}=0 & (2)\\
w &=\sum_{i=1}^{N} \lambda_{i} y_{i} x_{i} & (3)\\
\frac{\partial L}{\partial b} &=\sum_{i=1}^{N} \lambda_{i} y_{i}=0 & (4)\\
\frac{\partial L}{\partial \xi_{i}} &=C-\lambda_{i}-\mu_{i}=0 & (5)\\
C &=\lambda_{i}+\mu_{i} & (6)
\end{aligned}
$$

<center> 式 2 -5</center>

  将式 2 - 5 失去的后果带回拉格朗日函数中：

（1）拉格朗日函数

（2）带入后开展括号

（3）能够看到（2）式中第 2、5 项互相对消，第 3、8 项互相对消，第 7 项为零，合并后得

$$
\begin{aligned}
L(w, b, \xi, \lambda, \mu) &=\frac{1}{2} w^{T} w+C \sum_{i=1}^{N} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}\left(1-\xi_{i}-y_{i}\left(w^{T} x_{i}+b\right)\right)-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (1)\\
&=\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i}+\sum_{i=1}^{N} \lambda_{i} \xi_{i}+\sum_{i=1}^{N} \mu_{i} \xi_{i}+\sum_{i=1}^{N} \lambda_{i}-\sum_{i=1}^{N} \lambda_{i} \xi_{i}-\sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i}-\sum_{i=1}^{N} \lambda_{i} y_{i} b-\sum_{i=1}^{N} \mu_{i} \xi_{i} & (2)\\
&=\sum_{i=1}^{N} \lambda_{i}-\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{i} & (3)
\end{aligned}
$$

<center> 式 2 -6</center>

  在来看下拉格朗日乘子参数的条件：

（1）拉格朗日乘子 μ 的限度条件

（2）两边同时加上拉格朗日乘子 λ

（3）联合式 2 -5（6）的后果失去

$$
\begin{aligned}
\mu_i &\ge 0 & (1) \\
\lambda_i + \mu_i &\ge \lambda_i & (2) \\
C &\ge \lambda_i & (3)
\end{aligned}
$$

<center> 式 2 -7</center>

  同硬距离反对向量机一样，引入 KKT 条件如下：

$$
\left\{\begin{aligned}
\nabla_{w, b, \xi} L(w, b, \xi, \lambda, \mu) &=0 \\
\lambda_{i} & \geq 0 \\
y_{i}\left(w^{T} x_{i}+b\right)-1+\xi_{i} & \geq 0 \\
\lambda_{i}\left(y_{i}\left(w^{T} x_{i}+b\right)-1+\xi_{i}\right) &=0 \\
\mu_{i} & \geq 0 \\
\xi_{i} & \geq 0 \\
\mu_{i} \xi_{i} &=0
\end{aligned}\right.
$$

<center> 式 2 -8</center>

  通过使用拉格朗日乘子法后，失去了原模型的拉格朗日对偶模型并且须要满足如上所示的 KKT 条件：

$$
\begin{array}{}
\underset{\lambda}{\max} \sum_{i=1}^{N} \lambda_{i}-\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \lambda_{i} \lambda_{j} y_{i} y_{j} x_{i}^{T} x_{j} \\
\text {s.t.} \quad \sum_{i=1}^{N} \lambda_{i} y_{i}=0 \quad 0 \leq \lambda_{i} \leq C \quad i=1,2, \cdots, N
\end{array}
$$

<center> 式 2 -9</center>

三、算法步骤

  将下面式 2 - 9 的对偶模型与硬距离反对向量机的对偶模型互相比拟后，会发现两者的差异就在于拉格朗日乘子的约束条件不同，前者只多了一个惩办因子 C 的上界，所以其求解的算法与硬距离反对向量机的大致相同，只有如下两个中央不一样。

  算法中一个有改变的中央是新 λ_b 的上下界的计算：

（1）、（2）所有 λ 都须要大于等于零并且小于等于 C

（3）分状况探讨 λ_b 的限度条件

（4）综合（1）式失去最终变量 λ_b 的限度条件

$$
\begin{aligned}
&0 \leq \lambda_{b}^{\text {new}} \leq C & (1)\\
&0 \leq \lambda_{a}^{\text {old}}+y_{a} y_{b}\left(\lambda_{b}^{\text {old}}-\lambda_{b}^{\text {new}}\right) \leq C & (2)\\
\Rightarrow&\left\{\begin{array}{}
\lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}-C \leq \lambda_{b}^{\text {new}} \leq \lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}, \quad y_{a} y_{b}=+1 \\
\lambda_{b}^{\text {old}}-\lambda_{a}^{\text {old}} \leq \lambda_{b}^{\text {new}} \leq C-\lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}, \quad y_{a} y_{b}=-1
\end{array}\right. & (3)\\
\Rightarrow&\left\{\begin{array}{}
\max \left(0, \lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}-C\right) \leq \lambda_{b}^{\text {new}} \leq \min \left(C, \lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}\right), \quad y_{a} y_{b}=+1 \\
\max \left(0, \lambda_{b}^{\text {old}}-\lambda_{a}^{\text {old}}\right) \leq \lambda_{b}^{\text {new}} \leq \min \left(C, C-\lambda_{a}^{\text {old}}+\lambda_{b}^{\text {old}}\right), \quad y_{a} y_{b}=-1
\end{array}\right. & (4)
\end{aligned}
$$

<center> 式 3 -1</center>

  算法中另一个有改变的中央是 KKT 条件的判断：

（1）思考 λ 等于 0 的状况

（2）依据式 2 -5 中（6）的后果可知 μ 等于 C

（3）因为 μ 等于 C，依据式 2 - 8 的 KKT 条件可知 ξ 等于 0

（4）将 ξ 等于 0 带入 KKT 条件得

$$
\begin{aligned}
\lambda &=0 & (1)\\
\mu &=C & (2)\\
\xi &=0 & (3)\\
y_{i}\left(w^{T} x_{i}+b\right)-1 & \geq 0 & (4)
\end{aligned}
$$

<center> 式 3 -2</center>

（1）思考 λ 等于 C 的状况

（2）依据式 2 -5 中（6）的后果可知 μ 等于 0

（3）因为 λ 不等于 0，依据式 2 - 8 的 KKT 条件可知

（4）因为 ξ 大于等于 0 联合（3）式可得

$$
\begin{aligned}
\lambda &=C & (1)\\
\mu &=0 & (2)\\
y_{i}\left(w^{T} x_{i}+b\right)-1 + \xi &=0 & (3)\\
y_{i}\left(w^{T} x_{i}+b\right)-1 & \leq 0 & (4)
\end{aligned}
$$

<center> 式 3 -3</center>

（1）思考 λ 介于 0 与 C 之间的状况

（2）依据式 2 -5 中（6）的后果可知 μ 也介于 0 与 C 之间

（3）依据式 2 - 8 的 KKT 条件可知

（4）因为 λ 不等于 0，依据式 2 - 8 的 KKT 条件并联合（3）式可得

$$
\begin{aligned}
0 \lt \lambda &\lt C & (1)\\
0 \lt \mu &\lt C & (2)\\
\xi &=0 & (3)\\
y_{i}\left(w^{T} x_{i}+b\right)-1 & = 0 & (4)
\end{aligned}
$$

<center> 式 3 -4</center>

  综合下面式 3 -2、式 3 -3、式 3 - 4 可得如下的新限度条件：

$$
\left\{\begin{array}{}
y_{i}\left(w^{T} x_{i}+b\right)-1 \geq 0, & \lambda=0 \\
y_{i}\left(w^{T} x_{i}+b\right)-1=0, & 0 \lt \lambda \lt C \\
y_{i}\left(w^{T} x_{i}+b\right)-1 \leq 0, & \lambda=C
\end{array}\right.
$$

<center> 式 3 -5</center>

  算法步骤请参考硬距离反对向量机的第三小节和上面的代码实现。

四、代码实现

应用 Python 实现

import numpy as np

class SMO:
    """
    软距离反对向量机
    序列最小优化算法实现（Sequential minimal optimization/SMO）"""

    def __init__(self, X, y):
        # 训练样本特色矩阵（N * p）self.X = X
        # 训练样本标签向量（N * 1）self.y = y
        # 拉格朗日乘子向量（N * 1）self.alpha = np.zeros(X.shape[0])
        # 误差向量，默认为负的标签向量（N * 1）self.errors = -y
        # 偏移量 
        self.b = 0
        # 权重向量（p * 1）self.w = np.zeros(X.shape[1])
        # 代价值
        self.cost = -np.inf

    def fit(self, C = 1, tol = 1e-4):
        """
        算法来自 John C. Platt 的论文
        https://www.microsoft.com/en-us/research/uploads/prod/1998/04/sequential-minimal-optimization.pdf
        """
        # 更新变动次数
        numChanged = 0
        # 是否查看全副
        examineAll = True
        while numChanged > 0 or examineAll:
            numChanged = 0
            if examineAll:
                for idx in range(X.shape[0]):
                    numChanged += self.update(idx, C)
            else:
                for idx in range(X.shape[0]):
                    if self.alpha[idx] <= 0:
                        continue
                    numChanged += self.update(idx, C)
            if examineAll:
                examineAll = False
            elif numChanged == 0:
                examineAll = True
            # 计算代价值
            cost = self.calcCost()
            if self.cost > 0:
                # 当代价值的变动小于答应的范畴时完结算法
                rate = (cost - self.cost) / self.cost
                if rate <= tol:
                    break
            self.cost = cost

    def update(self, idx, C = 1):
        """对下标为 idx 的拉格朗日乘子进行更新"""
        X = self.X
        y = self.y
        alpha = self.alpha
        # 查看以后拉格朗日乘子是否满足 KKT 条件，满足条件则间接返回 0
        if self.checkKKT(idx, C):
            return 0
        if len(alpha[(alpha != 0)]) > 1:
            # 依照|E1 - E2|最大的准则寻找第二个待优化的拉格朗日乘子的下标
            jdx = self.selectJdx(idx)
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当胜利更新时间接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 当未更新胜利时，遍历不为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] == 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当胜利更新时间接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 当仍然没有没有更新胜利时，遍历为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] != 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当胜利更新时间接返回 1
            if self.updateAlpha(idx, jdx, C):
                return 1
        # 仍然没有更新时返回 0
        return 0

    def selectJdx(self, idx):
        """寻找第二个待优化的拉格朗日乘子的下标"""
        errors = self.errors
        if errors[idx] > 0:
            # 当误差项大于零时，抉择误差向量中最小值的下标
            return np.argmin(errors)
        elif errors[idx] < 0:
            # 当误差项小于零时，抉择误差向量中最大值的下标
            return np.argmax(errors)
        else:
            # 当误差项等于零时，抉择误差向量中最大值和最小值的绝对值最大的下标
            minJdx = np.argmin(errors)
            maxJdx = np.argmax(errors)
            if max(np.abs(errors[minJdx]), np.abs(errors[maxJdx])) - errors[minJdx]:
                return minJdx
            else:
                return maxJdx

    def calcB(self):
        """
        计算偏移量
        别离计算每一个拉格朗日乘子不为零对应的偏移量后取其平均值
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_gt = alpha[alpha > 0]
        # 拉格朗日乘子向量中不为零的数量
        alpha_gt_len = len(alpha_gt)
        # 全副为零时间接返回 0
        if alpha_gt_len == 0:
            return 0
        # b = y - Wx，具体算法请参考文章中的阐明
        X_gt = X[alpha > 0]
        y_gt = y[alpha > 0]
        alpha_gt_y = np.array(np.multiply(alpha_gt, y_gt)).reshape(-1, 1)
        s = np.sum(np.multiply(alpha_gt_y, X_gt), axis=0)
        return np.sum(y_gt - X_gt.dot(s)) / alpha_gt_len

    def calcCost(self):
        """
        计算代价值
        依照文章中的算法计算即可
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        cost = 0
        for idx in range(X.shape[0]):
            for jdx in range(X.shape[0]):
                cost = cost + (y[idx] * y[jdx] * X[idx].dot(X[jdx]) * alpha[idx] * alpha[jdx])
        return np.sum(alpha) - cost / 2

    def checkKKT(self, idx, C = 1):
        """
        查看下标为 idx 的拉格朗日乘子是否满足 KKT 条件
        1. alpha >= 0
        2. alpha <= C
        3. y * f(x) - 1 >= 0
        4. alpha * (y * f(x) - 1) = 0
        """
        y = self.y
        errors = self.errors
        alpha = self.alpha
        r = errors[idx] * y[idx]
        if (alpha[idx] > 0 and alpha[idx] < C and r == 0) or (alpha[idx] == 0 and r > 0) or (alpha[idx] == C and r < 0):
            return True
        return False

    def calcE(self):
        """
        计算误差向量
        E = f(x) - y
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        alpha_y = np.array(np.multiply(alpha, y)).reshape(-1, 1)
        errors = X.dot(X.T).dot(alpha_y).T[0] + self.b - y
        return errors

    def calcU(self, idx, jdx, C = 1):
        """
        计算拉格朗日乘子上界，使两个待优化的拉格朗日乘子同时大于等于 0
        依照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return max(0, alpha[jdx] + alpha[idx] - C)
        else:
            return max(0.0, alpha[jdx] - alpha[idx])

    def calcV(self, idx, jdx, C = 1):
        """
        计算拉格朗日乘子下界，使两个待优化的拉格朗日乘子同时大于等于 0
        依照文章中的算法计算即可
        """
        y = self.y
        alpha = self.alpha
        if y[idx] * y[jdx] == 1:
            return min(C, alpha[jdx] + alpha[idx])
        else:
            return min(C, C + alpha[jdx] - alpha[idx])

    def updateAlpha(self, idx, jdx, C = 1):
        """
        对下标为 idx、jdx 的拉格朗日乘子进行更新
        依照文章中的算法计算即可
        """
        if idx == jdx:
            return False
        X = self.X
        y = self.y
        alpha = self.alpha
        errors = self.errors
        # idx 的误差项
        Ei = errors[idx]
        # jdx 的误差项
        Ej = errors[jdx]
        Kii = X[idx].dot(X[idx])
        Kjj = X[jdx].dot(X[jdx])
        Kij = X[idx].dot(X[jdx])
        # 计算 K
        K = Kii + Kjj - 2 * Kij
        oldAlphaIdx = alpha[idx]
        oldAlphaJdx = alpha[jdx]
        # 计算 jdx 的新拉格朗日乘子
        newAlphaJdx = oldAlphaJdx + y[jdx] * (Ei - Ej) / K
        U = self.calcU(idx, jdx, C)
        V = self.calcV(idx, jdx, C)
        if newAlphaJdx < U:
            # 当新值超过上界时，批改其为上界
            newAlphaJdx = U
        if newAlphaJdx > V:
            # 当新值低于下界时，批改其为下界
            newAlphaJdx = V
        if oldAlphaJdx == newAlphaJdx:
            # 当新值与旧值相等时，判断为未更新，间接返回
            return False
        # 计算 idx 的新拉格朗日乘子
        newAlphaIdx = oldAlphaIdx + y[idx] * y[jdx] * (oldAlphaJdx - newAlphaJdx)
        # 更新权重向量
        self.w = self.w + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx] + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx]
        # 更新拉格朗日乘子向量
        alpha[idx] = newAlphaIdx
        alpha[jdx] = newAlphaJdx
        oldB = self.b
        # 从新计算偏移量
        self.b = self.calcB()
        # 从新计算误差向量
        eps = np.finfo(np.float32).eps
        newErrors = []
        for i in range(X.shape[0]):
            newError = errors[i] + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx].dot(X[i]) + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx].dot(X[i]) - oldB + self.b
            if np.abs(newError) < eps:
                newError = 0
            newErrors.append(newError)
        self.errors = newErrors
        return True

五、第三方库实现

scikit-learn⁴ 实现

from sklearn.svm import SVC

svc = SVC(kernel = "linear", C = 1)
# 拟合数据
svc.fit(X, y)
# 权重系数
w = svc.coef_
# 截距
b = svc.intercept_
print("w", w, "b", b)

六、动画演示

  下图展现了一个线性不可分的数据集，其中红色示意标签值为 -1 的样本点，蓝色代表标签值为 1 的样本点：

<center> 图 6 -1</center>

  下图展现了通过软距离反对向量机拟合数据后的后果，其中浅红色的区域为预测值为 -1 的局部，浅蓝色的区域则为预测值为 1 的局部：

<center> 图 6 -2</center>

七、思维导图

<center> 图 7 -1</center>

八、参考文献

1. https://en.wikipedia.org/wiki…
2. https://en.wikipedia.org/wiki…
3. https://en.wikipedia.org/wiki…
4. https://scikit-learn.org/stab…

残缺演示请点击这里

注：本文力求精确并通俗易懂，但因为笔者也是初学者，程度无限，如文中存在谬误或脱漏之处，恳请读者通过留言的形式批评指正

本文首发于——AI 导图 ，欢送关注