機器學習實戰6（SMO演算法）

georgefuxiaosong發表於2018-03-20

關於SVM的講解可以參閱《機器學習技法》的相關筆記：http://blog.csdn.net/u011239443/article/details/76572743

SMO演算法

實現

# -*- coding: utf-8 -*-
from numpy import *
# 載入資料
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

# 從0到m中產生一個不為i的整數
def selectJrand(i,m):
    j = i
    while(j == i):
        j = int(random.uniform(0,m))
    return j

# 使得aj 在邊界值[L,H]以內
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


#SMO 序列最小優化
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
    
    dataMat = mat(dataMatIn);labelMat = mat(classLabels).transpose();
    b = 0;m,n = shape(dataMat)
    alphas = mat(zeros((m,1)))
    iter = 0
    # alphaPairsChanged 用來更新的次數
    # 當遍歷  連續無更新 maxIter 輪，則認為收斂，迭代結束 
    while iter < maxIter:
        alphaPairsChanged = 0
        for i in range(m):
            # KKT 條件計算出
            fXi = float(multiply(alphas,labelMat).T * (dataMat*dataMat[i,:].T))+b
            # 誤差
            Ei = fXi - float(labelMat[i])
            # toler：容忍錯誤的程度
            # labelMat[i]*Ei < -toler 則需要alphas[i]增大，但是不能>=C
            # labelMat[i]*Ei > toler 則需要alphas[i]減小，但是不能<=0
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                # 從0到m中產生一個不為i的整數
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T * (dataMat*dataMat[j,:].T)+b)
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
               # 演算法講解“二變數優化問題”部分
                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

               # 見式 7.107
                eta = -2.0* dataMat[i,:] * dataMat[j,:].T + dataMat[i,:] *\
                                                            dataMat[i,:].T + dataMat[j,:] * dataMat[j,:].T
                if eta <= 0:print "eta <= 0";continue
                # 見式 7.106
                alphas[j] += labelMat[j]*(Ei-Ej)/eta
                # 見式 7.108 ，講alphas[j] 約束在 [L,H]
                alphas[j] = clipAlpha(alphas[j],H,L)
                if(abs(alphas[j] - alphaJold) < 0.00001):
                    print "j not moving"
               # 見式 7.109
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
               # 見式 7.114
                b1 = b - Ei -labelMat[i]*(alphas[i]-alphaIold)*dataMat[i,:]*dataMat[i,:].T\
                     - labelMat[j]*(alphas[j]-alphaJold)*dataMat[i,:]*dataMat[j,:].T
                b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMat[i, :] * dataMat[j, :].T \
                     - labelMat[j] * (alphas[j] - alphaJold) * dataMat[j, :] * dataMat[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: %d" %(iter,i,alphaPairsChanged)

        if alphaPairsChanged == 0:
            iter += 1
        else:
            iter = 0
        print "iter %d" %iter

    return b,alphas