SVM 采用smo算法计算 python

本文采用smo算法计算svm

程序有点问题，开始才用的libsvm的代码，准备将其java代码写成python的，后面发现用libsvm的数据格式老是出问题。就参考了

机器学习实战的代码。

程序有很多要优化的地方

1）核函数要完善，这里只写了线性核函数。但是整个程序中没有用核函数进行计算。

2）一些异常状况的处理。

整个迭代公式可以参考

http://blog.csdn.net/macyang/article/details/38782399/

个人觉得非常棒，就是后面的smo要各种计算，推导。

其实最后迭代也是比较简单的：

1）找出误差ei

2）找出2k(i,j)-k(i,i)-k(j,j)

3)在ai的约束条件下，更新ai

4）更新aj

5）更新b

6）返回结果

#coding=utf-8
#!/usr/bin/python
import pprint
from numpy import *
import matplotlib.pyplot as plt

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 load_data(path):
    '''
    :param path:这里将libsvm中的数据格式进行读取
    :return:
    '''
    data_set=[]
    label_set=[]
    file_object=open(path)
    for line in file_object.readlines():
        lineArr = line.strip().split(" ")
        label_set.append(int(lineArr[0])) #第一列作为lable
        tempdataline=[]
        for i in range(1,len(lineArr)):
            tempdataline.append(float(lineArr[i].split(':')[1]))
        data_set.append(tempdataline)
    print shape(mat(data_set))
    return data_set,label_set
class mySVM(object):
    def __init__(self):
        '''
        :这里先初始化数据和权重
        '''
        self.data_set=[]
        self.label_set=[]
        self.b=[]           #最后结果中的位移偏量
        self.a=[]           #最后结果中的拉格朗日乘子
        self.kernel=[]      #核函数
        self.c=0.6 #设置惩罚因子
        self.tol=0.01 #设置容忍极限
    def k(self,x1,x2):
        '''
        :param x1: 传入两个向量
        :param x2:
        :return: 返回核函数
        '''
        linek=mat(x1).transpose()*mat(x2) #这里先采用线性核函数,后面可以改成高斯核函数
        return linek[0][0] #因为矩阵乘法最后返回的是矩阵,所以要取第一个元素

    def train(self,data_set,label_set):
        '''
        :param data_set:传入数据集
        :param label_set: 传入标签集
        :return: 返回训练好的a和b
        '''
        data_set=mat(data_set)
        label_set=mat(label_set).transpose()
        m,n = shape(data_set)
        pprint.pprint(data_set)
        print(m,n)
        a=mat(zeros((m,1)))
        max_iter=100 #迭代次数
        b=0
        toler=0.1
        c=1
        print a
        print lable_set
        while(max_iter>=0):
            max_iter=max_iter-1
            for i in range(m):#遍历a
                #计算  ei=ui-yi
                tempu=multiply(a,label_set).T*data_set*data_set[i,:].T
                ui=tempu[0][0] + b
                ei=ui-label_set[i][0]
                '''
                 * 把违背KKT条件的ai作为第一个
                 * 满足KKT条件的情况是：
                 * yi*f(i) >= 1 and alpha == 0 (正确分类)
                 * yi*f(i) == 1 and 0<alpha < C (在边界上的支持向量)
                 * yi*f(i) <= 1 and alpha == C (在边界之间)
                 *
                 *  这里借用了另一个网友的话
                 *
                 * ri = y[i] * Ei = y[i] * f(i) - y[i]^2 >= 0
                 * 如果ri < 0并且alpha < C 则违反了KKT条件
                 * 因为原本ri < 0 应该对应的是alpha = C
                 * 同理，ri > 0并且alpha > 0则违反了KKT条件
                 * 因为原本ri > 0对应的应该是alpha =0
                '''
                if((label_set[i]*ei<-toler) and (a[i]<c)) or ((label_set[i]*ei>toler) and (a[i]>0)):
                    j=0 #这里应该是寻找MAX|E1 - E2|   ,为了简单点,就随机选择了一个.只是速度要慢些
                    while(j==i):
                        j=int(random.uniform(0, m))


                    uj=float(multiply(a, label_set).T*data_set*data_set[j, :].T)+b
                    ej=uj-float(label_set[j])
                    aj_old=a[j]
                    ai_old=a[i]
                    if(label_set[i]<>label_set[j]):  #开始计算L和H
                        L=max(0, label_set[j]-a[i])
                        H=min(c, c+a[j]-a[i])
                    else:
                        L=max(0, a[j]+a[i]-c)
                        H=min(c, a[j]+a[i])

                    #这里计算2k(i,j)-k(i,i)-k(j,j) ,这里应该用核函数,先这样将就用了.不先核化
                    eta=2.0*data_set[i, :]*data_set[j, :].T-data_set[i, :]*data_set[i, :].T- data_set[j, :]*data_set[j, :].T

                    if (eta >= 0):#如果eta等于0或者大于0 则表明a最优值应该在L或者U上
                        continue

                    '''
                    这里更新a[j]
                    '''
                    a[j]=a[j]-label_set[j]*(ei-ej)/eta
                    if(a[j]<L):
                        a[j]=L
                    elif(a[j]>H):
                        a[j]=H
                    #更新ai
                    a[i]+=label_set[j]*label_set[i]*(aj_old-a[j])
                    #更新b
                    '''
                    根据公式
                    bnew1 =bold −E1 −y1 (α1new −α1old)K(x1 ,x1 )−y2 (α2new −α2old)K(x2 ,x2 )
                    bnew2 =bold −E2 −y1 (α1new −α1old)K(x1 ,x1 )−y2 (α2new −α2old)K(x2 ,x2 )
                    '''
                    b1=b-ei-label_set[i]*(a[i]-ai_old)*data_set[i, :]*data_set[i, :].T- \
                       label_set[j]*(a[j]-aj_old)* data_set[i, :]*data_set[j, :].T
                    b2=b-ej-label_set[i]*(a[i]-ai_old)*data_set[i, :]*data_set[i, :].T- \
                       label_set[j]*(a[j]-aj_old)* data_set[i, :]*data_set[j, :].T
                    #这里选择合适的b

                    if(0<a[i]) and (c>a[i]):
                        b=b1
                    elif(0<a[j]) and (c>a[j]):
                        b=b2
                    else:
                        b=(b1+b2)/2.0
        return a,b

if __name__ == '__main__':
    print("-------start load data-----")
    c=0.6 #定义系数c
    path="./testSet.txt"
    data_set,lable_set=loadDataSet(fileName=path)
    mysvm= mySVM()
    a,b=mysvm.train(data_set=data_set,label_set=lable_set)
    print a,b