SVM算法

本文主要介绍支持向量机理论推导及其工程应用。

1 基本介绍

支持向量机算法是一个有效的分类算法，可用于分类、回归等任务，在传统的机器学习任务中，通过人工构造、选择特征，然后使用支持向量机作为训练器，可以得到一个效果很好的base-line训练器。

支持向量机具有如下的优缺点，

优点：

1. 高维空间有效；
2. 维度大于样本数量的情况下，依然有效；
3. 预测时使用训练样本的子集（也即支持向量），节省内存；
4. 可以使用不同的核函数用于决策；

缺点：

1. 如果特征的数目远远大于样本的数目，性能将会降低；
2. 不能直接提供概率估计，需要通过5-fold 交叉验证来获得；

2 理论推导

2.1 间隔与支持向量

划分超平面，(mathbf{w}^{T} * mathbf{x} + b = 0 (1))；

样本空间中任意一个样本 (mathbf{x}) 到超平面的距离为 (r = frac{||mathbf{w}^{T} * mathbf{x} + b||}{||mathbf{w}||} (2))；

假设超平面能将样本正确分类，即对样本 ((mathbf{x}_{i},y_{i}))，

若 (y_{i} = +1) ，则 (mathbf{w}^{T} * mathbf{x}_{i} + b > 0)；
若 (y_{i} = -1) ，则 (mathbf{w}^{T} * mathbf{x}_{i} + b < 0)；

令,

(mathbf{w}^{T} * mathbf{x}_{i} + b ge +1, y_{i} = +1 (3))

(mathbf{w}^{T} * mathbf{x}_{i} + b le -1, y_{i} = -1)

距离超平面最近的几个训练样本点，使上式的等号成立，则它们称为支持向量；

两个异类支持向量到超平面的距离之和为

(r = frac{2}{||mathbf{w}||} (4))

欲找到最大间隔的划分超平面，也就是找到 (mathbf{w}) 和b，使得r最大，即

(max_{mathbf{w},b} frac{2}{||mathbf{w}||} (5))

(s.t. y_{i}(mathbf{w}^{T} * mathbf{x}_{i} + b) ge 1, i = 1,2...,m)

等价于，

(min_{mathbf{w},b}frac{1}{2}||mathbf{w}||^{2} (6))

(s.t. y_{i}(mathbf{w}^{T} * mathbf{x}_{i} + b) ge 1, i = 1,2...,m)

希望求解式（6），来得到最大间隔划分超平面所对应的模型，

(f(mathbf{x})=mathbf{w}^{T} * mathbf{x} + b (7))

2.2 对偶问题

对上述公式使用拉格朗日乘子法，可得其对偶问题，具体就是对上述公式的每个约束添加拉格朗日乘子 (alpha_{i} ge 0) ,则拉格朗日函数可写为，

(L(mathbf{w},b,mathbf{alpha})=frac{1}{2}||mathbf{w}||^{2} + sum_{i=1}^{m}alpha_{i} * (1 - y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b)) (8))

其中， (mathbf{alpha} = (alpha_{1};alpha_{2};...;alpha_{m}))，令(L(mathbf{w},b,mathbf{alpha})) 分别对 (mathbf{w}) 和b求偏导，可得，

(mathbf{w} = sum_{i=1}^{m}alpha_{i} * y_{i} * mathbf{x}_{i} (9))

(0 = sum_{i=1}^{m}alpha_{i} * y_{i} (10))

将式（9）带入（8）中，可将(L(mathbf{w},b,mathbf{alpha})) 中的 (mathbf{w}) 和b消去，再考虑式（10）的约束，可得式（6）的对偶问题，

(max_{mathbf{alpha}}sum_{i=1}^{m}alpha_{i} - frac{1}{2}sum_{i=1}^{m}sum_{j=1}^{m}alpha_{i} * alpha_{j} * y_{i} * y_{j} * mathbf{x}_{i}^{T} * mathbf{x}_{j} (11))

(s.t. sum_{i=1}^{m}alpha_{i} * y_{i} = 0)

(alpha_{i} ge 0 i = 1,2...,m)

解出 (mathbf{alpha}) 后，求出 (mathbf{w}) 和b，即可得到模型，

(f(mathbf{x})=mathbf{w}^{T} * mathbf{x} + b = sum_{i=1}^{m}alpha_{i} * y_{i} * mathbf{x}_{i}^{T} * mathbf{x}_{j} + b (12))

式（11）中解出的 (alpha_{i}) 是式（8）中的拉格朗日乘子，恰好对应着样本 ((mathbf{x}_{i},y_{i}))，注意到式（6）中有不等式约束，因此上述过程满足KKT条件，即，

(alpha_{i} ge 0 (13))

(y_{i} * f(mathbf{x}_{i}) - 1 ge 0)

(alpha_{i} * (y_{i} * f(mathbf{x}_{i}) - 1) = 0)

对于训练样本 ((mathbf{x}_{i},y_{i}))，总有 (alpha_{i} = 0) 或者 (y_{i} * f(mathbf{x}_{i}) = 1)；

若 (alpha_{i} = 0) ，则该样本不会在式（12）中的求和中出现，也就不会对 (f(mathbf{x})) 有任何影响；

若 (alpha_{i} ge 0) ,则必有 (y_{i} * f(mathbf{x}_{i}) = 1) ，所对应样本点位于最大分割边界上，是一个支持向量；

2.3 核函数

如果训练样本不能线性可分，可将样本从原始空间映射到更高维的特征空间，使得样本在这个特征空间内线性可分；

令 (phi(mathbf{x})) 表示为 (mathbf{x}) 映射后的特征向量，于是，在特征空间中划分超平面所对应的模型可表示为 (f(mathbf{x})=mathbf{w}^{T}phi(mathbf{x}) + b (14))

其中， (mathbf{w}) 和b是模型参数，类似于式（6）有，

(min_{mathbf{w},b}frac{1}{2}||mathbf{w}||^{2} (15))

(s.t. y_{i}(mathbf{w}^{T} * phi(mathbf{x}_{i}) + b) ge 1, i = 1,2...,m)

对偶问题为，

(max_{mathbf{alpha}}sum_{i=1}^{m}alpha_{i} - frac{1}{2}sum_{i=1}^{m}sum_{j=1}^{m}alpha_{i} * alpha_{j} * y_{i} * y_{j} * phi(mathbf{x}_{i})^{T} * phi(mathbf{x}_{j}) (16))

(s.t. sum_{i=1}^{m}alpha_{i} * y_{i} = 0)

(alpha_{i} ge 0 i = 1,2...,m)

由于样本 (mathbf{x}_{i},mathbf{x}_{j}) 映射到特征空间之后，特征空间的维数可能很高，甚至是无穷维，因此直接计算 (phi(mathbf{x}_{i})^{T} * phi(mathbf{x}_{j})) 通常很困难，可以设想这样一个函数，

(k(mathbf{x}_{i},mathbf{x}_{j}) = <phi(mathbf{x}_{i}),phi(mathbf{x}_{j})> = phi(mathbf{x}_{i})^{T} * phi(mathbf{x}_{j}) (17))

则式（16）可重写为，

(max_{mathbf{alpha}}sum_{i=1}^{m}alpha_{i} - frac{1}{2}sum_{i=1}^{m}sum_{j=1}^{m}alpha_{i} * alpha_{j} * y_{i} * y_{j} * k(mathbf{x}_{i},mathbf{x}_{j}) (18))

(s.t. sum_{i=1}^{m}alpha_{i} * y_{i} = 0)

(alpha_{i} ge 0 i = 1,2...,m)

求解即可得到模型，

(f(mathbf{x})=mathbf{w}^{T} * phi(mathbf{x}) + b = sum_{i=1}^{m}alpha_{i} * y_{i} * phi(mathbf{x}_{i})^{T} * phi(mathbf{x}_{j}) + b = sum_{i=1}^{m}alpha_{i} * y_{i} * k(mathbf{x}_{i},mathbf{x}_{j}) + b (19))

其中，(k(.,.)) 就是核函数；

常见的核函数有，

线性核，(k(mathbf{x}_{i},mathbf{x}_{j}) = mathbf{x}_{i}^{T} * mathbf{x}_{j})；

多项式核，(k(mathbf{x}_{i},mathbf{x}_{j}) = (mathbf{x}_{i} * mathbf{x}_{j})^d,d ge 1) ,为多项式的次数；

高斯核，(k(mathbf{x}_{i},mathbf{x}_{j}) = exp(-frac{||mathbf{x}_{i} - mathbf{x}_{j}||^{2}}{2 * sigma^{2}}),sigma > 0)；

拉普拉斯核，(k(mathbf{x}_{i},mathbf{x}_{j}) = exp(-frac{||mathbf{x}_{i} - mathbf{x}_{j}||}{sigma}).sigma > 0)；

Sigmoid核，(k(mathbf{x}_{i},mathbf{x}_{j}) = tanh(eta * mathbf{x}_{i}^{T} * mathbf{x}_{j} + heta),eta > 0, heta < 0)；

2.4 软间隔

前面介绍的支持向量机形式要求所有样本满足约束条件（3）,即所有样本都必须划分正确，也即“硬划分”，而软间隔则是允许某些样本不满足约束

(y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b) ge 1 (20))

在最大化间隔的同时，不满足约束的样本应该尽可能少，于是，优化目标可写为，

(min_{mathbf{w},b} frac{1}{2}||mathbf{w}||^{2} + C * sum_{i=1}^{m}ell_{0/1}(y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b) - 1) (21))

其中C > 0, (ell_{0/1}) 是“0-1损失函数”,

(ell = 1, if z < 0 (22))

(ell = 0,otherwise)

如果C为无穷大，则式（21）迫使所有样本均满足约束（20），于是式（21）等价与式（6）；

如果C取有限值，式（21）允许一些样本不满足约束（20）；

由于 (ell_{0/1}) 非凸，非连续，因此可以选用其它函数代替它，如hinge-loss，exponential-loss，logistic-loss等损失函数。

若采用hinge-loss函数，则式（21）变为，

(min_{mathbf{w},b} frac{1}{2}||mathbf{w}||^{2} + C * sum_{i=1}^{m}max(0, 1 - y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b)) (23))

引入松弛变量 (xi_{i} ge 0),则式（23）重写为，

(min_{mathbf{w},b,xi_{i}} frac{1}{2}||mathbf{w}||^{2} + C * sum_{i=1}^{m}xi_{i} (24))

(s.t. y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b) ge 1 - xi_{i})

(xi_{i} ge 0, i = 1,2,...,m)

与式（8）类似，通过拉格朗日乘子法可得式（24）的拉格朗日函数。

(L(mathbf{w},b,xi,mu) = frac{1}{2}||mathbf{w}||^{2} + C * sum_{i=1}^{m}xi_{i} + sum_{i=1}^{m}alpha_{i} * (1 - xi_{i} - y_{i} * (mathbf{w}^{T} * mathbf{x}_{i} + b)) - sum_{i=1}^{m}mu_{i} * xi_{i} (25))

其中， (alpha_{i} ge 0, mu_{i} ge 0)是拉格朗日乘子；

令(L(mathbf{w},b,xi,mu)) 分别对 (mathbf{w},b,xi_{i}) 求偏导，可得，

(mathbf{w} = sum_{i=1}^{m}alpha_{i} * y_{i} * mathbf{x}_{i} (26))

(0 = sum_{i=1}^{m}alpha_{i} * y_{i} (27))

(C = alpha_{i} + mu_{i} (28))

将式（26）-（28）代入式（25）中，可得式（24）的对偶问题，

(max_{mathbf{alpha}} = sum_{i=1}^{m}alpha_{i} - frac{1}{2}sum_{i=1}^{m}sum_{j=1}^{m}alpha_{i} * alpha_{j} * y_{i} * y_{j} * mathbf{x}_{i}^{T} * mathbf{x}_{j} (29))

(s.t. sum_{j=1}^{m}alpha_{i} * y_{i} * mathbf{x}_{i} = 0)

(0 le alpha_{i} le C,i=1,2,...,m)

对比式（29）与式（11），唯一的差别就是在于对偶变量的约束不同，前者是 (0 le alpha_{i} le C)，后者是 (0 le alpha_{i})，解法与之前一样，引入核函数后能得到式（19）同样的支持向量展开式；

类似于式（13），对软间隔支持向量机，KKT条件要求。

(C ge alpha_{i} ge 0,mu_{i} ge 0 (30))

(y_{i} * f(mathbf{x}_{i}) - 1 + xi_{i} ge 0)

(alpha_{i} * (y_{i} * f(mathbf{x}_{i}) - 1 + xi_{i}) = 0)

(xi_{i} ge 0, mu_{i} * xi_{i} = 0)

对任意训练样本 ((mathbf{x}_{i},y_{i}))，总有 (alpha_{i} = 0) 或者 (y_{i} * f(mathbf{x}_{i}) = 1 - xi_{i})；

（1） 若 (alpha_{i} = 0) ，则该样本不会对(f(mathbf{x})) 有任何影响；

（2） 若 (alpha_{i} > 0) ，则必有 (y_{i} * f(mathbf{x}_{i}) = 1 - xi_{i})，则该样本为支持向量；

（2.1） 若 (alpha_{i} < C) ，则 (mu_{i} > 0, xi_{i} = 0)，必有 (y_{i} * f(mathbf{x}_{i}) = 1)，则该样本位于最大间隔上；

（2.2） 若 (alpha_{i} = C) ，则 (mu_{i} = 0, xi_{i} > 0)，

（2.2.1） 若 (xi_{i} le 1)，则样本落在最大分割内部；

（2.2.2） 若 (xi_{i} ge 1)，则样本被错误分类；

2.5 多分类

SVM算法最初是为二分类问题设计，当处理多分类问题时，需要构造合适的多类分类器。

1.one-versus-rest（一对多法）：

训练时，依次把某个类别的样本归为一类，其他剩余的样本归为另一类，有k个类别的样本就构造出k个SVM；预测时，将未知样本分类为具有最大分类函数值的那类。

缺陷：会存在数据倾斜；分类结果出现重叠（属于多个分类器）或者不可分类（不属于任何一个分类器）。

2.one-versus-one（一对一法）：

训练时，选择一个类的样本作为正样本，负样本则只选择一个类，又k个类别的样本，就构造出 (frac{k(k-1)}{2}) 个分类器，虽然分类器的数组增加了，但是训练阶段所用的总时间却比"one-versus-rest"方法少很多。预测时，每个分类器都会预测出一个结果，然后统计最后的预测结果。尽管这个方法也有分类重叠现象，但是不会有不可分类的现象，因为不可能所有类别的票数都是0。

3.DAG SVM：

类似于"one-versus-one"方法，只是在对一个样本进行分类之前，先按照下图的结构组织分类器（这是一个有向无环图，因此被称作DAG SVM），

在预测时，可以先问分类器"1 vs 5"，如果回答是5，就往左走；再问分类器"2 vs 5"，如果还回答5，就继续往左走，一直问下去，就可以得到分类结果，如果有k个类别，那么只调用k-1个，分类速度快，且没有分类重叠和不可分类现象。

缺陷：如果一开始分类器回答错误，那么后面的分类器是无法纠正的。

3 工程应用

使用Python机器学习开源库scikit-learn中提供的SVM算法来进行实验。scikit-learn中提供的SVM模型既可以支持稠密数据（numpy.ndarray），也支持稀疏数据（scipy.sparse）。如果要使用SVM模型来对稀疏数据进行预测，它必须符合这些数据类型，例如，为了获得最优的性能，在使用C-ordered numpy.ndarray（稠密）或者scipy.sparse.csr_matrix（稀疏）这些数组时，指定数据类型dtype=float64。

example 1，画出不同SVM分类器在iris数据集上的分界线；

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm,datasets

iris = datasets.load_iris()
X = iris.data[:,:2]
y = iris.target

h = 0.02
C = 1.0
svc = svm.SVC(kernel = 'linear',C=C).fit(X,y)
rbf_svc = svm.SVC(kernel = 'rbf',gamma = 0.7, C=C).fit(X,y)
poly_svc = svm.SVC(kernel = 'poly',degree = 3,C=C).fit(X,y)
lin_svc = svm.LinearSVC(C=C).fit(X,y)

x_min,x_max = X[:,0].min() - 1,X[:,0].max() + 1
y_min,y_max = X[:,1].min() - 1,X[:,1].max() + 1

xx,yy = np.meshgrid(np.arange(x_min,x_max,h),
                    np.arange(y_min,y_max,h))

titles = ['SVC with linear kernel',
          'LinearSVC (linear kernel)',
          'SVC with RBF kernel',
          'SVC with polynomial (degree 3) kernel']

for i ,clf in enumerate((svc,lin_svc,rbf_svc,poly_svc)):
    plt.subplot(2,2,i+1)
    plt.subplots_adjust(wspace = 0.4,hspace = 0.4)
    Z = clf.predict(np.c_[xx.ravel(),yy.ravel()])
    Z = Z.reshape(xx.shape)
    plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.8)
    plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)
    plt.xlabel('Sepal length')
    plt.ylabel('Sepal width')
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.xticks(())
    plt.yticks(())
    plt.title(titles[i])

plt.show()

各个分类器的分界面如下所示，

example 2，画出最大分割超平面；

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm

# we create 40 separable points
np.random.seed(0)
X = np.r_[np.random.randn(20, 2) - [2, 2], np.random.randn(20, 2) + [2, 2]]
Y = [0] * 20 + [1] * 20

# fit the model
clf = svm.SVC(kernel='linear')
clf.fit(X, Y)

# get the separating hyperplane
w = clf.coef_[0]
a = -w[0] / w[1]
xx = np.linspace(-5, 5)
yy = a * xx - (clf.intercept_[0]) / w[1]

# plot the parallels to the separating hyperplane that pass through the
# support vectors
b = clf.support_vectors_[0]
yy_down = a * xx + (b[1] - a * b[0])
b = clf.support_vectors_[-1]
yy_up = a * xx + (b[1] - a * b[0])

# plot the line, the points, and the nearest vectors to the plane
plt.plot(xx, yy, 'k-')
plt.plot(xx, yy_down, 'k--')
plt.plot(xx, yy_up, 'k-.')

plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1],
            s=80, facecolors='none')
plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired)

plt.axis('tight')
plt.show()

最大分割超平面如下所示，

example 3，针对异或类型数据，研究SVM（高斯核）中参数对分类的影响。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm

xx, yy = np.meshgrid(np.linspace(-3, 3, 500),
                     np.linspace(-3, 3, 500))
np.random.seed(0)
X = np.random.randn(300, 2)
Y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)

CArray = [0.1,1.0,10.0,100.0]
gammaArray = [0.1,1.0,10.0,100.0]

clfList = []
for C in CArray:
    clf = svm.SVC(kernel = 'rbf',gamma = gammaArray[1], C=C).fit(X,Y)
    clfList.append(clf)

for gamma in gammaArray:
    clf = svm.SVC(kernel = 'rbf',gamma = gamma, C=CArray[1]).fit(X,Y)
    clfList.append(clf)

titles = ['C = 0.1,gamma = 1.0',
          'C = 1.0,gamma = 1.0',
          'C = 10.0,gamma = 1.0',
          'C = 100.0,gamma = 1.0',
          'C = 1.0,gamma = 0.1',
          'C = 1.0,gamma = 1.0',
          'C = 1.0,gamma = 10.0',
          'C = 1.0,gamma = 100.0',
          ]

# fit the model
for i,clf in enumerate(clfList):
    print "i = ",i
    plt.subplot(2, 4, i + 1)
    plt.subplots_adjust(wspace=0.4, hspace=0.4)
    # plot the decision function for each datapoint on the grid
    Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)


    plt.imshow(Z, interpolation='nearest',
           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',
           origin='lower', cmap=plt.cm.PuOr_r)
    contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,
                       linetypes='--')
    plt.scatter(X[:, 0], X[:, 1], s=30, c=Y, cmap=plt.cm.Paired)
    plt.xticks(())
    plt.yticks(())
    plt.axis([-3, 3, -3, 3])
    plt.title(titles[i])
plt.show()

分类对比图，如下所示，

上面四个图是固定 (gamma) ，选取不同的C，可以发现C越大时，允许被错误分类的样本数量就越少；C越小时，允许被错误分类的样本数量就可能就会越多。

下面四个图是固定C，选取不同的 (gamma) ， (gamma = frac{1}{2*sigma^{2}}) 会影响高斯核函数的分布， (gamma) 越大，高斯核函数分布就会越陡峭； (gamma) 越小，高斯很函数分布就会越平缓； 因此当 (gamma) 越大时，截取等高面时，样本就只会在自己周围形成分布； (gamma) 越小时，截取等高面时，样本就可以在自己周围较大的范围内形成分布，因此同类样本就有可能连接在一起。

4 Reference

机器学习（周志华）

支持向量机（五）SMO算法

SVM-支持向量机详解（四）--多类分类器

SVM的两个参数 C 和 gamma

1.4. Support Vector Machines