【机器学习】svm

机器学习算法——SVM

目录

• 机器学习算法——SVM
  • 1. 背景
  • 2. SVM推导
    • 2.1 几何间隔和函数间隔
    • 2.2 SVM原问题
    • 2.3 SVM对偶问题
    • 2.4 SMO算法
      • 2.4.1 更新公式
      • 2.4.2 裁剪
      • 2.4.3 优化变量的选择
      • 2.4.4 偏移$b$和误差$E_i$的更新
  • 3. SVM的python实现
  • 4. 改进

1. 背景

​ 在线性分类任务中，对于同一个数据集，可能有多个分离超平面。例如在下图中，H2和H3都能够将白色点和黑色点分离开来，那么在这些分界面中，是否存在一个最优的分界面？一个直观的想法是，离所有点都比较远的分割面会是一个好的分割面。可以证明，这样的最优分割面是唯一的。因此SVM的目标就变成了寻找最大间隔分离超平面。

2. SVM推导

2.1 几何间隔和函数间隔

​ 对于数据集({x_i, y_i}_{i=1}^N, y_i in{-1, +1})，分类器(f(x)=w^Tx + b)，任意点(x_i)到分割面的几何间隔为：

[d_i = frac{y_i (w^T x_i + b)}{||w||}=y_i left(frac{w^T}{||w||}x + frac{b}{||w||} ight)， d=min_{iin{1, 2, ..., N}} d_i label{geo_dist} ag{1} ]

从公式可知，如果成比例地改变(w)和(b)，并不会影响(d_i)的值。定义函数间隔：

[hat{d}_i = ||w||d_i = y_i left(w^T x_i+b ight), hat{d}=min_{iin{1, 2, ..., N}} hat{d}_i label{func_dist} ag{2} ]

2.2 SVM原问题

​ 最大间隔分离超平面问题可以表述为：

[max_{w,b} d label{p1} ag{3}\ s.t quad d_i = y_i left(frac{w^T}{||w||}x + frac{b}{||w||} ight)geq d, forall i in{1, 2, ...,N} ]

将几何间隔替换为函数间隔，公式( ef{p1})可以写作：

[max_{w,b} frac{hat{d}}{||w||}label{p2} ag{4}\ s.t quad y_i left(w^T x_i+b ight)geq ||w||d = hat{d} , forall i in{1, 2, ...,N} ]

事实上函数间隔的取值并不会影响问题( ef{p2})的最优解（等比例放缩(w)和(b)不会影响不等式约束），取函数间隔(hat{d} = 1)，公式( ef{p2})变为：

[min_{w,b} frac{1}{2}w^Tw label{p3} ag{5}\ s.t quad y_i (w^Tx_i+b)geq 1, forall i in{1, 2, ...,N} ]

​ 然而在某些情况下，最小函数间隔是1并不能成立（非线性可分情况），为了处理这种情形，引入松弛变量(xi_igeq 0)，约束条件变为了函数间隔加上松弛变量要大于1，此时问题( ef{p3})变为：

[min_{w,b} frac{1}{2}w^Tw+Csum_{i=1}^Nxi_i label{p4} ag{6}\ s.t quad y_i (w^Tx_i+b)geq 1 - xi_i\ xi_i geq 0, quad forall i in{1, 2, ...,N} ]

2.3 SVM对偶问题

​ 写出问题( ef{p4})的拉格朗日函数：

[mathcal{L}(w, b, xi, alpha, mu) = frac{1}{2}w^Tw + Csum_{i=1}^Nxi_i + sum_{i=1}^N alpha_i(1-xi_i - y_i(w^T x_i + b)) - sum_{i=1}^N mu_i xi_i, alpha geq 0 label{Lagrangrian} ag{7} ]

原始问题为：( heta_p = min_{w,b,xi} max_{alpha} mathcal{L}(w, b,xi, alpha))，对偶问题为：( heta_d = max_{alpha} min_{w,b,xi}mathcal{L}(w,b,xi,alpha))。求解对偶问题：

[ abla_w mathcal{L}(w, b,xi, alpha) = w - sum_{i=1}^N alpha_i y_i x_i label{grad} ag{8}\ abla_b mathcal{L}(w, b,xi, alpha) = -sum_{i=1}^N alpha_i y_i \ abla_{xi_i} mathcal{L}(w, b,xi, alpha) = C-alpha_i-mu_i , forall iin {1, 2, .., N} ]

根据KKT条件得到：

[w = sum_{i=1}^N alpha_i y_i x_i label{kkt} ag{9}\ sum_{i=1}^N alpha_i y_i = 0 \ C-alpha_i - mu_i = 0 \ alpha_i geq 0 \ mu_i geq 0\ alpha_i (1-xi_i - y_i (w^Tx_i+b)) = 0 \ mu_i xi_i = 0 \ xi_i geq 0 \ y_i (w^Tx_i +b)geq 1 - xi_i ]

将KKT条件代入到对偶问题可得：

[max_alpha frac{1}{2}sum_{i=1}^{N}sum_{j=1}^N alpha_i alpha_j y_i y_j x_i^Tx_j - sum_{i=1}^N alpha_i label{p5} ag{10}\ s.t quad sum_{i=1}^N alpha_i y_i = 0\ 0leq alpha_i leq C , quad forall i in {1, 2, ..., N} ]

同时观察KKT条件可以得到：

• 如果(alpha_i < C)，一定有(xi_i = 0)（(mu_i = C -alpha_i eq 0 Longrightarrow xi_i = 0)），支持向量(x_i)恰好落在边界
• 如果(alpha_i = C, 0<xi_i < 1)，则(x_i)分类正确，在间隔边界与超平面之间
• 如果(alpha_i = C, xi_i = 1)，则(x_i)落在超平面上
• 如果(alpha_i = C, xi_i > 1)，则(x_i)分类错误

2.4 SMO算法

2.4.1 更新公式

​ 问题( ef{p5})仍然是一个二次规划问题，可以用一般的QP算法解决。但是对于SVM模型，有一种更加快速的优化算法。类似于坐标下降，SMO每次固定其他变量，只优化两个变量。

​ 假设在问题( ef{p5})中，选择优化变量(alpha_1, alpha_2)。等式约束可以改写为：

[alpha_1 y_1 + alpha_2 y_2 = -sum_{i=3}^N alpha_i y_i riangleq zeta label{alpha1_alpha2} ag{11}\ alpha_2 = zeta y_2 -alpha_1 y_1 y_2 ]

将待优化的变量分离出来，并将(alpha_2)替换掉：

[egin{aligned} g(alpha_1, alpha_2) =& frac{1}{2}alpha_1^2 K_{11} + frac{1}{2} alpha_2^2 K_{22} + alpha_1 alpha_2 y_1 y_2 K_{12} + alpha_1 y_1sum_{i=3}^N alpha_i y_i K_{1i} + alpha_2 y_2sum_{i=3}^N alpha_i y_i K_{2i} \ &+ frac{1}{2}sum_{i=3}^{N}sum_{j=3}^N alpha_i alpha_j y_i y_j K_{ij} -(alpha_1 + alpha_2) - sum_{i=3}^N alpha_i \ =& frac{1}{2}alpha_1^2 K_{11} + frac{1}{2}left( zeta y_2 - alpha_1 y_1 y_2 ight)^2 K_{22} + alpha_1 left( zeta y_2 - alpha_1 y_1 y_2 ight) y_1 y_2 K_{12}\ &+ alpha_1 y_1 v_1 + left( zeta y_2 - alpha_1 y_1 y_2 ight)y_2 v_2 -alpha_1 - zeta y_2 + alpha_1 y_1 y_2 + mathrm{const} \ =& frac{1}{2} left( K_{11} - 2K_{12} + K_{22} ight)alpha_1^2 + left[zeta y_1(K_{12} - K_{22}) + y_1(v_1 -v_2)+(y_1 y_2 - 1) ight]alpha_1 + mathrm{const} end{aligned}label{alpha1} ag{12} ]

因此有：

[frac{partial g}{partial alpha_1} = (K_{11} - 2K_{12}+K_{22})alpha_1 + zeta y_1(K_{12} - K_{22}) + y_1(v_1 -v_2)+(y_1 y_2 - 1) label{grad_alpha1} ag{13} ]

其中：

[v_1 = sum_{i=3}^N alpha_i y_i K_{1i} = sum_{i=1}^N alpha_i^{old} y_i K_{1i} - alpha_1^{old}y_1 K_{11} - alpha_2^{old}y_2 K_{12} = f(x_1) - alpha_1^{old}y_1 K_{11} - alpha_2^{old}y_2 K_{12} - b label{v} ag{14}\ v_2 = sum_{i=3}^N alpha_2 y_i K_{2i} = sum_{i=1}^N alpha_i^{old} y_i K_{2i} - alpha_1^{old}y_1 K_{21} - alpha_2^{old}y_2 K_{22} = f(x_2) - alpha_1^{old}y_1 K_{21} - alpha_2^{old}y_2 K_{22} - b ]

所以：

[egin{aligned} v_1 - v_2 &= f(x_1) - f(x_2) - alpha_1^{old}y_1(K_{11}-K_{21}) - (zeta y_1 - alpha_1 y_1 y_2)y_2(K_{12}-K_{22}) \ &= f(x_1) - f(x_2) - alpha_1^{old}y_1(K_{11} - 2K_{12} + K_{22}) - zeta y_1 y_2 (K_{12} - K_{22}) end{aligned}label{diff_v} ag{15} ]

带入到公式( ef{grad_alpha1})得到：

[egin{aligned} frac{partial g}{partial alpha_1} &= (K_{11} - 2K_{12}+K_{22})alpha_1^{new} + zeta y_1(K_{12} - K_{22}) \ &+ y_1( f(x_1) - f(x_2) - alpha_1^{old}y_1(K_{11} - 2K_{12} + K_{22}) - zeta y_1 y_2 (K_{12} - K_{22}))+(y_1 y_2 - 1)\ &= (K_{11} - 2K_{12}+K_{22})(alpha_1^{new}-alpha_1^{old}) + y_1 ((f(x_1) - y_1)-(f(x_2) -y_2)) \ &= (K_{11} - 2K_{12}+K_{22})(alpha_1^{new}-alpha_1^{old}) - y_1 (E_1-E_2) end{aligned}label{grad_alpha1_2} ag{16} ]

令该导数为0，得到：

[alpha_1^{new} = alpha_1^{old} + frac{y_1(E_1-E_2)}{K_{11} - 2K_{12}+K_{22}}label{alpha1_new} ag{17} ]

2.4.2 裁剪

​ 现在考虑对偶问题中的框约束。(alpha_1, alpha_2)的等式约束( ef{alpha1_alpha2})可以写作：

[alpha_1^{new} y_1 + alpha_2^{new} y_2 = zeta = alpha_1^{old} y_1 + alpha_2^{old} y_2 label{old_new} ag{18} ]

对(y_1, y_2)分类讨论：

• 当(y_1 eq y_2)时，(alpha_{1}^{new}-alpha_2^{new} = k = alpha_1^{old}-alpha_2^{old})

  • (k>0Longrightarrow alpha_1 in left[k, C ight])
  • (k < 0Longrightarrow alpha_1 in left[0, C+k ight])

  因此有(L = max{0, alpha_1^{old}-alpha_2^{old}}, H = min{C, C + alpha_1^{old}-alpha_2^{old}})

• 当(y_1 = y_2)时，(alpha_{1}^{new}+alpha_2^{new} = k = alpha_1^{old}+alpha_2^{old})，同上一种情况可得(L = max{0, alpha_1^{old}+alpha_2^{old}-C}, H = min{C, alpha_1^{old}+alpha_2^{old}})

对(alpha_1^{new})的裁剪过程如下：

[egin{equation}label{clip} ag{19} alpha_1^{new,cliped}= egin{cases} H&,alpha_1^{new}geq H \ alpha_1^{new}&, L<alpha_1^{new}< H\ L&, alpha_1^{new}leq L end{cases} end{equation} ]

在计算出(alpha_1^{new})之后，代入( ef{old_new})可以得到：

[alpha_2^{new} = alpha_2^{old}+y_1 y_2 (alpha_1^{old}-alpha_2^{old}) label{alpha2_new} ag{20} ]

2.4.3 优化变量的选择

​ (alpha_i, alpha_j)的选择。在选择第一个变量(alpha_i)时，找出违反KKT条件最严重的，这样能加快优化过程。KKT条件具体是：

[egin{aligned} alpha_{i} &=0 Leftrightarrow y_{i} fleft(x_{i} ight) geqslant 1 \ 0<alpha_{i} &<C Leftrightarrow y_{i} fleft(x_{i} ight)=1 \ alpha_{i} &=C Leftrightarrow y_{i} fleft(x_{i} ight) leqslant 1 end{aligned} ]

其中(fleft(x_{i} ight)=sum_{j=1}^{N} alpha_{j} y_{j} Kleft(x_{i}, x_{j} ight)+b)。在检验时，首先在支持向量中寻找，即(0<alpha_{i} <C)，如果支持向量都满足KKT条件，则在全部数据集中寻找。在给定了第一个变量时，第二个变量(alpha_j)的选择要使(alpha_j)有足够大的变化，即使(|E_i -E_j|)最大。

2.4.4 偏移(b)和误差(E_i)的更新

​ 根据KKT条件可知，当向量(x_i)是支持向量时有(y_{i} fleft(x_{i} ight)=1)。在计算完(alpha_1^{new})之后，如果(0<alpha_1^{new}<C)：

[sum_{i=1}^N y_i alpha_i K_{1i} + b= y_1 ]

于是有：

[b = y_1 - sum_{i=1}^N y_i alpha_i^{old} K_{1i} + y_1K_{11}(alpha_1^{old}-alpha_1^{new})+y_2 K_{12}(alpha_2^{old} - alpha_2^{new}) ]

代入(E_1 = y_1 - f(x_1) = y_1 - sum_{i=1}^N alpha_i^{old}y_i K_{1i}-b^{old})得到：

[b^{new} = E_1 + y_1K_{11}(alpha_1^{old}-alpha_1^{new})+y_2 K_{12}(alpha_2^{old} - alpha_2^{new}) + b^{old} ]

同样地当(0<alpha_2^{new}<C)，有：

[b_{2}^{ ext {new }}=E_{2}+y_{1} K_{12}left(alpha_{1}^{ ext {old }}-alpha_{1}^{ ext {new }} ight)-y_{2} K_{22}left(alpha_{2}^{ ext {old }}-alpha_{2}^{ ext {new }} ight)+b^{ ext {old }} ]

当(alpha_1^{new}, alpha_2^{new})是0或者C时，选择(frac{b_1^{new} + b_2^{new}}{2})作为更新值（因为在([b_1^{new}, b_2^{new}])中的数都能满足KKT条件，选中点作为近似）

​ 当更新完(alpha_1, alpha_2, b)之后，误差(E_i)更新为：

[E_{i}^{ ext {new }}=y_i - sum_{S} y_{j} alpha_{j} Kleft(x_{i}, x_{j} ight)+b^{ ext {new }} ]

其中(S)是支持向量的集合。

3. SVM的python实现

SVM有很多实现方式：

• hingle loss；可以证明，线性SVM问题和采用hingle loss的线性分类器是等价的，因此可以使用梯度下降方法进行优化
• SMO；简单的实现是选择一个符合条件的(alpha_i)，然后随机挑选一个(alpha_j)。为了加快训练速度，可以根据某些准则合理地选择
• 二次优化；直接使用QP求解器

本文实现了simple smo，代码如下：

#!/usr/bin/env python 
# -*- coding:utf-8 -*-
# -*- coding: utf-8 -*-

import numpy as np
# 裁剪函数
def clip(x, L, H):
    if(x < L):
        x = L
    elif(x > H):
        x = H
    return x
# 随机选择数字
def rand_select(i, m):
    j = i
    while(j == i): j = int(np.random.uniform(0, m))
    return j

class SVM():
    def __init__(self, kernel="linear", C = 1.0, sigma=1.0, **kwargs):
        # 只支持线性核和高斯核
        if (kernel not in ['linear', 'gaussian']):
            raise ValueError("Now only support linear and gaussian kernel")
        elif kernel == "linear":
            kernel_fn = Kernel.linear()
        elif kernel == "gaussian":
            kernel_fn = Kernel.gaussian(sigma)

        self.kernel_method = kernel
        self.kernel = kernel_fn

        self.sigma = sigma
        self.C = C

        self.w , self.b = None, None
        self.n_sv = -1
        self.sv_x, self.sv_y, self.alphas = None, None, None
        self.eps = 1e-6

    # 预测函数， y = sum_(i in S) alpha_i y_i K(x_i, x) + b
    def predict(self, x):
        wx = self.b
        for i in range(len(self.sv_x)):
            wx += self.alphas[i] * self.sv_y[i] * self.kernel(self.sv_x[i], x)
        return wx

    '''
    svm训练，常见的方法有以下几种:
    1. simple_smo
    2. smo
    3. cvopt
    4. hinge loss
    目前只实现了simple_smo
    '''
    def train(self, X, y, maxIter = 500):
        # 计算核矩阵
        if(self.kernel_method == "linear"):
            K = np.dot(X, X.T)
        elif (self.kernel_method == "gaussian"):
            d2 = np.sum(X ** 2, axis = -1, keepdims= True) - 2 * np.dot(X, X.T) + np.sum(X ** 2, axis = -1,)
            K = np.exp(- 0.5 * d2 / self.sigma)

        # 优化
        alphas, b = self._SMO(K, y, maxIter = maxIter)

        # 计算支持向量
        idx = [i for i in range(len(alphas)) if alphas[i] > self.eps]
        self.n_sv = len(idx)
        self.sv_x = X[idx,:]
        self.alphas = alphas[idx]
        if self.kernel_method == "linear":
            self.w = np.sum((self.alphas * y[idx]).reshape(-1, 1) * X[idx,:], axis = 0)
        self.b = b

    def _SMO(self, K, y, tol = 0.001, maxIter = 40):
        print("begin SMO...")
        m = len(y)
        self.alphas, self.b = np.zeros(m), 0.
        y = y.reshape(1, -1)

        niter = 0
        while (niter < maxIter):
            # select alpha_i, alpha_j
            changed = 0
            for i in range(m):
                # 在外层循环中，选择一个违反了KKT条件的alpha_i
                yi, ai_old = y[0][i], self.alphas[i].copy()
                Ei = np.dot(self.alphas * y , K[i].T) + self.b - yi
                if ((yi * Ei < -tol) and (ai_old < self.C)) or ((yi * Ei > tol) and (ai_old > 0)):
                    # 随机找内层的alpha_j
                    j = rand_select(i, m)
                    yj, aj_old = y[0][j], self.alphas[j].copy()
                    Ej = np.dot(self.alphas * y, K[j].T) + self.b - yj

                    # 更新 alpha_j
                    eta = K[i, i] + K[j, j] - 2 * K[i, j]
                    if(eta <= 0): continue;
                    aj_new = aj_old + yj * (Ei - Ej) / eta

                    # 裁剪 alpha_j
                    if yi != yj:
                        L, H = max(0, aj_old - ai_old), min(self.C, self.C + aj_old - ai_old)
                    else:
                        L, H = max(0, aj_old + ai_old - self.C), min(self.C, aj_old + ai_old)
                    if H - L == 0:
                        continue
                    aj_new = clip(aj_new, L, H)

                    # 更新alpha_i
                    s = yi * yj
                    delta_j = aj_new - aj_old
                    if(abs(delta_j) < self.eps): continue;
                    ai_new = ai_old - s * delta_j
                    delta_i = ai_new - ai_old

                    # 更新b
                    bi = self.b - Ei - yi * delta_i * K[i, i] - yj * delta_j * K[i, j]
                    bj = self.b - Ej - yi * delta_i * K[i, j] - yj * delta_j * K[j, j]

                    if 0 < ai_new < self.C:
                        self.b = bi
                    elif 0 < aj_new < self.C:
                        self.b = bj
                    else:
                        self.b = (bi + bj) / 2
                    self.alphas[i], self.alphas[j] = ai_new, aj_new
                    changed += 1
            if (changed == 0):
                niter +=1;
            else:
                niter = 0

        print("Finish SMO...")
        return self.alphas, self.b

'''
核函数的单独实现，后期可以在里面添加
'''
class Kernel(object):
    # 线性核
    @staticmethod
    def linear():
        return lambda X, y: np.inner(X, y)

    # 高斯核
    @staticmethod
    def gaussian(sigma):
        return lambda X, y: np.exp(-np.sqrt(np.linalg.norm(X - y) ** 2 / (2 * sigma)))

测试代码：

#!/usr/bin/env python 
# -*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

# generate data

anchors = np.array([[-1, -1], [1, 1]]) * 0.5
n1 = np.random.random((30, 2))
n2 = np.random.random((30, 2))
X1 = anchors[0] - n1
X2 = anchors[1] + n2
y1 = np.array([-1]*30)
y2 = np.array([1]*30)

print("trianing...")
X = np.vstack((X1, X2))
y = np.hstack((y1, y2))

idx = np.array(range(len(y)))
np.random.shuffle(idx)
X = X[idx,:]
y = y[idx]

from svm import *
model = SVM(C = 0.6, kernel="linear")
model.train(X, y, maxIter = 40)

print(model.n_sv)
print(model.sv_x)
print(model.sv_y)

xx = np.array([-1, 1])
yy = -model.w[0] / model.w[1] * xx - model.b / model.w[1]
fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(len(y)):
    if y[i] == -1:
        ax.scatter(X[i][0], X[i][1], marker = "o", color = "red")
    else:
        ax.scatter(X[i][0], X[i][1], marker="x", color="green")

for sv_x in model.sv_x:
    cir1 = Circle(xy = (sv_x[0], sv_x[1]), radius=0.1, alpha=0.5)
    ax.add_patch(cir1)

ax.plot(xx, yy, color = "orange")
plt.show()

测试结果如下：

4. 改进

• LRU缓存；在SVM+SMO中的实现中，核矩阵的计算成为了很大的开销。如果预先将核矩阵计算好，空间复杂度为(O(N^2))，如果边用边计算，又会因为重复计算增加开销。考虑到在实际计算中，用到的样本仅为支持向量附近的一些数据点，因此用两个cache保存核和误差，在对参数更新之后更新缓存即可。

• 冷热数据分离；在更新参数时，优先更新支持向量（热数据）对应的(alpha)（即(0<alpha<C)），在没有这样的点时，再全局（冷数据）寻找进行更新。