机器学习:支持向量机(SVM)

优点:泛化错误率低,计算开销不大,结果易解释
缺点:对参数调节和核函数的选择敏感,原始分类器不加修改仅适用于处理二分类问题

线性可分

  • 给定一个二维数据集,如果可以用一条直线,将数据分成两类,在直线的一边都是一种分类,另一边的都是另一种分类,我们说这个数据集线性可分
  • 该直线称为分隔超平面(Separating Hyperplane),也就是分类的决策边界
  • 同理对 N 维数据集,如果存在可对数据进行分类的 N-1 维超平面,我们说该数据集线性可分
  • 采用这种方式来构建分类器,如果数据点离决策边界越远,那么其最后的预测结果也就越可信

支持向量

  • 数据集可能存在多个可用于分类的超平面
  • 我们希望找到离分隔超平面最近的点,确保它们离分隔面的距离尽可能远
  • 符合这个条件的超平面鲁棒性(Robust)最好,对未知数据的分类准确性高
  • 在这里距离又被称为间隔(Margin)
  • 离分隔超平面最近的那些点就是支持向量(Support Vector)
  • SVM(Support Vector Mechine)就是要最大化支持向量到分隔面的距离

寻找最大间隔

下图以有两个特征的二维数据为例子
  
     
  
离分割超平面最近的点就是支持向量,对于最佳分割超平面而言,两边的支持向量到分割超平面的距离是相等的,并且该距离等于过支持向量且与分割面平行的面到分割面的距离
  
分隔超平面可表示为
  W1X1+...+WjXj+...+WmXm+b=0
  
可表示为
  WX+b=0  或  WTX+b=0
  
平行面的 W 系数有相同的比例,所以两个平行面可表示为
  
  α(W1X1+...+WjXj+...+WmXm)+c=0
  β(W1X1+...+WjXj+...+WmXm)+d=0
  
对平面的 Wb 系数进行同比例的缩放,表示的依然是同一个平面,所以可改写为
  
  W1X1+...+WjXj+...+WmXm+c/α=0
  W1X1+...+WjXj+...+WmXm+d/β=0
  
由于平面等距,所以有
  
  bc/α=d/βb=e
  
进一步改写为
  
  W1X1+...+WjXj+...+WmXm+be=0
  W1X1+...+WjXj+...+WmXm+b+e=0
  
三个方程都除以 e 得到
  
  (W1X1+...+WjXj+...+WmXm+b)/e=0
  (W1X1+...+WjXj+...+WmXm+b)/e=1
  (W1X1+...+WjXj+...+WmXm+b)/e=1
  
W/eb/e 改用 Wb 表示,最终得到
  
  WTX+b=0
  WTX+b=1
  WTX+b=1
  
之所以取 1 是为了方便计算距离,平行线 ax+by+c=0ax+by+d=0 的距离为
  
  |cd|a2+b2
  
所以这三个面的距离可表示为
  
  1W12+...+Wj2+...+Wm2=1||W||
    
满足 WTX+b>0 的分类记为 y=1
满足 WTX+b<0 的分类记为 y=1
  
由于样本数据最近的点在 WTX+b=±1 上,样本数据可以进一步表示为
  
满足 WTX+b>1 的分类记为 y=1
满足 WTX+b<1 的分类记为 y=1
    
  
这样 SVM 的训练过程就是要寻找合适的 Wb 以满足
  
  max1||W||
  
  s.t.  yi(WTxi+b)1, i=1,2,...,n  (s.t. 是约束条件的意思)


max1||W|| 等价于 min12||W||2 (取为 12 和平方的目的是为了后面方便求解)


所以优化目标又可以写为
  
  min12||W||2
  
  s.t.  yi(WTxi+b)1, i=1,2,...,n


拉格朗日乘数法
  
  给定函数 y=f(X) 和约束条件 hi(X)=0,为求解函数极值,构建拉格朗日函数
  
    L=f(X)+i=1kαihi(X)
  
  对 Xαi 求导并令导数等于 0 即
  
    xL=0
    αL=0
  
  求出的 X 带入 f(X) 可得到极值
  
  如果约束条件是不等式 hi(X)0
  只要满足 KKT 条件(这里不推导了)也可以用拉格朗日数法


于是前面的优化目标的拉格朗日函数为
  
  L=12||W||2+i=1nαi(1yi(WTxi+b))
  
Wb 求导并令导数等于 0
  
  LWj=12(W12+...+Wj2+...+Wm2)Wji=1nαiyiWjxijWj
  
     =Wji=1nαiyixij
  
     =0


   Lb=i=1nαiyibb
  
     =i=1nαiyi
  
     =0
  
将结果带入 L 可得
  
  L=12||W||2+i=1nαi(1yi(WTxi+b))
  
    =12(W12+...+Wm2)+i=1nαii=1nαiyi(WTxi)i=1nαiyib
  
    =12(W12+...+Wm2)+i=1nαii=1nαiyi(WTxi)
  
    =i=1nαi+j=1m12Wj2i=1nαiyi(WTxi)
  
    =i=1nαi+j=1m12(i=1nαiyixij)(i=1nαiyixij)i=1nαiyij=1m(Wjxij)
  
    =i=1nαi+j=1m12(i=1nαiyixij)(i=1nαiyixij)i=1nαiyij=1m((i=1nαiyixij)xij)
  
    =i=1nαi12i=1nj=1nαiαjyiyjxiTxj
  
  s.t.  i=1nαiyi=0, i=1,2,...,n
  
  s.t.  αi0, i=1,2,...,n
  
求出 α 后,将 Wj=i=1nαiyixij 代入 f(X)=WTX+b 得到点 x 的值为
  
  f(X)=WTX+b
  
      =(i=1nαiyixi1)x1+...+(i=1nαiyixij)xj+...+(i=1nαiyixnm)xm+b
  
      =i=1nαiyi(xi1x1+...+xijxj+...+ximxm)+b
  
      =i=1nαiyi(xiTx)+b
  
  或者
  
  f(X)=WTX+b
  
      =(i=1nαiyixi1)x1+...+(i=1nαiyixij)xj+...+(i=1nαiyixnm)xm+b
  
      =a1y1j=1mx1jxj+...+aiyij=1mxijxj+...+anynj=1mxnjxj+b
  
      =[a1y1 ... aiyi ... anyn]×([x11...x1j...x1mxi1...xij...ximxn1...xnj...xnm]×[x1xjxm])+b
  
      =(ay)T×(X×xT)+b
  
  其中
    ay(n,1) 矩阵
    X(n,m) 矩阵,是样本数据集
    x(1,m) 矩阵,是要求解的一条数据
  
到此为止都是假设数据线性可分,考虑到数据不一定完全线性可分的时候,可以引入一个常数限制 α,允许部分数据划分错误,即 α 的约束条件变为
  
  s.t.  i=1nαiyi=0, i=1,2,...,n
  
  s.t.  Cαi0, i=1,2,...,n
  
可以使用通用的二次规划算法求解 α 的值,但这会造成很大的开销,为了简化计算,人们提出了很多高效算法,SMO 是一个著名代表

SMO(Sequential Minimal Optimization,序列最小优化)

SMO 将大优化问题分解为多个小优化问题,这些小优化问题往往很容易求解,并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是一致的,而 SMO 算法的求解时间短很多
  
基本思路是先取两个值 αiαj 做变量,然后固定其他 α 值,由于只有两个变量所以很容易优化,然后再继续取下一对 α 值,照同样的方法优化,就这样不断迭代,直到无法再优化,或者达到最大迭代次数
  
i=1nαiyi=0 可以得到,当固定一组 α 以外的其他值时,有
  
  αiyi+αjyj=e
其中
  e=ki,jαkyk
  e 为常数
  
α1α2 为例子,因为其他 α 值被当做常数,可以将优化目标中除 α1α2 外的项去掉(因为常数无法优化)
  
所以原优化目标
  
  L=i=1nαi12i=1nj=1nαiαjyiyjxiTxj
  
等价于
  
  L=α1+α212α12y12x1Tx112α22y22x2Tx2α1α2y1y2x1Tx2
  
     α1y1i=3nαiyix1Txiα2y2i=3nαiyix2Txi
  
   =α1+α212α12x1Tx112α22x2Tx2α1α2y1y2x1Tx2
  
     α1y1i=3nαiyix1Txiα2y2i=3nαiyix2Txi
  
再用 α1=eα2y2y1 代入得到
  
  L=eα2y2y1+α212(eα2y2)2x1Tx112α22x2Tx2(eα2y2)α2y2x1Tx2
  
     (eα2y2)i=3nαiyix1Txiα2y2i=3nαiyix2Txi
  
注意有 1=y2 以及 y=1y,对 α2 求导得到

  Lα2=y1y2+y22+(eα2y2)x1Tx1y2α2x2Tx2ey2x1Tx2+2α2y22x1Tx2
  
      +y2i=3nαiyix1Txiy2i=3nαiyix2Txi
  
     =y1y2+y22+ex1Tx1y2α2x1Tx1α2x2Tx2ey2x1Tx2+2α2x1Tx2
  
      +y2i=3nαiyix1Txiy2i=3nαiyix2Txi

令导数为 0 得到

  α2(x1Tx1+x2Tx22x1Tx2)
  
  =y2(y2y1+ex1Tx1ex1Tx2+i=3nαiyix1Txii=3nαiyix2Txi)
  
  =y2(y2y1+ex1Tx1ex1Tx2
  
     +(f(x1)i=12αiyix1Txib)(f(x2)i=12αiyix2Txib))
  
  =y2((f(x1)y1)(f(x2)y2)+ex1Tx1ex1Tx2i=12αiyix1Txi+i=12αiyix2Txi)
  
  =y2((f(x1)y1)(f(x2)y2)+(α1y1+α2y2)x1Tx1(α1y1+α2y2)x1Tx2
  
     α1y1x1Tx1α2y2x1Tx2+α1y1x2Tx1+α2y2x2Tx2)
  
  =y2((f(x1)y1)(f(x2)y2)+α2y2(x1Tx1+x2Tx22x1Tx2))
  
  =y2((f(x1)y1)(f(x2)y2))+α2(x1Tx1+x2Tx22x1Tx2)

这过程中,代入 α1y1+α2y2=ef(x) 时用到的 α (即等式右边用的)是旧的

可以得到优化后的值为
  
  α2new=α2old+y2((f(x1)y1)(f(x2)y2))x1Tx1+x2Tx22x1Tx2

        =α2old+y2(E1E2)x1Tx1+x2Tx22x1Tx2
  
然后计算优化后的 α1
  
  α1new=α1old+y1y2(α2oldα2new)


得到优化值后,继续取下一对 α 值进行优化
  
  
考虑到 Cα0α1y1+α2y2=e 的限制,在优化的过程中,α 值是受到限制的
  
y1y2 异号,比如 y1=1y2=1 时,有 α1α2=e
  
  当 e>Ce<C
  
    由于有 Cα0,该式子无解
  
  当 Ce0
  
    α2 的取值范围在 (0,Ce)
  
  当 0eC
  
    α2 的取值范围在 (e,C)
  
  可以看到 α2 的取值范围是
  
    L=max(0,e)H=min(C,Ce)
  
  即有 α2new  的取值范围是
  
    max(0, α2oldα1old)α2newmin(C, C+α2oldα1old)
  
同理可求得 y1y2 同号时
  
  max(0, α2old+α1oldC)α2newmin(C, α2old+α1old)
  
  
接下来求 b
  
  由 f(x)=i=1nαiyi(xiTx)+b 可以反推

  If: 0<a1new<C ,通过第一个点计算
  
    b1new=y1i=3n(αiold)yi(xiTx1)(α1new)y1(x1Tx1)(α2new)y2(x2Tx1)
  
    又有 E1=i=3n(αiold)yi(xiTx1)+(α1old)y1(x1Tx1)+(α2old)y2(x2Tx1)+boldy1
  
    b1new=y1E1+(α1old)y1(x1Tx1)+(α2old)y2(x2Tx1)+boldy1
  
          (α1new)y1(x1Tx1)(α2new)y2(x2Tx1)
  
         =E1y1(x1Tx1)(α1newα1old)y2(x2Tx1)(α2newα2old)+bold
  
  elif: 0<a2new<C,通过第二个点计算

    b2new=E2y2(x1Tx2)(α1newα1old)y2(x2Tx2)(α2newα2old)+bold
  
  else :
    a1newa2new 或者是 0 或者是 C
    
    bnew=b1new + b2new2

分类预测

前面已经给出了公式,这里再总结一下
  
  Wj=i=1nαiyixij
  
  f(X)=WTX+b
  
    或者
  
  f(x)=i=1nαiyi(xiTx)+b
  
    或者
    
  f(x)=(ay)T×(X×xT)+b
  
  
  αi!=0 对应的 xi 就是支持向量(??)

简化版 SMO 代码

# coding=utf-8
import numpy as np


def selectJrand(i, n):
    """
    在 (0, n) 之间随机选择一个不等于 i 的 值
    """
    j = i
    while j == i:
        j = int(np.random.uniform(0, n))
    return j


def clipAlpha(aj, H, L):
    """
    限制 aj 的值在 L 和 H 之间
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    """
    dataMatIn   - 样本数据集,(n, m) 维数组,n 是样本数,m 是特征数
    classLabels - 标签,(1, n) 维数组,取值为 1 或 -1
    C           - 使 C > a > 0,以处理不 100% 线性可分的情况,加上 C 的限制后不会一直增加 a 值,同时允许部分数据错分
    toler       - 能容忍的误差范围
    maxIter     - 迭代次数
    """

    # 样本转为 numpy 矩阵
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()

    # 样本数,特征数
    n, m = np.shape(dataMatrix)

    # 初始化 a 和 b 为 0
    b = 0
    alphas = np.mat(np.zeros((n, 1)))

    # 如果连续 maxIter 迭代 a 值都没改变就跳出
    # 如果 a 值有变化要重新计算迭代次数
    iter = 0
    while iter < maxIter:
        alphaPairsChanged = 0

        # 遍历每个样本
        for i in range(n):
            # f = WX + b
            #   = (a . y)^T * (X * x^T) + b
            fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b

            # 计算误差
            Ei = fXi - float(labelMat[i])

            # 误差偏大,看能不能优化 a 值
            if (labelMat[i] * Ei < -toler and alphas[i] < C) or (labelMat[i] * Ei > toler and alphas[i] > 0):
                # 随机选择另一个样本
                j = selectJrand(i, n)

                # 计算误差
                fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])

                # 保存对应的 a 值
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()

                # 计算优化后的 a2 的取值范围
                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

                # a2-new = a2-old - y2(E1-E2)/(2 * x1 * x2^T - x1 * x1^T - x2 * x2^T)
                #        = a2-old - y2(E1-E2)/eta
                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

                # 优化 aj
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta

                # 限制 aj
                alphas[j] = clipAlpha(alphas[j], H, L)

                if abs(alphas[j] - alphaJold) < 0.00001:
                    # 没优化
                    print "j not moving enough"
                    continue

                # 优化 ai
                # a1-new = a1-old + y1*y2*(a2-old - a2-new)
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])

                # 计算 b 值
                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):
    """
    计算 W 系数
    """
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(X)
    w = np.zeros((n, 1))

    for i in range(m):
        w += np.multiply(alphas[i]*labelMat[i], X[i, :].T)

    return w

完整版 SMO

  • 简化版 SMO 运算速度慢
  • 接收参数一样
  • 对 alphas 的更改运算一样
  • 对第二个 alphas 的选择不一样,简化版是随机的,完整版复杂一点
  • 外循环最多循环一次,不会出现重置 iter 的情况,而且两次迭代如果 alphas 都没有优化就会跳出外循环

核函数

  • 有的数据是非线性可分的,比如某些二维数据无法用直线划分而是用圆划分
  • 可以将数据从原始空间映射到更高维度的特征空间(比如讲二维数据映射到三维空间),使得数据在这个新的特征空间内线性可分
  • 对于有限维空间,一定存在这样一个高维空间使得数据线性可分
  • 用于映射的函数就被称为核函数,很多其他的机器学习算法也都用到核函数
  • 假设映射后的特征向量为 Φ(x),则划分超平面表示为f(x)=wTΦ(x)+b
  • 参数的推导过程和 SMO 一样

高斯核函数

径向基函数是 SVM 中常用的一个核函数。
径向基函数是一个采用向量作为自变量的函数,能够基于向量距离运算输出一个标量。
径向基函数的高斯版本的高斯版本,函数为
  
  k(xi,xj)=exp(||xixj||22σ2)
  
其中
  exp() 表示自然数 e 的指数
  σ是函数值跌落到 0 的速度参数
  
该函数将数据映射到更高维的空间
具体来说是映射到一个无穷维的空间(不理解,貌似是按泰勒展开的话有无穷维)
  
优化的目标函数由
  
  L=i=1nαi12i=1nj=1nαiαjyiyjxiTxj
  
改写为
  
  L=i=1nαi12i=1nj=1nαiαjyiyjk(xi, xj)
  
训练过程和前面的 SMO 一样,推导过程中的 xiTxj 改为 k(xi, xj) 就可以
  
高斯核函数计算

def kernel(X, omiga):
    """
    用高斯核函数对样本数据 X 转换

    X 是二维矩阵

    返回矩阵 K

    K(i, j) = exp(-(||Xi - Xj||**2)/(2 * (omiga**2)))
    """

    m, n = np.shape(X)

    K = np.mat(np.zeros((m, m)))

    for i in range(m):
        for j in range(m):
            # 计算 ||Xi - Xj||**2
            deltaRow = X[j, :] - X[i, :]
            d = deltaRow * deltaRow.T

            K[i, j] = np.exp(d / (-2 * omiga ** 2))

    return K
 
 


posted @   moon~light  阅读(469)  评论(0编辑  收藏  举报
编辑推荐:
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
阅读排行:
· 10年+ .NET Coder 心语 ── 封装的思维:从隐藏、稳定开始理解其本质意义
· 提示词工程——AI应用必不可少的技术
· 地球OL攻略 —— 某应届生求职总结
· 字符编码:从基础到乱码解决
· SpringCloud带你走进微服务的世界
点击右上角即可分享
微信分享提示