机器学习｜SVM支持向量机代码实现解读

作者: 蓝白绛 | 来源:发表于2019-02-20 19:39 被阅读0次

前言

如果你能找到这里，真是我的幸运~这里是蓝白绛的学习笔记，本集合主要针对《统计学习方法》这本书，主要是基本机器学习算法代码实现的解读。
本篇的代码整理来自github：wepe-Basic Machine Learning and Deep Learning
本篇整理SVM支持向量机实现，github地址为：wepe-SVM

第七章支持向量机

1、模型

支持向量机(Support Vector Machine，SVM)是一种二类分类模型。它的基本模型是定义在特征空间上的间隔最大的线性分类器；支持向量机还包括核技巧，使它成为实质上的非线性分类器。
支持向量机学习方法包含构建由简至繁的模型：线性可分支持向量机(linear support vector machine in linearly separable case)、线性支持向量机(linear support vector machine)、非线性可分支持向量机(non-linear support vector machine)。

线性可分支持向量机：当训练数据线性可分时，通过硬间隔最大化，学习一个线性分类器，即硬间隔支持向量机。
线性支持向量机：当训练数据近似线性可分时，引入松弛变量 $\xi$ ，通过软间隔最大化，学习一个线性分类器，即软间隔支持向量机。
非线性支持向量机：当训练数据线性不可分时，通过使用核技巧 $K(x,z)=\phi(x)\cdot\phi(z)$ 及软间隔最大化，学习一个非线性分类器。

函数间隔：超平面 $(w,b)$ 关于样本点 $(x_i,y_i)$ 的函数间隔为 $\hat{\gamma}_i=y_i(w\cdot x_i+b)$ 超平面 $(w,b)$ 关于训练数据集 $T$ 的函数间隔为 $\hat{\gamma}=\min_{i=1,2...N}\hat{y}_i.$
几何间隔：超平面 $(w,b)$ 关于样本点 $(x_i,y_i)$ 的几何间隔为 $\gamma_i=\frac{y_i(w\cdot x_i+b)}{\|w\|}$ 超平面 $(w,b)$ 关于训练数据集 $T$ 的函数间隔为 $\gamma=\min_{i=1,2...N}y_i.$
拿最简单的线性可分支持向量机来说，原始问题为： $\begin{align} &\max_{w,b}\quad \frac{1}{2}\|w\|^2\\ & \begin{array}{r@{\quad}r@{}l@{\quad}l} s.t.&y_i(w\cdot x_i+b)-1 &\geq0 &i=1,2,...,N \end{array} \end{align}$ 经过构建拉格朗日函数，对 $w$ 和 $b$ 取偏导，令其等于 $0$ 并带回拉格朗日函数，转化为对偶问题为： $\begin{aligned} &\min_{\alpha}\quad \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i x_j-\sum_{i=1}^N\alpha_i \\ & \begin{array}{r@{\quad}r@{}l@{\quad}l} s.t.& \sum\limits_{i=1}^N\alpha_iy_i & =0 \\ & \alpha_i & \geq0 & i=1,2,...,N \end{array} \end{aligned}$
KKT条件为： $\begin{aligned} & \nabla_wL(w^*,b^*,\alpha^*)=w-\sum_{i=1}^N\alpha_i^*y_ix_i=0 \\ & \nabla_bL(w^*,b^*,\alpha^*)=\sum_{i=1}^N\alpha^*_iy_i=0 \\ & \alpha_i^*(y_i(w^*\cdot x_i+b^*)-1)=0,\quad i=1,2,...,N \\ & y_i(w^*\cdot x_i+b^*)-1\geq0,\quad i=1,2,...,N \\ & \alpha_i^*\geq0,\quad i=1,2,...,N \end{aligned}$ 最终可求得 $w^*=\sum_{i=1}^N\alpha_i^*y_ix_i$ 选取其中一个 $\alpha_j^*>0$ 可求得 $b^*=y_j-\sum_{i=1}^N\alpha_i^*y_i(x_i\cdot x_j)$
常用核函数

多项式核函数 $K(x,z)=(x\cdot z+1)^p$
高斯核函数 $K(x,z)=\exp(-\frac{\|x-z\|^2}{2\sigma^2})$ 函数 $K(x,z)$ 成为核函数的充分必要条件为 $K(x,z)$ 对应的Gram矩阵是半正定矩阵。

将原始问题转化为对偶问题的优点：

对偶问题更容易求解。
可以自然引入核函数，进而推广到非线性分类问题。

2、求解凸二次规划问题

序列最小最优化算法(sequential minimal optimization，SMO)算法用来求解对偶问题： $\begin{aligned} &\min_{\alpha}\quad \frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\alpha_i\alpha_jy_iy_jx_i x_j-\sum_{i=1}^N\alpha_i \\ & \begin{array}{r@{\quad}r@{}l@{\quad}l} s.t.& \sum\limits_{i=1}^N\alpha_iy_i & =0 \\ & \alpha_i & \geq0 & i=1,2,...,N \end{array} \end{aligned}$ SMO算法是一种启发式算法，包括两个部分：
(1) 求解两个变量二次规划的解析方法；
(2) 选择变量的启发式方法。

基本思想是固定 $\alpha_i$ 和 $\alpha_j$ 两个参数以外的所有 $\alpha$ 参数，而 $y_i\alpha_i+y_j\alpha_j$ 为定值， $\alpha_i$ 可以由 $\alpha_j$ 表示。固定其他参数后，就变成了求两个变量的二次规划问题。
第一个变量一般选取违反KKT条件最严重的样本点；
第二个变量一般选取使目标函数增长最快的变量。

3、代码实现

本篇的SVM支持向量机实现原始代码的github地址为：wepe-SVM，本文仅做一些注释。

class SVCSMO():
    """
        Simple implementation of a Support Vector Classification using the
        Sequential Minimal Optimization (SMO) algorithm for training.
    """
    def __init__(self, max_iter=10000, kernel_type='linear', C=1.0, epsilon=0.001, sigma=5.0):
        """
        :param max_iter: maximum iteration
        :param kernel_type: Kernel type to use in training.
                        'linear' use linear kernel function. 线性核函数
                        'quadratic' use quadratic kernel function. 二次核函数
                        'gaussian' use gaussian kernel function. 高斯和函数
        :param C: Value of regularization parameter. C软间隔最大化中的惩罚系数
        :param epsilon: Convergence value. 停机条件中的epsilon
        :param sigma: parameter for gaussian kernel高斯核函数中的sigma
        """
        self.kernels = {
            'linear' : self.kernel_linear,
            'quadratic' : self.kernel_quadratic,
            'gaussian' : self.kernel_gaussian
        }
        self.max_iter = max_iter
        self.kernel_type = kernel_type
        self.C = C
        self.epsilon = epsilon
        self.sigma = sigma

    def fit(self, X, y):
        n, d = X.shape[0], X.shape[1]# n为数据条数，d为数据维数
        alpha = np.zeros((n))# 初始化alpha，全为0，维数等于数据条数n
        kernel = self.kernels[self.kernel_type]# 根据字典获得对应的核技巧
        count = 0# 记录迭代次数，后面如果迭代次数>=max_iter就退出循环
        while True:# 每循环一次count加一，记录遍历次数。
            count += 1
            alpha_prev = np.copy(alpha)# 每次迭代将alpha拷贝一份，这里的copy()为深拷贝
            for j in range(0, n):# SMO算法中的外层循环，本来应该选取违反KKT条件最严重的alpha_j
                # 但是这里对所有alpha_j，j=1,2,...,n-1都进行了计算
                i = self.get_rnd_int(0, n-1, j) # 获得一个[0, n-1]之间的，并不等于j的数
                # alpha_i的选取本来应该选取使目标函数有足够多下降的alpha_i，但是这里是随机选取了一个不等于j的数。
                x_i, x_j, y_i, y_j = X[i,:], X[j,:], y[i], y[j]# 筛出第i条数据和第j条数据，及其对应的y
                k_ij = kernel(x_i, x_i) + kernel(x_j, x_j) - 2 * kernel(x_i, x_j)# 计算eta
                if k_ij == 0:# 如果计算出得eta为0则跳过。eta在对alpha_j的剪辑中需要用到，做分母。
                    continue
                alpha_prime_j, alpha_prime_i = alpha[j], alpha[i]# alpha_j和alpha_i的原始值。
                (L, H) = self.compute_L_H(self.C, alpha_prime_j, alpha_prime_i, y_j, y_i)# 计算alpha_j的边界

                # 计算w和b
                self.w = self.calc_w(alpha, y, X)
                self.b = self.calc_b(X, y, self.w)

                # 计算E_i和E_j
                E_i = self.E(x_i, y_i, self.w, self.b)
                E_j = self.E(x_j, y_j, self.w, self.b)

                # 计算未剪辑的alpha_j
                alpha[j] = alpha_prime_j + float(y_j * (E_i - E_j))/k_ij
                # 根据alpha_j的边界计算剪辑的alpha_j
                alpha[j] = max(alpha[j], L)
                alpha[j] = min(alpha[j], H)
                # 通过alpha_j计算alpha_i
                alpha[i] = alpha_prime_i + y_i*y_j * (alpha_prime_j - alpha[j])
                
            # 如果迭代完一轮之后，alpha的变化的范数小于精度epsilon，则停止迭代。
            diff = np.linalg.norm(alpha - alpha_prev)
            if diff < self.epsilon:
                break
                
            # 如果到达最大迭代次数，则停止迭代。
            if count >= self.max_iter:
                print("Iteration number exceeded the max of %d iterations" % (self.max_iter))
                return

        # 迭代完毕，计算模型最终的w和b
        self.b = self.calc_b(X, y, self.w)
        if self.kernel_type == 'linear':
            self.w = self.calc_w(alpha, y, X)

        # alpha_i>0对应的数据即为支持向量
        alpha_idx = np.where(alpha > 0)[0]
        support_vectors = X[alpha_idx, :]
        return support_vectors, count

    # 预测
    def predict(self, X):
        return self.h(X, self.w, self.b)

    def calc_b(self, X, y, w):
        b_tmp = y - np.dot(w.T, X.T)
        return np.mean(b_tmp)

    def calc_w(self, alpha, y, X):
        return np.dot(alpha * y, X)

    # Predict函数调用h函数
    def h(self, X, w, b):
        return np.sign(np.dot(w.T, X.T) + b).astype(int)

    # 计算预测值与实际值的差
    def E(self, x_k, y_k, w, b):
        return self.h(x_k, w, b) - y_k

    # 计算L、H边界
    def compute_L_H(self, C, alpha_prime_j, alpha_prime_i, y_j, y_i):
        if(y_i != y_j):# alpha_j - alpha_i = k
            return (max(0, alpha_prime_j - alpha_prime_i), min(C, C - alpha_prime_i + alpha_prime_j))
        else:# alpha_j + alpha_i = k
            return (max(0, alpha_prime_i + alpha_prime_j - C), min(C, alpha_prime_i + alpha_prime_j))

    # 产生一个[a,b]之间的数，并且不等于z，用于选取第二个变量
    def get_rnd_int(self, a,b,z):
        i = z
        while i == z:
            i = rnd.randint(a,b)
        return i

    # 线性核函数
    def kernel_linear(self, x1, x2):
        return np.dot(x1, x2.T)
    # 多项式核函数
    def kernel_quadratic(self, x1, x2):
        return (np.dot(x1, x2.T) ** 2)
    # 高斯核函数
    def kernel_gaussian(self, x1, x2, sigma=5.0):
        if self.sigma:
            sigma = self.sigma
        return np.exp(-linalg.norm(x1-x2)**2 / (2 * (sigma ** 2)))

总的来说，SVM就是找到一个超平面使间隔最大化，我们将其转化为对应的对偶问题，方便求解。

结尾

如果您发现我的文章有任何错误，或对我的文章有什么好的建议，请联系我！如果您喜欢我的文章，请点喜欢~*我是蓝白绛，感谢你的阅读！

机器学习｜SVM支持向量机代码实现解读

前言

第七章支持向量机

1、模型

2、求解凸二次规划问题

3、代码实现

结尾

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

机器学习｜SVM支持向量机代码实现解读

前言

第七章 支持向量机

1、模型

2、求解凸二次规划问题

3、代码实现

结尾

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

第七章支持向量机