支持向量机：SVM

  SVM 是一种监督式的机器学习算法，可用于分类或回归问题。它使用一种称为核函数的技术来变换数据，然后基于这种变换，算法找到预测可能的两种分类之间的最佳边界。通俗来讲，它是一种二类分类模型，其基本模型定义为特征空间上的间隔最大的线性分类器，即支持向量机的学习策略便是间隔最大化，最终可转化为一个凸二次规划问题的求解。

1. 最大间隔分隔

  上图中绿色和蓝色分别表示不同的两个类别，可以看出数据是线性可分，第二张图中间的直线就是一个分类函数，它可以将两类样本完全分开。这种可以将数据集分隔开来的直线称为分隔超平面。在上面给出的数据点都在二维平面上，所以此时分隔超平面就只是一条直线。但是，如果所给的数据集是三维的，那么此时用来分隔数据的就是一个平面。显而易见，更高维的情况可以依此类推。
  支持向量就是离分隔超平面最近的那些点。接下来要试着最大化支持向量到分隔面的距离，要找到这些支持向量。
  我们定义分隔超平面的表达式为$\text{y}(x)=w^{T}x+b$，分类的结果为${\text{y}}_i=\pm 1$(正例 $label$ 为1，负例 $label$ 为 -1)。如果数据点处于正方向（即正例）并且离分隔超平面很远的位置时， $w^{T}x+b$会是一个很大的正数，同时$label\ast (w^{T}x+b)$也会是一个很大的正数。而如果数据点处于负方向（负例）并且离分隔超平面很远的位置时，此时由于类别标签为-1，则$label\ast (w^{T}x+b)$仍然是一个很大的正数。
  现在的目标就是找出分类器定义中的$w$和$b$。为此，我们必须找到具有最小间隔的数据点,一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。这就可以写作： arg ⁡ max ⁡ w , b { min ⁡ n ( l a b e l ∗ ( w T x + b ) ) ∗ 1 ∣ ∣ w ∣ ∣ } \arg \mathop {\max }\limits_{w,b} \{ \mathop {\min }\limits_n (label*({w^T}x + b))*\frac{1}{{||w||}}\} argw,bmax​{nmin​(label∗(wTx+b))∗∣∣w∣∣1​}
  注：
    1. $label\ast (w^{T}x+b)$被称为点到超平面的函数间隔
    2. $label\ast (w^{T}x+b))\ast \frac{1}{||w||}$为点到超平面的几何间隔
    3. $\underset{n}{\min }(label\ast (w^{T}x+b))\ast \frac{1}{||w||}$为支持向量，即上图中红色的点
  直接求解上述问题相当困难，但是如果我们令所有支持向量的$label\ast (w^{T}x+b)$都为1，那么就可以通过求$\frac{1}{||w||}$的最大值来得到最终解，但是并不是所有点的$label\ast (w^{T}x+b)$都为1，只有离超平面最近的点才为1，越远的点值越大，所以这里我们加个约束条件：如果$label\ast (w^{T}x+b)>=1$则$label\ast (w^{T}x+b)=1$，那么到超平面最近的距离为$\frac{1}{||w||}$，即最大的间隔分隔为$\frac{1}{||w||}$。

2. 最优化求解
  上面我们已经分析出最大的间隔为$\frac{1}{||w||}$，$||w||$的意思是$w$的二范数，可以写成$\sqrt{w.w}$，所以$$\max (\frac{1}{||w||})<=>\min (\frac{1}{2}||w|{|}^2)$$

  到这里，最优化问题可以使用拉格朗日乘子法去解，使用了KKT条件的理论，这里直接作出这个式子的拉格朗日目标函数：$$L(w,b,a)=\frac{1}{2}||w|{|}^2-\sum _{i=1}^n{a}_i(y_i(w{x}_i+b)-1)$$  首先让$L$关于$w$，$b$最小化，分别令$L$关于$w$，$b$的偏导数为$0$，得到: $$\frac{\partial L}{\partial w}=0=>w=\sum _{i=1}^n{a}_i{x}_i{y}_i$$$$\frac{\partial L}{\partial b}=0=>\sum _{i=1}^n{a}_i{y}_i=0$$  带入上面的公式得到：
$$L(w,b,a)=\frac{1}{2}||w||2-\sum _{i=1}^n{a}_i{x}_i{y}_{i}w+\sum _i^n{a}_i=\sum _i^n{a}_i-\frac{1}{2}||w|{|}^2=\sum _i^n{a}_i-\frac{1}{2}\sum _{i,j=1}^n{a}_i{a}_j{y}_i{y}_j{x}_i{x}_j$$其中约束条件为：
          $a_i>=0$和$\sum _{i=1}^n{a}_i{y}_i=0$
  这个就是我们需要最终优化的公式。

3. 松弛变量
  我们知道数据都不那么干净，所以通过引入松弛变量来允许数据点可以处于分隔面错误的一侧。这样我们的优化目标就能保持仍然不变，但是新的约束条件则变为：
          $C>a_i>=0$和$\sum _{i=1}^n{a}_i{y}_i=0$
  C用于控制“最大化间隔和保证大部分点的函数间隔小于1.0这两个目标的权重，C值越大，表示离群点影响越大，就越容易过度拟合，反之有可能欠拟合。例如：正例有10000个，而负例只给了100个，C越大表示100个负样本的影响越大，就会出现过度拟合，所以C决定了负样本对模型拟合程度的影响。

4. 简单$SMO$ 算法
  在优化算法的实现代码中，常数C是一个参数，因此我们就可以通过调节该参数得到不同的结果。一旦求出了所有的$alpha$(即$a_i$)，那么分隔超平面就可以通过这些$alpha$来表达。
  $SMO$算法是将大优化问题分解为多个小优化问题来求解的，目标是求出一系列$alpha$和$b$，一旦求出了这些$alpha$，就很容易计算出权重向量$w$并得到分隔超平面。其工作原理是：每次循环中选择两个$alpha$进行优化处理。一旦找到一对合适的$alpha$，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个$alpha$必须要符合一定的条件，条件之一就是这两个$alpha$必须要在间隔边界之外，而其第二个条件则是这两个$alpha$还没有进行过区间化处理或者不在边界上。之所以要同时改变2个$alpha$，原因是前面的约束条件$\sum _{i=1}^n{a}_i{y}_i=0$，如果只是修改一个$alpha$，很可能导致约束条件失效。
  到目前为止，我们的 SVM 还比较弱，只能处理线性的情况，下面我们将引入核函数，进而推广到非线性分类问题

5. 优缺点
  优点
    1. 可以解决高维问题，即大型特征空间
    2. 解决小样本下机器学习问题
    3. 能够处理非线性特征的相互作用
    4. 无局部极小值问题
    5. 无需依赖整个数据
    6. 泛化能力比较强

  缺点
    1. 当观测样本很多时，效率并不是很高
    2. 对非线性问题没有通用解决方案，有时候很难找到一个合适的核函数
    3. 对于核函数的高维映射解释力不强，尤其是径向基函数
    4. 常规SVM只支持二分类
    5. 对缺失数据敏感

6. 主要应用
  文本分类、图像识别（主要二分类领域）

7. 实例分析
  上面看了那么多理论，现在我们开始动手实现一组数据处理。

• 准备数据

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 				#返回特征集与标签集

• $SMO$算法
    有了数据之后，我们利用简单SMO算法获取$\text{y}(x)=w^{T}x+b$的常量$b$和拉格朗日乘子$alpha$（$SMO$）。

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):   #简单SMO算法
    dataMatrix = mat(dataMatIn) 							#特征值转为矩阵
    labelMat = mat(classLabels).transpose() 				#标签值转为矩阵，且转置
    m, n = shape(dataMatrix)    							#特征值维度
    b = 0
    alphas = mat(zeros((m, 1))) 							#初始化m行1列的alpha向量全为0

    iter = 0
    while iter < maxIter:       							#迭代次数
        alphaPairsChanged = 0   							#记录alphas是否优化
        for i in range(m):
            #预测的类别 y[i] = w^T*x[i]+b; 其中因为 w = Σ(1~n) a[n]*label[n]*x[n]
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b   #预测结果
            Ei = fXi - float(labelMat[i])   				#误差  预测结果 - 真实结果
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):   #误差较大，需要优化，正常值在(0~C)
                j = selectJrand(i, m)   					#选择第二个随机alpha[j]
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b   #预测结果
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()

                #利用L和H将alpha调整到0~C
                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

                #序列最小优化算法计算alpha[j]的最优值
                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]值
                alphas[j] = clipAlpha(alphas[j], H, L)      #对L和H进行调整
                if (abs(alphas[j] - alphaJold) < 0.00001):  #如果改变幅度较小，无需继续优化
                    print("j not moving enough")
                    continue
                    
				#设置常数b
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                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):	#检查alpha值是否做了更新，如果在更新则将iter设为0后继续运行程序，否则退出
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)

    return b, alphas			#返回模型常量值和拉格朗日因子

  算法中使用到的方法：

def selectJrand(i, m):  #选取随机数，i为alpha下标，m为alpha数目
    j = i
    while (j == i):
        j = int(random.uniform(0, m))    #产生一个0~m的随机数
    return j

def clipAlpha(aj, H, L):	#调整目标值
    if aj > H:  			#目标值大于最大值
        aj = H
    if L > aj:  			#目标值小于最小值
        aj = L
    return aj   			#返回目标值

1. 可视化
     我们定义的超平面表达式为$\text{y}(x)=w^{T}x+b$，所以需要先把$w^T$计算出来。上面推导我们已经算出$w=\sum _{i=1}^n{a}_i{x}_i{y}_i$。

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

  有了$w^T$和$b$，现在就可以在图表中将我们的超平面显示出来。

def plotfig_SVM(xMat, yMat, ws, b, alphas):
    xMat = mat(xMat)
    yMat = mat(yMat)

    # b原来是矩阵，先转为数组类型后其数组大小为（1,1），所以后面加[0]，变为(1,)
    b = array(b)[0]
    fig = plt.figure()
    ax = fig.add_subplot(111)

    ax.scatter(xMat[:, 0].flatten().A[0], xMat[:, 1].flatten().A[0])

    x = arange(-1.0, 10.0, 0.1)     #x最大值，最小值根据原数据集dataArr[:, 0]的大小而定
    y = (-b-ws[0, 0]*x)/ws[1, 0]    #根据x.w + b = 0 得到，其式子展开为w0.x1 + w1.x2 + b = 0, x2就是y值
    ax.plot(x, y)

    for i in range(shape(yMat[0, :])[1]):
        if yMat[0, i] > 0:
            ax.scatter(xMat[i, 0], xMat[i, 1], color='blue')
        else:
            ax.scatter(xMat[i, 0], xMat[i, 1], color='green')

    # 找到支持向量，并在图中标红
    for i in range(100):
        if alphas[i] > 0.0:
            ax.scatter(xMat[i, 0], xMat[i, 1], color='red')
    plt.show()

6. 完整$SMO$算法
  从上图可以看到，我们已经成功将蓝绿两种点区分开，并且超平面在到两边的支持向量上都是最大的间隔。但是在几百个点组成的小规模数据集上，这种简单的SMO算法是没有什么问题，但是在更大的数据集上的运行速度就会变慢。下面我们就讨论一下完整的$SMO$算法，在这两个版本中，实现$alpha$的更改和代数运算的优化环节一模一样。但完整$SMO$选择alpha的方式却不同。
  完整$SMO$算法是通过一个外循环来选择第一个$alpha$值的，并且其选择过程会在两种方式之间进行交替：一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描。而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。对整个数据集的扫描相当容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历。同时，该步骤会跳过那些已知的不会改变的alpha值。在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中，我们会在选择j之后计算错误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。
  分析$SMO$

8. 核函数
  $SMO$算法对线性可分的数据具有非常好的效果，但是如果特征多了，维度更多，它也没有更好的方法，这里我们就需要引入一种叫做核函数的方法，它可以将数据从一个特征空间映射到另一个特征空间。所以当我们碰到线性不可分的数据时，可以将其映射到高维空间。
  $SVM$优化中一个特别好的地方就是，所有的运算都可以写成内积的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以把内积运算替换成核函数，而不必做简化处理。
  假设$X$是输入空间，$H$是特征空间，存在一个映射$\varphi$，使得$X$中的点$x$能够计算得到$H$空间中的点，即$h=\varphi (x)$。同理如果$x$，$z$是$X$空间中的点。函数$k(x,z)$满足条件：
$$k(x,z)=\varphi (x)\cdot \varphi (z)$$  则称k为核函数，而ϕ为映射函数。

  常见核函数：
    线性核函数  ：$K(x,y)=x^{T}y+c$
    多项式核函数 ：$K(x,y)=(a{x}^{T}y+c{)}^d$
    高斯核函数  ：$K(x,z)=\exp (-\frac{||x-z|{|}^2}{2{\sigma }_{}^2})$
    拉普拉斯核函数 ：$K(x,z)=\exp (-\frac{||x-z||}{{\sigma }_{}^{}})$
    $Sigmoid$核函数：$K(x,z)=\tanh (\beta x_i^{T}z+\theta )$
  其中高斯核函数最常用，可以将数据映射到无穷维，也叫做径向基函数，是某种沿径向对称的标量函数。这里我们主要讨论高斯核函数。如果$x$和$z$很相近，$||x-z||$趋近于0，那么核函数值为1，如果$x$和$z$相差很大，$||x-z||$趋近于无穷大，那么核函数值约等于0，因此它能够把原始特征映射到无穷维，当数据被映射到高维之后，就可以直接使用高斯核了。
  高斯核函数转换：

def kernelTrans(X, A, kTup):
    m, n = shape(X)
    K = mat(zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T
    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)) 	#径向基函数的高斯版本
    else:
        raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
        
    return K

未完待续

参考资料
  支持向量机通俗导论
  机器学习实战