SVM之简化版SMO算法

创建一个alpha向量并将其初始化为零向量
当迭代次数小于最大迭代次数时（外循环）
    对数据集中的每个特征向量（内循环）：
        如果该特征向量可以被优化：
            随机选择另外一个特征向量
            同时优化这两个向量
            如果两个向量都不能被优化，退出内循环
    如果所有向量都没有被优化，进行下一次迭代

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
    
def selectJrand(i,m):
    #随机选择一个不等于i的j值
    j=i
    while (j==i):
        j = int(random.uniform(0,m))
    return j
    
def clipAlpha(aj,H,L): #切一切alpha
    #将alpha 限制在 L和 H之间
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj
    
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    #SVM的一种实现，最小序列法（SMO)
    #SMO 简化版
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iter_ = 0
    while (iter_ < maxIter): # 迭代次数小于设置的最大迭代次数
        alphaPairsChanged = 0 #用于记录alpha是否已经进行优化
        for i in range(m): #对数据集遍历
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b #预测的yi值（类别）
            Ei = fXi - float(labelMat[i])#样本i的预测值与真值的误差
            #检测样本是否违反KKT条件
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b #预测的yj值（类别）
                Ej = fXj - float(labelMat[j]) #样本j的预测值与真值的误差
                alphaIold = alphas[i].copy()#记录原始的alphai 的值
                alphaJold = alphas[j].copy()#记录原始的alphaj 的值
                #用于将aj调整到 0到 C之间
                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:
                    continue
                #eta 是 alphaj的最优修改量
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: #作为简化的算法，这里省略了eta=0（很少发生）时求 alphaj的过程。
                    #print ("eta>=0")
                    continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    #print ("j not moving enough")
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#alphai 变化量和 alphaj 变化量 之和为0， 即保持 alphai+alphaj 不变
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[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 changed %d" % (iter_, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):#成功的更新了一对alpha
            iter_ += 1
        else: iter_ = 0
        #print ("iteration number: %d"% iter_)
    return b,alphas
    
def calcWs(alphas,dataArr,classLabels):#计算权重系数
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w
    
def classfy(Xi, w, b): #做分类预测
    y_predict = Xi*w + b
    if y_predict>0:
        return 1
    else:
        return -1
        
        
        
dataArr, labelArr = loadDataSet("testSet.txt")
#print(labelArr)
b, alphas = smoSimple(dataArr, labelArr, C=0.6, toler =0.001, maxIter=40)
w = calcWs(alphas, dataArr, labelArr)
print("偏置b: \n", b,"\n")
#print("alphas:\n", alphas,"\n")
print("权重矩阵w:\n", w,"\n") #预测值Y = Xi *w + b
print("支持向量:")
for i in range(len(alphas)):
    if alphas[i]> 0: #打印支持向量
        print(dataArr[i], labelArr[i])
print("第一个样本预测的分类：")
print(classfy(mat(dataArr)[0], w,b))