机器学习05-支持向量机_1

支持向量机

分隔超平面

如下图是一组线性可分的二维点，我们可以通过我们可以使用梯度下降算法找到将测试数据进行区分的直线，由于测试数据区分度很高，因此我们可能找到不止一条合适的分割线即多组权值((w0,w1,w2),(w0’,w1’,w2’)….)，这个时候我们如何判定最合适的分割线？同理扩展到3维，n维都存在这个问题。

因此我们引入分隔超平面的概念，通过确认一个分隔训练集的平面，使得所有训练集合数据尽量的远离该平面，这个平面称之为分隔超平面，对于二维数据超平面为一维的一条直线，对于三维数据超平面为二维平面，对于n维数据超平面为n-1维的面。确定超平面后对于给定数据点如果离超平面的距离越远则最终得到的预测结果越可信。
通过利用支持向量来确定分隔超平面，支持向量指离分隔超平面最近的点，使得支持向量到分隔面的距离最大化，就能找到最优的分隔超平面。

分隔超平面最优化问题

依据前两章的做法定义超平面表达式 $w^T x + b$，其中b相当于w0，要计算点A到超平面的距离，就需要知道平面的法线，超平面的法线值为 $|w^T + b|/ ||w||$。在梯度下降中的logistics回归我们使用了Sigmoid函数定义我们的二分类算法 $h_\theta(x) = g(\theta^Tx)=\frac{1}{1+e^-\theta^Tx}$ 用来表示函数的特征属于y=1的概率：

$P(y=1|x;\theta) = h_\theta(x)$ $P(y=0|x;\theta) = 1-h_\theta(x)$

当P大于0.5时，特征属于y=1的分类，反之属于y=0的分类。在支持向量机中我们将logistics做了一个变形，定义$h_w,_b(x) = g(w^Tx+b)$ 并且定义：

$g(z) = \begin{cases} 1,&\text{z >= 0}\\ -1, & \text{z<0}\end{cases}$

对于超平面$w^Tx + b$当>0时表示数据点在超平面上方分类为1，当<0时表示数据点在超平面下方表示分类为-1，当=0时表示数据点在超平面上这是我们可以默认此时分类为1。

间隔函数

需要求数据点到超平面的距离，我们定义 $y\cdot(w^Tx + b)$为我们的间隔函数，其中$y$表示目标分类 -1 或 1，这样我们的间隔函数值就恒大于等于0。函数的目标是为了找出分类器的参数 $w^Tx$ 和 $b$,我们需要先确定我们的支持向量(也叫最小间隔)，然后使得该最小间隔最大，可以写成：

$$\arg\, \max_{w,b}\lbrace \min_{n}(y\cdot(w^Tx+b)\cdot\frac{1}{||w||})\rbrace$$
公式中$\min_{n}(y\cdot(w^Tx + b )\cdot\frac{1}{||w||})$表示的是超平面的支持向量，$\arg \max_w,_b$使得支持向量距离最大,令所有支持向量$label\cdot(w^Tx +b)$值为1，方程就转换为求 $||w||^{-1}$的最小值。我们又知道只有那些离超平面最近的点的距离才等于1，而其它点都大于1，因此在求值过程中我们需要给函数一个约束条件最终该函数转换为:

$$\begin{cases}\min||w||^2,&\text{式(1)}\\ st. y_i(w^Tx_i + b) -1 >=0 \end{cases}$$

上面的公式比较难于理解，这里有另外一种理解方式：定义两个平面

$$\begin{cases} w^Tx + b = 1,&\text{for y=1}\\ w^Tx + b = -1, &\text{for y=-1}\end{cases}$$ 这两个平面分别表示类别y = 1 和 -1时，离超平面最近的点所在的平面。式子统一转换为 $y\cdot(w^T + b) >=1$如下图中的H1平面 H2平面：

图片来源
我们要实现的目标是使得H1、H2到超平面的距离尽量远即： $Margin=\frac{2}{||w||}$的最大值,于是就能构造出求解方程组：
$$\begin{cases}\min||w||^2,&\text{式(2)}\\ st. y_i(w^Tx_i + b) -1 >=0 \end{cases}$$
这里的式（1）和式（2）为了方便计算可以统一写成 $\min\frac{||w||^2}{2}$。

拉格朗日方法

为了求解我们的约束方程，需要引入拉格朗日函数，拉格朗日函数是用于求解多元函数在收到一个或多个约束条件时的极值问题的方法。使用拉格朗日函数可将一个n个变量和k个约束条件的最优化问题转换为n+k个变量的n+k个变量方程组。比如：要求$f(x,y)$ 在 $g(x,y)=c$时的最大值时，我们可以引入新变量拉格朗日乘数$\alpha$,这是我们只需要下列拉格朗日函数的极值:

$$L(x,y,\alpha) = f(x,y) + \alpha\cdot\left( g(x,y) -c \right)$$
拉格朗日乘法所得的极值会包含原问题的所有极值点，但并不保证每个极值点都是原问题的极值点。

使用拉格朗日方法求解我们的方程得到：

$$L(w,b,\alpha) = \min\frac{1}{2}\cdot||w||^2 - \sum_{i=0}^n \alpha_i\left(y_i(w^Tx_i + b) - 1)\right)$$
要求解 $L(w,b,\alpha)$的最小值并且使之等价于 $\min||w||^2$， 又因为 $y_i(w^Tx_i + b) >=1$，所有我们需要先求解 $\max_{\alpha>=0}L(w,b,\alpha)$,因此原始方程解变成求解：
$$\min_{w,b}\,\max_{\alpha_i>=0} L(w,b,\alpha) = \max_{\alpha_i>=0}\min_{w,b}L(w,b,\alpha)$$

我们先求$\min_{w,b}L(w,b,\alpha)$,先对w，b分别求偏导，令其等于0：

$$\frac{d_L(w,b,\alpha)}{d_w} = ||w|| - \sum_{i=0}^n\alpha_iy_ix_i$$
$$\frac{d_L(w,b,\alpha)}{d_b} = - \sum_{i=0}^n\alpha_iy_i$$ 得到： $$||w|| = \sum_{i=0}^n\alpha_iy_ix_i\, ,\,\,\,\sum_{i=0}^n\alpha_iy_i = 0$$ 带入 $L$ 得： $$\begin{cases}L = \sum_{i}^n\alpha_i - \frac{1}{2}\sum_{i,\,j}^ny_iy_j\alpha_i\alpha_j<x_i\cdot x_j>,&\text{式(2)}\\ st. \alpha_i >=0,\\ st.\sum_{i=0}^n\alpha_iy_i = 0 \end{cases}$$

推导公式可以参见 参考链接

关于离群问题

很多情况下给定的训练数据不是完全线性可分的，可能存在离群点，即y=-1的点对应的特征值很接近超平面甚至超出了超平面进入到y=1的类别中，这些点会影响超平面的移动甚至会导致超平面最终无解。我们引入了一个松弛变量$\xi$来表示点的偏移程度,目的是我们允许H1 和 H2之间存在数据点。
我们的超平面转变为：$subject \,\,to ：y_i(w\cdot x_i - b - 1 + \xi) >= 0$ 并且 $\xi >=0$
最优值：$\min_{w,b,\xi}\frac{1}{2}w^Tw + C\sum_{i}\xi_i$
这里的C表示离群点对结果影响的权重值（当C=$\infty$时退化到原来的理想线性可分场景）,对上面的公式采用之前的拉格朗日方法推导得到：

$$\begin{cases}L(w,b,\xi,\alpha,u) = \frac{1}{2}w^Tw + C\sum_{i=1}^n\xi_i - \sum_{i=1}^n\alpha_i[y_i(w^Tx_i - b) + \xi_i -1] - \sum_{i=1}^n u_i\xi_i\\ st.\alpha_i(y_i(w^Tx_i -b) + \xi_i -1) =0 \\ st. u_i\xi_i = 0\end{cases}$$

最终得到：

$$\begin{cases}L = \sum_{i}^n\alpha_i - \frac{1}{2}\sum_{i,\,j}^ny_iy_j\alpha_i\alpha_j<x_i\cdot x_j>\\ st. 0<= \alpha_i <=C\\ st.\sum_{i=0}^n\alpha_iy_i = 0 \end{cases}$$

这里可以看到增加松弛变量和惩罚系数后我们的解表达式不变，只是$\alpha_i$ 的取值范围有原来的 $0<= \alpha_i < \infty$变成了 $0<= \alpha_i <=C$

==这里的推导详见底部的引用链接 smo.pdf 文档==

优化过程

这里的目的是在代码实现前，先看一下算法的优化问题和回归问题。
超平面$w^Tx + b =$,算法的最终目的是获取最优的$w^T 、b$,通过前面的推导公式 $w=\sum_{i=0}^n\alpha_iy_ix_i$ 可以转而得到 $w$因此我们的目标转而获取 $\alpha_i 、 b$, 在优化过程中我们每次选两个$\alpha_1 、\alpha_2$求最优化，然后再根据$\alpha_1、、alpha_2$优化我们的 $b$, 这里只写一下回归公式结果，具体推导过程都在==参考的smo.pdf里面写的很详细==。

$$\begin{cases} \eta = - ||x_2 - x_1||^2 \\ w = \sum_{i=1}^N\alpha_iy_ix_i \\ E_2^{old} = x_2^Tw_{old} + b_{old} - y_2 \\ \alpha_2^{new} = \alpha_2^{old} + \frac{y_2(E_2^{old}-E_1^{old})}{\eta} \\ \alpha_1^{new} = \alpha_1^{old} + y_1y_2\Delta\alpha_2\\ \Delta b = E(x,y)^{old} + \Delta \alpha_1y_1x_1^Tx + \Delta \alpha_2y_2x_2^Tx \end{cases}$$

算法实现

支持向量机的实现都是参见<<机器学习实战>>里面的代码。测试数据详见参考里的 svnTestSet.txt
简单版本的实现代码：

'''
Created on Nov 4, 2010
Chapter 5 source file for Machine Learing in Action
@author: Peter
'''
from numpy import *
from time import sleep

def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
    b = 0; m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            #下面两行源自公式 aixiyi x + b - y >= 0 使得点能够在H1 和 H2之外
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
            #判断是否在误差区间,如果超过误差区间则需要作调整
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                ## Two Lagrange Multiplier 问题 smo.pdf 里面也有详细解释
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: print "L==H"; continue
                # (x1 - x2)^2
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta == 0: print "eta>=0"; continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
                                                                         #the update is in the oppostie direction
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
                print "iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
        print "iteration number: %d" % iter
    return b,alphas

def calcWs(alphas,dataArr,classLabels):
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

dataMat, labelMat = loadDataSet("svnTestSet.txt")
b,result = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)

result:
>>> b
matrix([[-3.73850318]])
>>> alphas[alphas>0]
matrix([[ 0.1516457 ,  0.14693547,  0.06334441,  0.00565379,  0.35627179]])
#这里可以看到，一共找到五个支持向量

核函数

考虑如下分布图，在这种case下我们无法通过上面的松弛变量和惩罚参数来进行优化，这个时候我们就需要将这种无法区分的训练数据通过我们核函数(kernel函数)将数据转换为易于分类器理解的形式。通过核函数我们将我们特征数据从一个特征空间转换到另外一个特征空间，在新的空间下我们能很方便的利用smo算法对数据进行处理（另外一种理解是从低维度转换到高纬度）。

当一个不可分的低维数据通过Kernel转到到高维之后的视角：

==径向基核函数==是目前使用最广泛的Kernel函数，当然还有其他类型的Kernel函数（也可以自己指定自己的核函数）例如：多项式核函数，线性核函数

径向基核函数

径向基函数能够基于变量距离运算输出一个标量，这个距离可以是<0,0>向量与其他向量距离，径向基函数的高斯版本：

$$k(x,y) = exp(\frac{-||x-y||^2}{2\sigma})$$
这个公式可以理解为将原始空间映射为无穷维空间。

核函数代码：

def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    elif kTup[0]=='rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

参考：
http://blog.csdn.net/v_july_v/article/details/7624837
http://blog.pluskid.org/?p=685
http://blog.sina.com.cn/s/blog_4298002e010144k8.html
https://zh.wikipedia.org/wiki/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E4%B9%98%E6%95%B0
ftp://www.ai.mit.edu/pub/users/tlp/projects/svm/svm-smo/smo.pdf 或 http://oe7d0gss7.bkt.clouddn.com/smo.pdf
http://cs229.stanford.edu/notes/cs229-notes3.pdf
http://oe7d0gss7.bkt.clouddn.com/svnTestSet.txt