本文采用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