《统计学习方法》之支持向量机代码实现

    先回忆下感知机，感知机是一种线性分类器，对于任意的线性可分的数据集，感知机算法总是能够找到一个超平面将数据集正确分割。对于一个线性可分的数据集，存在无数多的分离超平面正确分割数据集，例如下图中的$L_1,L_2$，它们都能正确分割数据集，但是它们间是否存在优劣呢？也就是说在众多分离超平面中是否存在一个最好的分割超平面呢？

    例如上图中的$L_1和L_2$，显然$L_2$比$L_1$更好，因为它离两边的数据集都比较远，这也就意味着即使需预测的数据点和当前已知的数据集的分布略有偏差模型也能正确预测，也就是选择$L_2$能够使模型地泛化性能更好。这就引出的支持向量机的核心思想：在众多分割超平面中选出最优的超平面，使得该超平面到不同类别的数据点的最小距离最大，通俗地讲就是离两边都足够远。
    中间的推导就不再赘述了，有疑问的可参考大佬博客，最后引出需要解决的问题如下：

式子中的$\alpha$为拉格朗日乘子，实现过程中需要特别注意的一点就是判断KKT条件是否成立时，对等式的判断需要在误差范围内进行，比如判断A==0，不能直接写if A == 0，而是应该写为 if abs(A) < eps，eps表示允许的最大误差值。
代码：

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
"""
Created on 2020/2/15 15:50
@author: phil
"""
import numpy as np


def default_kernel_function(x, z):
    # 默认的核函数，即计算两个向量的内积
    x = x.flatten()
    z = z.flatten()
    return np.sum(x * z)


class SVM:
    # 支持向量机
    def __init__(self):
        self.alpha = None  # 拉格朗日乘子
        self.b = 0  # 参数：偏置b
        self.kerner_func = None  # 核函数
        self.N = None  # 样本数量
        self.X = None
        self.y = None
        self.C = None  # 松弛变量的惩罚因子
        self.K = None  # K[i][j]表示x_i,x_j在的内积
        self.E = None  # 使用数组记录Ei值
        self.eps = 0.0001  # 精度，如果误差值在精度内则可认为两者相等

    def gx(self, x):
        # 式子7.104定义的g(x)
        result = 0.0
        for i in range(self.N):
            result += self.alpha[i] * self.y[i] * self.kerner_func(self.X[i], x)
        result += self.b
        return result

    def Ex(self, i):
        # 式子7.105定义的E(i)
        return self.gx(self.X[i]) - self.y[i]

    def choose_alpha1(self):
        # 选择第一个不满足KKT条件的拉格朗日乘子
        for i in range(self.N):
            alphai = self.alpha[i]
            gxi = self.gx(self.X[i])
            yi = self.y[i]
            # 注意这里判断KKT条件时不能使用等号判断，因为浮点运算容易产生误差
            # 也就是书中所说的 KKT条件的检验是在eps的范围内进行的
            # 也就是说如果A瞒住abs(A-B) < eps则可认为A==B
            if abs(alphai) < self.eps and yi * gxi < 1:
                # 违反KKT条件 alpha_i == 0 and y_i * g(x_i) >= 1
                return i
            elif abs(alphai-self.C) < self.eps and yi * gxi > 1:
                # 违反KKT条件 alpha_i = C and y_i*g(x_i) <= 1
                return i
            elif 0 < alphai < self.C and abs(yi*gxi-1) > self.eps:
                # 违法KKT条件 0 < alpha_i < C and y_i*g(x_i) == 1
                return i

        # 如果没有alpha违反KKT条件则返回None
        return None

    def choose_alpha2(self, E1):
        # 选择alpha2使得|E1-E2|最大
        alpha2 = -1
        maxE2 = -1
        maxE1_E2 = -1
        for i in range(len(self.X)):
            E2 = self.E[i]
            if maxE1_E2 < abs(E1 - E2):
                maxE1_E2 = abs(E1 - E2)
                alpha2 = i
                maxE2 = E2
        return alpha2, maxE2

    def fill_K(self):
        # 填充K矩阵，K[i][j]表示kernel_function(X[i], X[j])
        self.K = np.zeros((self.N, self.N))
        for i in range(self.N):
            Xi = self.X[i]
            for j in range(i, self.N):
                self.K[i][j] = self.K[j][i] = self.kerner_func(Xi, self.X[j])

    def fill_E(self):
        # 填充E列表，E[i]表示式子7.105定义的Ei
        self.E = np.zeros((self.N, 1))
        for i in range(self.N):
            self.E[i] = self.Ex(i)

    def fit(self, X, y, C, max_iter=10, kernel_func=default_kernel_function):
        # 默认核函数是线性核函数，也就是范围两个向量的内积
        self.X, self.y, self.C = X, y, C
        self.kerner_func = kernel_func
        self.N = X.shape[0]

        # 初始化拉格朗日乘子
        self.alpha = np.zeros((self.N, 1))

        # 预先计算K矩阵
        self.fill_K()

        # 预先计算E列表
        self.fill_E()

        for i in range(max_iter):
            # 首先选出不满住KKT条件的拉格朗日乘子
            alpha1 = self.choose_alpha1()
            if alpha1 is None:
                # 表示所有变量都满足KKT条件
                break
            E1 = self.E[alpha1]
            # 根据选出了的alpha1选择alpha2
            alpha2, E2 = self.choose_alpha2(E1)

            # 记录当前选出的alpha1和alpha2对应的拉格朗日乘子的值
            alpha1_old = self.alpha[alpha1]
            alpha2_old = self.alpha[alpha2]

            # 计算alpha2对应的边界
            if self.y[alpha1] != self.y[alpha2]:
                L = max(0, alpha2_old - alpha1_old)
                H = min(C, C + alpha2_old - alpha1_old)
            else:
                L = max(0, alpha2_old + alpha1_old - C)
                H = min(C, alpha2_old + alpha1_old)

            # 计算没有边界限制的情况下alpha2的取值，也就是二次函数函数值的最低点对应的自变量的取值
            eta = self.K[alpha1][alpha1] + self.K[alpha2][alpha2] - 2*self.K[alpha1][alpha2]
            alpha2_new_uncut = alpha2_old + self.y[alpha2] * (E1 - E2) / eta

            # 根据最低点计算更新后的alpha2
            if alpha2_new_uncut > H:
                alpha2_new = H
            elif L <= alpha2_new_uncut <= H:
                alpha2_new = alpha2_new_uncut
            else:
                alpha2_new = L
            # 根据更新后的alpha2计算出更新后的alpha1
            alpha1_new = alpha1_old + self.y[alpha1] * self.y[alpha2] * (alpha2_old - alpha2_new)

            # 计算b的更新值
            b1_new = -E1 - self.y[alpha1] * self.K[alpha1][alpha1] * (alpha1_new - alpha1_old) \
                     - self.y[alpha2] * self.K[alpha2][alpha1] * (alpha2_new - alpha1_old) + self.b
            b2_new = -E2 - self.y[alpha1] * self.K[alpha1][alpha2] * (alpha1_new - alpha1_old) \
                     - self.y[alpha2] * self.K[alpha2][alpha2] * (alpha2_new - alpha1_old) + self.b

            if 0 < alpha1_new < C:
                self.b = b1_new
            elif 0 < alpha2_new < C:
                self.b = b2_new
            else:
                self.b = (b1_new + b2_new) / 2

            # 更新拉格朗日乘子
            self.alpha[alpha1] = alpha1_new
            self.alpha[alpha2] = alpha2_new

            # 更新列表E
            self.E[alpha1] = self.Ex(alpha1)
            self.E[alpha2] = self.Ex(alpha2)

            # 只有在使用默认核函数的情况下可以使用这种方式计算出w
            if self.kerner_func == default_kernel_function:
                W = np.zeros_like(self.X[0].T)*1.0
                for j in range(self.N):
                    W += self.alpha[j]*self.y[j]*self.X[j].T
                print("W", W)
                print("b", self.b)

    def predict(self, X):
        preds = []
        for Xi in X:
            gxi = self.gx(Xi)
            if gxi >= 0:
                preds.append(1)
            else:
                preds.append(-1)
        return np.array(preds)


if __name__ == "__main__":
    # 使用书上的简单例子做测试
    X = np.array([[3, 3], [4, 3], [1, 1]])
    y = np.array([[1], [1], [-1]])
    model = SVM()
    model.fit(X, y, 999999)     # 将C设置为一个比较大的数
    print("predict", model.predict(X))

参考链接：

1. https://www.pkudodo.com/2018/12/16/1-8/
2. https://github.com/WenDesi/lihang_book_algorithm