SVM-SMO

关于𝛼向量的函数

• SVM优化目标函数
  
  KKT条件：
  𝑦𝑖(𝑤∙𝜙(𝑥𝑖)+𝑏)≥1−𝜉𝑖
  𝜇𝑖≥0,𝛼𝑖≥0
  𝑤=∑𝛼𝑖𝑦𝑖𝜙(𝑥𝑖)，∑𝛼𝑖𝑦𝑖=0， 𝐶−𝛼𝑖−𝜇𝑖=0（0≤𝛼𝑖≤𝐶）
  𝛼𝑖(𝑦𝑖(𝑤𝑇𝜙(𝑥𝑖)+𝑏)−1+𝜉𝑖)=0，𝜇𝑖𝜉𝑖=0（(𝐶−𝛼𝑖)𝜉𝑖 = 0）

• 则有：
  𝛼𝑖=0 ⇒𝜉𝑖=0 ⇒𝑦𝑖(𝑤∙𝜙(𝑥𝑖)+𝑏)≥1 ⇒𝑦𝑖𝑔(𝑥𝑖) ≥ 1
  0<𝛼𝑖<𝐶 ⇒ 𝑦𝑖(𝑤𝑇𝜙(𝑥𝑖)+𝑏)−1+𝜉𝑖=0， 𝜉𝑖=0 ⇒ 𝑦𝑖(𝑤∗∙𝜙(𝑥𝑖)+𝑏)=1 ⇒ 𝑦𝑖𝑔(𝑥𝑖)=1
  𝛼𝑖=𝐶 ⇒𝑦𝑖(𝑤𝑇𝜙(𝑥𝑖)+𝑏)−1+𝜉𝑖=0 ⇒ 𝑦𝑖(𝑤∗∙𝜙(𝑥𝑖)+𝑏)≤1 ⇒𝑦𝑖𝑔(𝑥𝑖) ≤ 1

SMO算法

1. 每次只优化两个变量，将其他的变量都视为常数。由于∑𝛼𝑖𝑦𝑖=0.假如将𝛼3,𝛼4,…,𝛼𝑚　固定，那么𝛼1,𝛼2之间的关系也确定了。

2. 定义𝐾𝑖𝑗=𝜙(𝑥𝑖)∙𝜙(𝑥𝑗)，去除常量后目标函数：
   
   C是协调软边界的正则化惩罚系数。

3. 根据约束条件𝛼1𝑦1+𝛼2𝑦2=k；0≤𝛼𝑖≤𝐶 𝑖=1,2。又由于𝑦1,𝑦2均只能取值1或者-1, 这样𝛼1,𝛼2在[0,C]和[0,C]形成的盒子里面，并且两者的关系直线的斜率只能为1或者-1。
   

4. 我们采用启发式的迭代法，假设我们上一轮迭代得到的解是𝛼𝑜𝑙𝑑1,𝛼𝑜𝑙𝑑2，假设求解得到的𝛼2是𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2，本轮迭代完成后的解为𝛼𝑛𝑒𝑤1,𝛼𝑛𝑒𝑤2。
   由于𝛼𝑛𝑒𝑤2必须满足上图中的约束。假设L和H分别是上图中𝛼𝑛𝑒𝑤2所在的线段的边界。
   𝐿≤𝛼𝑛𝑒𝑤2≤𝐻
   左图：
   𝐿=𝑚𝑎𝑥(0,𝛼𝑜𝑙𝑑2−𝛼𝑜𝑙𝑑1) 𝐻=𝑚𝑖𝑛(𝐶,𝐶+𝛼𝑜𝑙𝑑2−𝛼𝑜𝑙𝑑1)
   右图：
   𝐿=𝑚𝑎𝑥(0,𝛼𝑜𝑙𝑑2+𝛼𝑜𝑙𝑑1−𝐶) 𝐻=𝑚𝑖𝑛(𝐶,𝛼𝑜𝑙𝑑2+𝛼𝑜𝑙𝑑1)
   我们通过求导得到的𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2，

   1. 𝛼𝑛𝑒𝑤2 = 𝐻 （𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2 >𝐻）
   2. 𝛼𝑛𝑒𝑤2 = 𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2（𝐿< 𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2 <𝐻）
   3. 𝛼𝑛𝑒𝑤2 = 𝐿（𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2 <𝐿）

5. 求解𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2 ，令
   
   则目标优化函数：
   
   𝛼1𝑦1+𝛼2𝑦2=𝜍 ⇒ 𝛼1=𝑦1(k−𝛼2𝑦2)

   带入目标优化函数，消除𝛼1，并对𝛼2求偏导值为0化简：

   (𝐾11+𝐾22−2𝐾12)𝛼2=𝑦2(𝑦2−𝑦1+𝜍𝐾11−𝜍𝐾12+𝑣1−𝑣2)

   令
   将𝜍=𝛼1𝑦1+𝛼2𝑦2带入上式化简得：

   (𝐾11+𝐾22−2𝐾12)𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2 = (𝐾11+𝐾22−2𝐾12)𝛼𝑜𝑙𝑑2 + 𝑦2(𝐸1−𝐸2)

   𝛼𝑛𝑒𝑤,𝑢𝑛𝑐2=𝛼𝑜𝑙𝑑2 + 𝑦2(𝐸1−𝐸2) / (𝐾11+𝐾22−2𝐾12)

6. 得出𝛼𝑛𝑒𝑤2，根据𝛼𝑛𝑒𝑤1+𝛼𝑛𝑒𝑤2 = 𝛼𝑜𝑙𝑑1 + 𝛼𝑜𝑙𝑑2 求得 𝛼𝑛𝑒𝑤1

两个变量的选择

第一个变量的选择

• SMO算法称选择第一个变量为外层循环，这个变量需要选择在训练集中违反KKT条件最严重的样本点。
  一般来说，首先选择违反0<𝛼∗𝑖<𝐶 ⇒ 𝑦𝑖𝑔(𝑥𝑖)=1这个条件的点。如果这些支持向量都满足KKT条件，再选择违反𝛼∗𝑖=0⇒𝑦𝑖𝑔(𝑥𝑖)≥1 和 𝛼∗𝑖=𝐶⇒𝑦𝑖𝑔(𝑥𝑖)≤1的点。

第二个变量的选择

• SMO算法称选择第二个变量为内层循环，假设在外层循环已经找到𝛼1, 第二个变量𝛼2的选择标准是让|𝐸1−𝐸2|有足够大的变化。
• 由于𝛼1定了的时候,𝐸1也确定了，所以要想|𝐸1−𝐸2|最大，只需要在𝐸1为正时，选择最小的𝐸𝑖作为𝐸2， 在𝐸1为负时，选择最大的𝐸𝑖作为𝐸2，可以将所有的𝐸𝑖保存下来加快迭代。
• 如果内存循环找到的点不能让目标函数有足够的下降， 可以采用遍历支持向量点来做𝛼2,直到目标函数有足够的下降， 如果所有的支持向量做𝛼2都不能让目标函数有足够的下降，可以跳出循环，重新选择𝛼1

b和差值𝐸𝑖

1. 在每次完成两个变量的优化之后，需要重新计算b
2. 当0<𝛼𝑛𝑒𝑤1<𝐶时
3. 合并上面两式：
   
   同理，当0<𝛼𝑛𝑒𝑤2<𝐶时:
   
4. 𝑏𝑛𝑒𝑤 =（𝑏𝑛𝑒𝑤1+𝑏𝑛𝑒𝑤2）/ 2
   
   S是所有支持向量𝑥𝑗的集合(0< 𝛼𝑖 ≤ 𝐶)。

SMO算法小结

输入是m个样本(𝑥1,𝑦1),(𝑥2,𝑦2),…,(𝑥𝑚,𝑦𝑚)，其中x为n维特征向量。
y为二元输出，值为1或-1；精度e，输出为近似解𝛼

1. 取初值𝛼0=0，𝑘=0(第k次迭代)
2. 选择𝛼𝑘1，𝛼𝑘2，求出新的𝛼𝑛𝑒𝑤(k+1)
3. 计算𝑏𝑘+1和𝐸𝑖
4. 在精度e范围内检查是否满足终止条件：
   ∑𝛼𝑖𝑦𝑖=0
   0≤𝛼𝑖≤𝐶，𝑖=1,2…𝑚
   𝛼𝑖=0 ⇒𝑦𝑖𝑔(𝑥𝑖)≥1
   0<𝛼𝑖<𝐶 ⇒ 𝑦𝑖𝑔(𝑥𝑖)=1
   𝛼𝑖=𝐶 ⇒ 𝑦𝑖𝑔(𝑥𝑖)≤1
5. 如果满足则结束，返回𝛼𝑘+1,否则转到步骤2）

#选择ai，aj
def selecta(i,m):
    j=i
    while j==i:
        j=random.randint(0,m-1)
    return j
#a上下限
def clipalpha(aj,h,l):
    if aj>h:
        aj=h
    if aj<l:
        aj=l
    return aj
#smo算法获取alpha和b
def smo(xmat,ymat,c,toler,maxiter):
    b=0
    m,n=xmat.shape
    alpha=np.mat(np.zeros((m,1)))
    iters=0
    while iters < maxiter:
        alpha_=0
        for i in range(m):
            #计算i误差
            fxi=(np.multiply(alpha,ymat).T*xmat)*xmat[i,:].T+b
            ei=fxi-ymat[i]
            #选择违反kkt条件(当0<alpha<c)的点进行优化
            if ((ei*ymat[i]<-toler) and (alpha[i]<c)) or ((ei*ymat[i]>toler) and (alpha[i]>0)):
                j=selecta(i,m)
                #计算j误差
                fxj=(np.multiply(alpha,ymat).T*xmat)*xmat[j,:].T+b
                ej=fxj-ymat[j]
                #保存old a
                alphaiold=alpha[i].copy()
                alphajold=alpha[j].copy()
                #计算上下限
                if ymat[i] != ymat[j]:
                    high=min(c,c+alpha[j]-alpha[i])
                    low=max(0,alpha[j]-alpha[i])
                if ymat[i] == ymat[j]:
                    high=min(alpha[j]+alpha[i],c)
                    low=max(alpha[j]+alpha[i]-c,0)
                if high==low:
                    continue
                eta=xmat[i]*xmat[i].T+xmat[j]*xmat[j].T-2*xmat[i]*xmat[j].T
                if eta<=0:
                    continue
                #更新aj
                alpha[j]=alphajold+(ymat[j]*(ei-ej))/eta
                alpha[j]=clipalpha(alpha[j],high,low)
                if abs(alpha[j]-alphajold)<0.00001:
                    continue
                #更新ai
                alpha[i]=alphaiold+ymat[i]*ymat[j]*(alphajold-alpha[j])
                #
                b1=b-ei-ymat[i]*(xmat[i]*xmat[i].T)*(alpha[i]-alphaiold)-ymat[j]*(xmat[i]*xmat[j].T)*(alpha[j]-alphajold)
                b2=b-ej-ymat[i]*(xmat[i]*xmat[j].T)*(alpha[i]-alphaiold)-ymat[j]*(xmat[j]*xmat[j].T)*(alpha[j]-alphajold)
                if (alpha[i]>0) and (alpha[i]<c):
                    b=b1
                elif (alpha[j]>0) and (alpha[j]<c):
                    b=b2
                else:
                    b=(b1+b2)/2
                alpha_+=1
        if alpha_==0:
            iters+=1
        else:
            iters=0
    return b,alpha
#获取支持向量
def get_sv(xmat,ymat,alpha):
    m=xmat.shape[0]
    sv_x=[]
    sv_y=[]
    for i in range(m):
        if alpha[i]>0:
            sv_x.append(xmat[i])
            sv_y.append(ymat[i])
    sv_x=np.array(sv_x).T
    sv_y=np.array(sv_y).T
    return sv_x,sv_y