阅读本文需要的背景知识点：拉格朗日乘子法、KKT条件、一丢丢编程知识

一、引言

  前面一节我们介绍了一种分类算法——朴素贝叶斯分类器算法，从概率分布的角度进行分类。下面我们会花几节来介绍另一种在分类问题中有着重要地位的算法——支持向量机¹ （Support Vector Machine/SVM）。
  SVM从基础到复杂可以分成三种分别为线性可分支持向量机（也就是硬间隔支持向量机）、线性支持向量机（软间隔支持向量机）、非线性支持向量机（核函数支持向量机），这一节先来介绍第一种最基础的算法——硬间隔支持向量机算法（Hard-margin Support Vector Machine）。

二、模型介绍

原始模型

  先来看下图，展示了一组线性可分的数据集，其中红点表示 -1、蓝叉表示 1 ，可以看到将该数据集分开的直线有无限多条，那么如何选择一条相对合适的直线呢？如下图中的两条直线 A、B，直观上似乎直线 A 比直线 B 更适合作为决策边界，因为对于直线 B，离某些样本点过于近，虽然依然能正确分类样本点，但当在该样本点附近预测时，会得到完全不同的结果，似乎不符合直觉，其泛化能力更差一些，而反观直线 A，其更具有鲁棒性。

图2-1

  根据上面的推断，我们的目标为找到一条直线（多维时为超平面），使得任意样本点到该超平面的距离的最小值最大。超平面表达式如下：

式2-1

  同时使得每个样本点都在正确的分类下面，即满足下式：

式2-2

  观察下图，可以将问题转化为找到两个超平面 A、C，使得两个超平面之间的距离最大，同时每个样本点分类正确。

图2-2

  不妨假设超平面 A 为 wx + b = 1，超平面 B 为 wx + b = -1（由于 w、b 都可以成比例缩放，所以必然可以找到对应的值使得上式成立），同时分类正确的条件如下：

式2-3

  任意一点到超平面 wx + b = 0 的距离为（推导见四）：

式2-4

  由上式可知超平面 A：wx + b = 1，超平面 B：wx + b = -1 之间的间距为（推导见四）：

式2-5

  那么问题转化为了如下形式：

式2-6

  对向量的取模操作可以转化为向量内积的形式，同时取倒数，将最大值问题变成求最小值的问题，得到最后的模型如下：

式2-7

对偶模型

  对原模型用拉格朗日乘子法转换一下：
（1）使用拉格朗日乘子法，即在后面加上 N 个拉格朗日乘子的条件，得到对应的拉格朗日函数
（2）拉格朗日函数对 w 求偏导数并令其为零
（3）得到 w 的解析解
（4）拉格朗日函数对 b 求偏导数并令其为零
（5）得到新的条件

式2-8

  将上面 w、b 的解析解带回到拉格朗日函数：
（1）拉格朗日函数
（2）展开括号
（3）观察（2）式最后一项，必定为零，再带入 w 的解析解
（4）化简整理得到

式2-9

（1）原始问题
（2）使用拉格朗日乘子法转换
（3）转换为对偶问题，要想使得原始问题与对偶问题等价，需要引入 KKT 条件
（4）带入拉格朗日函数

式2-10

KKT（Karush-Kuhn-Tucker）条件²

  前面在推导原始问题的对偶问题中，使得两个问题等价的条件被称为 KKT 条件（Karush–Kuhn–Tucker conditions）。
  对偶问题需满足的 KKT 条件如下：

式2-11

  由上面 KKT 条件可以看出，对于任意一个样本点，其 λ = 0 或者 y(wx + b) = 1，为前者时，其对最后的函数结果没有影响，即我们可以忽略这些样本点。当为后者时，该样本点在两个最大间隔超平面内（图2-2中的虚线），被称为支持向量，这也是支持向量机名字的由来，最后的模型只与这些支持向量有关。
  经过运用拉格朗日乘子法后，得到了原模型的拉格朗日对偶模型并且需要满足如上所示的 KKT 条件：

式2-12

三、算法步骤

  由上面的模型可以看出这是一个二次规划问题³（Quadratic programming），能够通过一些特殊的方法来求解，例如内点法（Interior Point）等等，后面有机会再来介绍求解方法。下面介绍一种专门用于支持向量机优化问题的算法——序列最小优化算法⁴（Sequential minimal optimization/SMO）。
  SMO 算法的思路在笔者看来有点类似前面讲到的坐标下降算法，既然需要优化的参数过多，那就一次优化一个。在 SMO 中则是一次选择两个待优化变量，其余的当作常量，求解上面的对偶问题的极值，更新变量直至收敛。

算法推导

  注意这里的下标 a、b 即为上面描述的两个待优化变量的下标值，其中下标 b 与偏移量 b 不是同一个值，这里没有合适的字母表示了，所以用了 a、b，切勿与偏移量b混淆
  定义一些符号来方便表述，为后面推导做准备：
（1）两个样本点的点积用 K_ab 来表示，后面一节会介绍 SVM 中的一个很重要的特征——核函数，这里的K即代表的 Kernal
（2）用核函数来表示决策边界
（3）、（4）用 v_a 来表示去除了下标为 a、b 的样本点和偏移量

式3-1

  计算表达式，为后面推导做准备：
（1）计算下面表达式
（2）带入 f(x)
（3）化简消去偏移量 b
（4）将括号展开
（5）移项后得到只包含 λ_a 的结果
（6）同上得到只包含 λ_b 的结果

式3-2

  将多个变量的问题转换为两个变量的问题：
（1）原拉格朗日函数
（2）将上面的结果带入函数中，由于只有下标为 a、b 的两个变量，只需要处理包含这两个变量的项、其他的项统一用常数 C 来表示

式3-3

  定义符号为后面的推导做准备：
（1）原问题的条件
（2）只处理下标为 a、b 的两个变量，其他用常量 c 表示（注意这里为 c，切勿与式3-3中的 C 混淆）
（3）移项后两边同时乘以 y_a
（4）、（5）定义符号简化书写
（6）带入符号后 λ_a 的表达式

式3-4

  将两个变量的问题转换为单变量的问题：
（1）将上面得到的 λ_a 的表达式带入拉格朗日函数，得到只有 λ_b 为变量的函数
（2）拉格朗日函数对 λ_b 求偏导数
（3）注意到 s 的平方必然为 1，化简后得到
（4）合并同类项，整理后得到
（5）另偏导数为零，求得 λ_b 的解
（6）合并同类项，整理后得到

式3-5

  推导变量 λ_b 的更新公式：
（1）γ 的表达式
（2）带入 γ
（3）展开括号并带入 v
（4）合并同类项
（5）化简
（6）预测值与实际值之差记为 E
（7）分母表示为 K
（8）带入 E、K 后得到 λ_b 变量的更新公式

式3-6

  推导变量 λ_a 的更新公式：
（1）由式3-4中的（2）式，更新前更新后该条件不变
（2）两边同时乘以 y_a 后移项
（3）合并同类项

式3-7

  变量 λ_b 的限制条件：
（1）、（2）所有 λ 都需要大于等于零
（3）分情况讨论 λ_b 的限制条件
（4）综合（1）式得到最终变量 λ_b 的限制条件

式3-8

  偏移量 b 的计算：
（1）定义集合 S 为所有 λ 大于零的元素，即所有支持向量点对应的拉格朗日乘子
（2）偏移量 b 的计算公式即 y - wx，选取任意一个支持向量点进行计算
（3）更具鲁棒性的做法为计算所有支持向量点的结果取平均值

式3-9

算法具体步骤

1、初始化拉格朗日乘子向量、权重向量、偏移量
2、初始化误差向量
3、选取两个待优化的拉格朗日乘子参数
4、按照式3-6中的（8）式更新 λ_b
5、检查 λ_b 的上下界，超过上界取上界，低于下界则取下界
6、按照式3-7中的（3）式更新 λ_a
7、重新计算偏移量b、权重向量w和误差向量E
8、重复3～7直至所有拉格朗日乘子满足KKT条件或者代价函数没有明显的变化

待优化参数的选择

  根据 John C. Platt 的论文《Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines》⁵ ，给出了一个启发式的选择参数方法：第一个参数选择不满足 KKT 条件的拉格朗日乘子参数，第二个参数则选择使得｜E_a - E_b｜最大的拉格朗日乘子参数。

四、原理证明

任一点到超平面的距离

  假设 x_0 为任意点，x_t 为 x_0 到超平面 wx + b = 0 的投影：
（1）x_0x_t 表示 x_0 到 x_t 的向量，w 为超平面的法向量，其方向与 x_0x_t 一致，所以夹角为零
（2）cos 0 = 1
（3）计算两个向量的内积，将 x_0x_t 写成向量差的形式
（4）展开括号
（5）由于 x_t 是超平面上的点，所以 wx + b = 0
（6）化简得到
（7）由（2）式与（6）式得
（8）求得任一点到超平面的距离

式4-1

  ｜x_0x_t｜即为点到超平面的距离，既得证。

两个超平面的间距

  两个平行的超平面之间的距离等于点到两个超平面的距离之差的绝对值。根据式2-3 的限制条件，下面分两种情况讨论：
（1）当 wx + b 大于等于 1 时，其内部的绝对值可以去掉
（2）化简后得到结果
（3）当 wx + b 小于等于 -1 时，其内部的绝对值可以去掉
（4）化简后得到结果

式4-2

  综合上面两种情况得到两个超平面的距离公式，既得证。

五、代码实现

使用 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, tol = 1e-6):
        """
        算法来自 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)
            else:
                for idx in range(X.shape[0]):
                    if self.alpha[idx] <= 0:
                        continue
                    numChanged += self.update(idx)
            if examineAll:
                examineAll = False
            elif numChanged == 0:
                examineAll = True
            # 计算代价值
            cost = self.calcCost()
            # 当代价值的变化小于容许的范围时结束算法
            if cost - self.cost <= tol:
                break
            self.cost = cost

    def update(self, idx):
        """
        对下标为 idx 的拉格朗日乘子进行更新
        """
        X = self.X
        y = self.y
        alpha = self.alpha
        # 检查当前拉格朗日乘子是否满足KKT条件，满足条件则直接返回 0
        if self.checkKKT(idx):
            return 0
        if len(alpha[(alpha != 0)]) > 1:
            # 按照｜E1 - E2｜最大的原则寻找第二个待优化的拉格朗日乘子的下标
            jdx = self.selectJdx(idx)
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx):
                return 1
        # 当未更新成功时，遍历不为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] == 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx):
                return 1
        # 当依然没有没有更新成功时，遍历为零的拉格朗日乘子进行更新
        for jdx in range(X.shape[0]):
            if alpha[jdx] != 0:
                continue
            # 对下标为 idx、jdx 的拉格朗日乘子进行更新，当成功更新时直接返回 1
            if self.updateAlpha(idx, jdx):
                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):
        """
        检查下标为 idx 的拉格朗日乘子是否满足 KKT 条件
        1. alpha >= 0
        2. y * f(x) - 1 >= 0
        3. 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 r == 0) or (alpha[idx] == 0 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] + b - y
        return errors

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

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

    def updateAlpha(self, idx, jdx):
        """
        对下标为 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)
        V = self.calcV(idx, jdx)
        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.b = self.calcB()
        # 更新权重向量
        self.w = self.w + y[idx] * (newAlphaIdx - oldAlphaIdx) * X[idx] + y[jdx] * (newAlphaJdx - oldAlphaJdx) * X[jdx]
        # 更新拉格朗日乘子向量
        alpha[idx] = newAlphaIdx
        alpha[jdx] = newAlphaJdx
        # 重新计算误差向量
        self.errors = self.calcE()
        return True

smo = SMO(X, y)
smo.fit()
print("w", smo.w, "b", smo.b)

六、第三方库实现

scikit-learn⁶ 实现

from sklearn.svm import SVC

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

七、动画演示

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

图7-1

  下图展示了通过支持向量机拟合数据后的结果，其中浅红色的区域为预测值为 -1 的部分，浅蓝色的区域则为预测值为 1 的部分。两条虚线分别为 wx + b = 1 和 wx + b = -1，可以看到图中的支持向量一共有三个，分布在两条虚线上，符合前面说的支持向量机的特性。

图7-2

八、思维导图

图8-1

九、参考文献

1. en.wikipedia.org/wiki/Suppor…
2. en.wikipedia.org/wiki/Karush…
3. en.wikipedia.org/wiki/Quadra…
4. en.wikipedia.org/wiki/Sequen…
5. www.microsoft.com/en-us/resea…
6. scikit-learn.org/stable/modu…

完整演示请点击这里

注：本文力求准确并通俗易懂，但由于笔者也是初学者，水平有限，如文中存在错误或遗漏之处，恳请读者通过留言的方式批评指正
本文首发于——AI导图，欢迎关注