机器学习实战SMO算法源码解析

from numpy import *
from time import sleep

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

# <editor-fold desc="随机取不等于i的整数j">
def selectJrand(i,m): # 在0至m内随便取一个不等于i的整数, i是第一个alpha的下标, m则是所有alpha的数目
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j
# </editor-fold>


# <editor-fold desc="门限函数">
def clipAlpha(aj,H,L): # a2太大或太小时要进行调整
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj
# </editor-fold>

# <editor-fold desc="简化版SMO">
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    '''

    :param dataMatIn: 数据集
    :param classLabels: 类别标签
    :param C: 参数C
    :param toler: 容错率（松弛变量）
    :param maxIter:取消前最大循环次数
    :return:返回b 和alpha
    '''
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose() # 转换成矩阵，并转置，标记成为一个列向量，每一行和数据矩阵对应
    b = 0
    m,n = shape(dataMatrix)  # 行，列
    alphas = mat(zeros((m,1))) # 创建alpha向量，初始化为0向量
    iter = 0
    while (iter < maxIter): # 当迭代次数小于最大迭代次数
        alphaPairsChanged = 0 # 标记位，记录alpha在该次循环中，有没有完成优化，若完成优化，则迭代次数+1，否则不加1
        for i in range(m): # 遍历样本
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b # 样本i的预测值,multiply:对应元素相乘
            Ei = fXi - float(labelMat[i])  # 样本i的误差
            #if checks if an example violates KKT conditions
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): # 第一次筛选alpha_i，依据:原论文判断条件
                '''
                labelMat[i] = yi
                Ei = g(xi) - yi
                g(xi) = w*xi + b
                所以：labelMat[i]*Ei = yi*(g(xi) - yi) = yi*g(xi) - 1
                yi*g(xi) = labelMat[i]*Ei + 1
                0 <= alpha_i <= C
                当alpha_i = 0或C时，后面调整时将受门限函数制约，调整为L或H
                加入罚参数C的KKT条件为 alpha_i*(yi*g(xi) - 1) = 0
                所以选KKT条件 ：0 < alpha_i < C，等价于 yi*g(xi) = 1 ，等价于labelMat[i]*Ei + 1 = 1 即 labelMat[i]*Ei = 0

                **KKT条件比较苛刻，需要一个容忍值toler， 那么上式 labelMat[i]*Ei = 0 的条件放宽为 abs(labelMat[i]*Ei) < toler

                所以KKT条件为： abs(labelMat[i]*Ei) < toler 且 0 < alpha_i < C
                至于if conditions 是否等价于 不满足KKT的条件， 目前没想明白
                原论文给出的判断条件 ((r2 < -tol && alph2 < C) || (r2 > tol && alph2 > 0))
                判断条件可以结合原问题的拉格朗日函数来看L(w,b,alpha) = 1/2 * ||w||**2 -∑ alpha * (yi*g(xi)-1)
                '''
                j = selectJrand(i,m) # 随机选择第二个alpha
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b # 样本j的预测值
                Ej = fXj - float(labelMat[j]) # 样本j的误差

                # alpha_j_old = alphaJold
                alphaIold = alphas[i].copy() # 浅拷贝，分配新的内存，二者id不同，若用赋值的话id相同，a.append的话b也会变
                alphaJold = alphas[j].copy() # 对于非容器类型（如数字、字符串、和其他'原子'类型的对象）没有拷贝这一说

                # 书上p126页关于alpha_2_new的取值范围，分为y1==y2与y1!=y2两种情况
                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

                # K11+K22-2K12，K11 = X1*X1.T，K12 = X1*X2.T, 这里是原结论的负数形式，那么理论上eta应该小于0
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0:
                    print("eta>=0")
                    continue

                alphas[j] -= labelMat[j]*(Ei - Ej)/eta # alpha_2 未经剪辑时的解
                alphas[j] = clipAlpha(alphas[j],H,L)  # 门限函数阻止alpha_2的修改量过大(alpha_2剪辑后的解)
                # 这里不更新误差

                # 新加的一步优化：若更新步长太小，则放弃
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    print("j not moving enough")
                    continue

                # 书上原公式：alpha_1_new = alpha_1_old + y1*y2(alpha_2_old - alpha_2_new),没毛病
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])

                # b1_new = b1_old - E1 - y1* K11 *(alpha_1_new - alpha_1_old) - y2* K21 *(alpha_2_new - alpha_2_old)
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                # b2_new = b2_old - E2 - y1* K12 *(alpha_1_new - alpha_1_old) - y2* K22 *(alpha_1_new - alpha_1_old)
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T

                # 如果 alpha_1_new，alpha_2_new 同时在(0,C),则 b1_new = b2_new，此时直接赋值,否则取中点
                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))

        # 如果迭代完成，迭代次数+1
        if (alphaPairsChanged == 0):
            iter += 1
        else:
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas
# </editor-fold>

# <editor-fold desc="核转换函数，实际上就是计算K11,K22">
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
    m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin':
        K = X * A.T   #linear kernel
    elif kTup[0]=='rbf': # 高斯核
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T # ||x-y||**2
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab,kTup[1]为高斯函数方差
    else:
        raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K
# </editor-fold>

# <editor-fold desc="新建一个数据结构">
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        '''

        :param dataMatIn: 数据集
        :param classLabels: 类别标签
        :param C: 参数C
        :param toler: 容错率(松弛变量)
        :param kTup: 核函数类型及其参数
        '''
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag,误差E缓存,第一列是是否有效的标志位，第二列才是实际值
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
# </editor-fold>

# <editor-fold desc="计算误差E，k为下标">
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
# </editor-fold>

# <editor-fold desc="第2个变量的选择">
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1
    maxDeltaE = 0 #储存最大的|Ei-Ej|
    Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E，将第一列置1，表示有效，第二列存储误差值Ei
    validEcacheList = nonzero(oS.eCache[:,0].A)[0] # 返回第一列非0的坐标值，即Ei的i
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E 循环找出最大的|Ei-Ej|
            if k == i:
                continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej
# </editor-fold>

# <editor-fold desc="更新误差缓存">
def updateEk(oS, k):#after any alpha has changed update the new value in the cache 当alpha_k的值改变，对应的Ek误差值也变
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
# </editor-fold>

# <editor-fold desc="两个变量alpha_i和alpha_j更新">
def innerL(i, oS):
    Ei = calcEk(oS, i)
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):# KKT条件
        j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H:
            print("L==H")
            return 0
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0:
            print("eta>=0")
            return 0
        # 先更新第二个变量alpha_j
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        # 更新完alpha_j，更新误差缓存中的Ej，oS.b没有更新，如何更新的E？
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("j not moving enough")
            return 0

        # 更新第一个变量alpha_i
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        # 更新完alpha_i，更新误差缓存中的Ei，oS.b没有更新，如何更新的E？
        updateEk(oS, i) #added this for the Ecache   #the update is in the oppostie direction

        # 这里的Ei和Ej是alpha更新前的误差
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]

        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
            oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0
# </editor-fold>

# <editor-fold desc="完全体SMO">
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup) # 把数据丢进oS数据结构里
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        # 书中外层循环先遍历所有满足“alpha大于零小于C的样本”， 即在间隔边界上的支持向量点
        # 本程序外层循环是支持向量点和所有点交替选择
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):# 遍历样本，第一个变量选择
                alphaPairsChanged += innerL(i,oS) # 变量更新次数
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] # alpha小于零大于C的样本下标
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS) # 变量更新次数
                print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:
            entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
        # 循环终止条件：遍历所有样本点，alpha无更新，或者 迭代次数达设定值
    return oS.b,oS.alphas
# </editor-fold>

# <editor-fold desc="计算参数w">
def calcWs(alphas,dataArr,classLabels): # 计算 w*x + b 中的 w
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m,n = shape(X) # m个样本，n个特征
    w = zeros((n,1)) # W为n行1列，X*W = L.T
    for i in range(m):
        # w = ∑ alpha_i * xi * yi
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w
# </editor-fold>