SVM Python实现

Python实现SVM的理论知识

• SVM原始最优化问题:

minw,b,ξ12||w||2+C∑i=1mξ(i)minw,b,ξ12||w||2+C∑i=1mξ(i)

s.t.  y(i)(wTx(i)+b),i=1,2,...,mξ(i)≥0,i=1,2,...ms.t.  y(i)(wTx(i)+b),i=1,2,...,mξ(i)≥0,i=1,2,...m

• 原始问题转为对偶问题

minα12∑i=1m∑j=1mα(i)α(j)y(i)y(j)K(x(i),x(j))−∑i=1mα(i)minα12∑i=1m∑j=1mα(i)α(j)y(i)y(j)K(x(i),x(j))−∑i=1mα(i)

s.t.  ∑i=1mα(i)y(i)=00≤α(i)≤C,i=1,2,...,ms.t.  ∑i=1mα(i)y(i)=00≤α(i)≤C,i=1,2,...,m

• 对于第2点, 设α∗=(α∗1,α∗2,...,α∗m)α∗=(α1∗,α2∗,...,αm∗), 若存在α∗α∗的一个分量α∗j,0<α∗j<Cα∗j,0<αj∗<C, 则原始问题的解w∗,b∗w∗,b∗为

w∗=∑i=1mα∗iy(i)x(i)b∗=y(j)−∑i=1my(i)α(i)K(x(i),x(j))w∗=∑i=1mαi∗y(i)x(i)b∗=y(j)−∑i=1my(i)α(i)K(x(i),x(j))

• SMO算法中使用到的公式(公式用使用下标1和2不是指在样本中的第1个样本和第2个样本, 而是指第1个参数和第2个参数, 在编程中我们使用i和j替代, i和j是在X输入样本中的样本下标)

  E(i)=g(x(i))−y(i)其中g(x(i))=∑i=1mα(i)y(i)K(x(i),x(j))+b所以E(i)=(∑i=1mα(i)y(i)K(x(i),x(j))+b)−y(i)   where  i=1,2E(i)=g(x(i))−y(i)其中g(x(i))=∑i=1mα(i)y(i)K(x(i),x(j))+b所以E(i)=(∑i=1mα(i)y(i)K(x(i),x(j))+b)−y(i)   where  i=1,2

  
  

  
  

  αnew,unc2=αold2+y2(E(i)−E(j))η其中η=K11+K22−2K12注意:K11指的是在使用核函数映射之后的输入样本中的第i行与第j行样本,同理K22指的是在使用核函数映射之后的输入样本中的第i行与第j行样本...注意:η不能为小于等于0,如果出现这种情况,则在迭代函数中直接返回0α2new,unc=α2old+y2(E(i)−E(j))η其中η=K11+K22−2K12注意:K11指的是在使用核函数映射之后的输入样本中的第i行与第j行样本,同理K22指的是在使用核函数映射之后的输入样本中的第i行与第j行样本...注意:η不能为小于等于0,如果出现这种情况,则在迭代函数中直接返回0

  
  

  if (yi * Ei > toler and alphai > 0) or (yi * Ei < -toler and alphaj < 0):    # 违反了
      passelse:    # 没有违反KKT, 直接返回0
      pass

  • 如果y(1)≠y(2)y(1)≠y(2)则
    

    L=max(0,αold2−αold1)H=min(C,C+αold2−αold1)L=max(0,α2old−α1old)H=min(C,C+α2old−α1old)

    
    

  • 如果y(1)=y(2)y(1)=y(2)则
    

    L=max(0,αold2+αold1−C)H=min(C,αold2+αold1)L=max(0,α2old+α1old−C)H=min(C,α2old+α1old)

    
    

  • 注意: 得到的L与H不能相同, 如果相同则直接返回0

  • 定义函数clip_alpha(alpha, L, H)
    ```py

    def calc_alpha(alpha, L, H):
    if alpha > H:
    alpha = H
    elif alpha < L:
    alpha = L
    return alpha
    ```

  • 计算αnew1α1new

  • 计算出αnew1α1new之后, 比较abs(αnew1−αold1)abs(α1new−α1old)与我们规定的精度的值(一般零点几), 如果小于精度则直接返回0

  • 裁剪αnew,unc1α1new,unc-> 在clip_alpha函数中

    
    

    αnew1=αold1+y(1)y(2)(αold2−αnew2)α1new=α1old+y(1)y(2)(α2old−α2new)

    
    

  • 选择第2个α2α2, 选择依据就是让abs((i)−E(j))abs((i)−E(j))最大, 那个jj就是我们期望的值, 从而α2=αjα2=αj, 定义函数select_j(Ei, i, model)
py   # model封装了整个SVM公式中的参数(C, xi)与数据(X, y)   # 其中model还有一个E属性, E为一个mx2的矩阵, 初始情况下都为0, 如果E对应的alpha被访问处理过了, 就会在赋予model.E[i, :] = [1, Ei]   def select_j(Ei, i, model):       j = 0       max_delta_E = 0       Ej = 0       # 查找所有已经被访问过了样本       nonzeros_indice = nonzeros(model.E[:, 0])[0]       if len(nonzeros_indice) > 1:           # 在for循环中查找使得abs(Ei-Ej)最大的Ej和j           for index in nonzeros_indice:               # 选择的两个alpha不能来自同一个样本               if index == i:                    continue               E_temp = calc_E(i, model)               delta_E = abs(E_temp - Ei)               if delta_E > max_delta_E:                   delta_E = max_delta_E                   j = index                   Ej = E_temp       else:           j = i           while j = i:               Ej = int(random.uniform(0, model.m))       return j, Ej

  • α1α1是否违反了KKT条件

  • 计算αnew,unc2α2new,unc-> 在iter函数中

  • 计算E(i)E(i)-> 在calc_E函数中

• 使用SMO算法对对偶问题求解, 求出αα, 从而得出w,bw,b, 大致思路如下

  • 初始化αα向量为m×1m×1, 元素为0的向量, 一开始α(1)α(1)的选择没有之前的依据, 所以使用从第一个alphaalpha开始选取

  • 如果选入的αα没有违反KKT条件则跳过, 迭代下一个αα

  • 将选出的α(1)α(1)代入iter函数, 该函数的工作是根据当前α(1)α(1)选择出第二个α(2)α(2), 接着根据公式更新α(2),α(1)α(2),α(1), 公式在上面已经给出, 注意选出来的α2α2的在输入样本中不能与α1α1是同一个

  • 迭代完所有αα, 下一步就是找出满足支持向量条件的αα, 即0≤α≤C0≤α≤C, 再将他们迭代

Python实现SVM的代码

#!/usr/bin/env pythonfrom numpy import *class Model(object):

    def __init__(self, X, y, C, toler, kernel_param):
        self.X = X
        self.y = y
        self.C = C
        self.toler = toler
        self.kernel_param = kernel_param
        self.m = shape(X)[0]
        self.mapped_data = mat(zeros((self.m, self.m)))        for i in range(self.m):
            self.mapped_data[:, i] = gaussian_kernel(self.X, X[i, :], self.kernel_param)
        self.E = mat(zeros((self.m, 2)))
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0def load_data(filename):
    X = []
    y = []    with open(filename, 'r') as fd:        for line in fd.readlines():
            nums = line.strip().split(',')
            X_temp = []            for i in range(len(nums)):                if i == len(nums) - 1:
                    y.append(float(nums[i]))                else:
                    X_temp.append(float(nums[i]))
            X.append(X_temp)    return mat(X), mat(y)def gaussian_kernel(X, l, kernel_param): 
    sigma = kernel_param 
    m = shape(X)[0]
    mapped_data = mat(zeros((m, 1)))    for i in range(m):
        mapped_data[i] = exp(-sum((X[i, :] - l).T * (X[i, :] - l) / (2 * sigma ** 2)))    return mapped_datadef clip_alpha(L, H, alpha):
    if alpha > H:
        alpha = H    elif alpha < L:
        alpha = L    return alphadef calc_b(b1, b2):
    return (b1 + b2) / 2def calc_E(i, model):
    yi = float(model.y[i])
    gxi = float(multiply(model.alphas, model.y).T * model.mapped_data[:, i] + model.b)
    Ei = gxi - yi    return Eidef select_j(Ei, i, model):
    nonzero_indices = nonzero(model.E[:, 0].A)[0]
    Ej = 0
    j = 0
    max_delta = 0
    if len(nonzero_indices) > 1:        for index in nonzero_indices:            if index == i:                continue
            E_temp = calc_E(index, model)
            delta = abs(E_temp - Ei)            if delta > max_delta:
                max_delta = delta
                Ej = E_temp
                j = index    else:
        j = i        while j == i:
            j = int(random.uniform(0, model.m))
        Ej = calc_E(j, model)    return j, Ejdef iterate(i, model):
    yi = model.y[i]
    Ei = calc_E(i, model)
    model.E[i] = [1, Ei]    # 如果alpahi不满足KKT条件, 则进行之后的操作, 选择alphaj, 更新alphai与alphaj, 还有b
    if (yi * Ei > model.toler and model.alphas[i] > 0) or (yi * Ei < -model.toler and model.alphas[i] < model.C):        # alphai不满足KKT条件
        # 选择alphaj
        j, Ej = select_j(Ei, i, model)
        yj = model.y[j] 
        alpha1old = model.alphas[i].copy()
        alpha2old = model.alphas[j].copy()
        eta = model.mapped_data[i, i] + model.mapped_data[j, j] - 2 * model.mapped_data[i, j]   
        if eta <= 0:            return 0
        alpha2new_unclip = alpha2old + yj * (Ei - Ej) / eta        if yi == yj:
            L = max(0, alpha2old + alpha1old - model.C)
            H = min(model.C, alpha1old + alpha2old)        else:
            L = max(0, alpha2old - alpha1old)
            H = min(model.C, model.C - alpha1old + alpha2old)        if L == H:            return 0
        alpha2new = clip_alpha(L, H, alpha2new_unclip)        if abs(alpha2new - alpha2old) < 0.00001:           return 0
        alpha1new = alpha1old + yi * yj * (alpha2old - alpha2new)
        b1new = -Ei - yi * model.mapped_data[i, i] * (alpha1new - alpha1old) \
                - yj * model.mapped_data[j, i] * (alpha2new - alpha2old) + model.b
        b2new = -Ej - yi * model.mapped_data[i, j] * (alpha1new - alpha1old) \
                - yj * model.mapped_data[j, j] * (alpha2new - alpha2old) + model.b
        model.b = calc_b(b1new, b2new)
        model.alphas[i] = alpha1new
        model.alphas[j] = alpha2new
        model.E[i] = [1, calc_E(i, model)]
        model.E[j] = [1, calc_E(j, model)]        return 1
    return 0def smo(X, y, C, toler, iter_num, kernel_param):
    model = Model(X, y.T, C, toler, kernel_param)
    changed_alphas = 0
    current_iter = 0
    for i in range(model.m):
        changed_alphas += iterate(i, model)
        print("iter:%d i:%d,pairs changed %d"
              %(current_iter, i, changed_alphas))
    current_iter += 1
    print('start...') 
    while current_iter < iter_num and changed_alphas > 0:
        changed_alphas = 0
        # 处理支持向量
        alphas_indice = nonzero((model.alphas.A > 0) * (model.alphas.A < C))[0]        for i in alphas_indice:
            changed_alphas += iterate(i, model)
            print("iter:%d i:%d,pairs changed %d"
                  %(current_iter, i, changed_alphas))
        current_iter += 1
    return model.alphas, model.b

• 注意: 在测试SVM的时候, 我们也需要把测试集通过核函数转为m个特征的输入, 所以我们需要将训练集和测试集代入核函数中

作者：Andrew_Chan

原文出处：  https://www.cnblogs.com/megachen/p/10491461.html