zoukankan html css js c++ java

机器学习6—SVM学习笔记

（一）关于拉格朗日乘子法

首先来了解拉格朗日乘子法，那么为什么需要拉格朗日乘子法？记住，有拉格朗日乘子法的地方，必然是一个组合优化问题。那么带约束的优化问题很好说，就比如说下面这个：

m i n

这是一个带等式约束的优化问题，有目标值，有约束条件。那么想想假设没有约束条件这个问题是怎么求解的呢？是不是直接f对各个x求导等于0,，解x就可以了，可以看到没有约束的话，求导为0，那么各个x均为0吧，这样f=0了，最小。但是x都为0不满足约束条件呀，那么问题就来了。这里在说一点的是，为什么上面说求导为0就可以呢？理论上多数问题是可以的，但是有的问题不可以。如果求导为0一定可以的话，那么f一定是个凸优化问题，什么是凸的呢？像下面这个左图：
这里写图片描述

凸的就是开口朝一个方向（向上或向下）。更准确的数学关系就是：

f ( x 1 ) + f ( x 2 ) 2 > f ( x 1 + x 2 2 ) 或 者

注意的是这个条件是对函数的任意x取值。如果满足第一个就是开口向上的凸，第二个是开口向下的凸。可以看到对于凸问题，你去求导的话，是不是只有一个极点，那么他就是最优点，很合理。类似的看看上图右边这个图，很明显这个条件对任意的x取值不满足，有时满足第一个关系，有时满足第二个关系，对应上面的两处取法就是，所以这种问题就不行，再看看你去对它求导，会得到好几个极点。然而从图上可以看到，只有其中一个极点是最优解，其他的是局部最优解，那么当真实问题的时候你选择那个？说了半天要说啥呢，就是拉格朗日法是一定适合于凸问题的，不一定适合于其他问题，还好我们最终的问题是凸问题。

回头再来看看有约束的问题，既然有了约束不能直接求导，那么如果把约束去掉不就可以了吗？怎么去掉呢？这才需要拉格朗日方法。既然是等式约束，那么我们把这个约束乘一个系数加到目标函数中去，这样就相当于既考虑了原目标函数，也考虑了约束条件，比如上面那个函数，加进去就变为：

m i n

这里可以看到与 $α_{1}, α_{2}$

\partial f \partial x 1 = 4 x 1 + 2 α 1 = 0 \Rightarrow x 1 = - 0.5 α 1

把它在带到约束条件中去，可以看到，2个变量两个等式，可以求解，最终可以得到 $α_{1} = - 0.39, α_{2} = - 1.63$

（二）关于KKT条件

继续讨论关于带等式以及不等式的约束条件的凸函数优化。任何原始问题约束条件无非最多3种，等式约束，大于号约束，小于号约束，而这三种最终通过将约束方程化简化为两类：约束方程等于0和约束方程小于0。再举个简单的方程为例，假设原始约束条件为下列所示：

m i n

$α_{1}, α_{2}$

s . t .

为什么都变成等号与小于号，方便后面的，反正式子的关系没有发生任何变化就行了。

现在将约束拿到目标函数中去就变成：

L (x, α) = f (x) + α 1 g 1 (x) + α 2 g 2 (x)

$α_{1}, α_{2}$

L (x, α, β) = f (x) + \sum α i g i (x) + \sum β i h i (x)

$α_{1}, α_{2}$

(1) L对各个x求导为零；
(2) h(x)=0;
(3) $\sum α_{i} g_{i} (x) = 0 ， α_{i} \geq 0$

这三个式子前两个好理解，重点是第三个式子不好理解，因为我们知道在约束条件变完后，所有的g(x)<=0，且 $α_{i} \geq 0$

假设有一个目标函数，以及它的约束条件，形象的画出来就如下：
这里写图片描述
假设就这么几个吧，最终约束是把自变量约束在一定范围，而函数是在这个范围内寻找最优解。函数开始也不知道该取哪一个值是吧，那就随便取一个，假设某一次取得自变量集合为x1*，发现一看，不满足约束，然后再换呀换，换到了x2*,发现可以了，但是这个时候函数值不是最优的，并且x2*使得g1(x)与g2(x)等于0了，而g3(x)还是小于0。这个时候，我们发现在x2的基础上再寻找一组更优解要靠谁呢？当然是要靠约束条件g1(x)与g2(x)，因为他们等于0了，很极限呀，一不小心，走错了就不满足它们两了，这个时候我们会选择g1(x)与g2(x)的梯度方向往下走，这样才能最大程度的拜托g1(x)与g2(x)=0的命运，使得他们满足小于0的约束条件对不对。至于这个时候需不需要管g3(x)呢？正常来说管不管都可以，如果管了，也取g3在x2*处的梯度的话，因为g3已经满足了小于0的条件，这个时候在取在x2*处的梯度，你能保证它是往好的变了还是往差的变了？答案是都有可能。运气好，往好的变了，可以更快得到结果，运气不好，往差的变了，反而适得其反。那么如果不管呢？因为g1(x)与g2(x)已经在边缘了，所以取它的梯度是一定会让目标函数变好的。综合来看，这个时候我们就不选g3。那么再往下走，假设到了自变量优化到了x3*，这个时候发现g2(x)与g3(x)等于0，也就是走到边了，而g1(x)小于0，可变化的空间绰绰有余，那么这个时候举要取g2(x)与g3(x)的梯度方向作为变化的方向，而不用管g1(x)。那么一直这样走呀走，最终找到最优解。可以看到的是，上述如果g1(x)、g2(x)=0的话，我们是需要优化它的，又因为他们本身的条件是小于0的，所以最终的公式推导上表明，是要乘以一个正系数 $α$

比如上面例子的目标值与约束：

m i n

$α_{1}, α_{2}$

L (x, α) = x 2 1 - 2 x 1 + 1 + x 2 2 + 4 x 2 + 4 +

$α_{1}, α_{2}$

\partial L \partial x 1 = 2 x 1 - 2 - α 1 + 10 α 2 = 0 \Rightarrow x 1 = 0.5 ( α 1 - 10 α 2 + 2

$α_{1}, α_{2}$

α 1 * g 1 (x) = α 1 * (10 - x 1 - 10 x 2) = 0

$α_{1}, α_{2}$

说回来，这里有四种情况，正好两个 $α$

（2） $g_{1} (x) = g_{2} (x) = 0$

（3）其他两种情况再去讨论发现是不行的。

可以看到像这种简单的讨论完以后就可以得到解了。
x1=110/101=1.08;x2=90/101=0.89,那么它得到结果对不对呢？这里因为函数简单，可以在matlab下画出来，同时约束条件也可以画出来，那么原问题以及它的约束面画出来就如下所示：
这里写图片描述
这是截取下来的符合约束要求的目标面

可以看到最优解确实就是上面我们求的那个解。既然简单的问题可以这样解，那么复杂一点的只需要简单化，照样可以解，至此KKT条件解这类约束性问题就是这样，它对后续的SVM求解最优解至关重要。

【机器学习详解】SMO算法剖析

回顾初高中三角函数转换知识：

1、诱导公式：

sin(-α) = -sinα； cos(-α) = cosα； sin(π/2-α) = cosα； cos(π/2-α) = sinα；　　 sin(π/2+α) = cosα； cos(π/2+α) = -sinα；

sin(π-α) = sinα； cos(π-α) = -cosα；　　 sin(π+α) = -sinα； cos(π+α) = -cosα；

tanA= sinA/cosA；

tan(π/2+α)＝－cotα； tan(π/2-α)＝ cotα； tan(π-α)＝－tanα； tan(π+α)＝tanα

2、两角和差公式：

sin(AB) = sinAcosBcosAsinB cos(AB) = cosAcosBsinAsinB tan(AB) = (tanAtanB)/(1tanAtanB) cot(AB) = (cotAcotB1)/(cotBcotA)

3、倍角公式　　

sin2A=2sinA•cosA cos2A=cosA2-sinA2=1-2sinA2=2cosA2-1 tan2A=2tanA/（1-tanA2）=2cotA/(cotA2-1)

4、半角公式　　

tan(A/2)=(1-cosA)/sinA=sinA/(1+cosA); cot(A/2)=sinA/(1-cosA)=(1+cosA)/sinA.

sin^2(a/2)=(1-cos(a))/2 cos^2(a/2)=(1+cos(a))/2

tan(a/2)=(1-cos(a))/sin(a)=sin(a)/(1+cos(a))
5、和差化积　

sinθ+sinφ = 2 sin[(θ+φ)/2] cos[(θ-φ)/2] sinθ-sinφ = 2 cos[(θ+φ)/2]sin[(θ-φ)/2]

cosθ+cosφ = 2 cos[(θ+φ)/2]cos[(θ-φ)/2] cosθ-cosφ = -2 sin[(θ+φ)/2]sin[(θ-φ)/2]

tanA+tanB=sin(A+B)/cosAcosB=tan(A+B)(1-tanAtanB) tanA-tanB=sin(A-B)/cosAcosB=tan(A-B)(1+tanAtanB)

6、积化和差　　

sinαsinβ = -1/2*[cos(α-β)-cos(α+β)] cosαcosβ = 1/2*[cos(α+β)+cos(α-β)]

sinαcosβ = 1/2*[sin(α+β)+sin(α-β)] cosαsinβ = 1/2*[sin(α+β)-sin(α-β)]万能公式

两条直线相互垂直的条件:

两条直线在同一平面内
1、如果斜率为k1和k2，那么这两条直线垂直的充要条件是k1·k2=－1 ， tanαtan(α+π/2) = tanα(－cotα) = -1

2、如果一直线不存在斜率，则两直线垂直时，一直线的斜率必然为零.

3、两直线(A₁x+B₁y+C₁=0, A₂x+B₂y+C₂=0)垂直的充要条件是：A1A2+B1B2=0.

向量1 (x1,y1),长度 L1 =√(x1²+y1²)
向量2 (x2,y2),长度 L2 =√(x2²+y2²)
（x1,y1）到(x2,y2)的距离：D=√[(x1 - x2)² + (y1 - y2)²]
两个向量垂直,根据勾股定理：L1² + L2² = D²
∴ (x1²+y1²) + (x2²+y2²) = (x1 - x2)² + (y1 - y2)²
∴ x1² + y1² + x2² + y2² = x1² -2x1x2 + x2² + y1² - 2y1y2 + y2²
∴ 0 = -2x1x2 - 2y1y2
∴ x1x2 + y1y2 = 0
该定理还可以扩展到三维向量：x1x2 + y1y2 + z1z2 = 0,那么向量(x1,y1,z1)和（x2,y2,z2）垂直
甚至扩展到更高维度的向量,两个向量L1,L2垂直的充分必要条件是:L1×L2=0

test6.py

#-*- coding：utf-8

import sys
sys.path.append("svmMLiA.py")

import svmMLiA
from numpy import *

dataArr, labelArr = svmMLiA.loadDataSet("testSet.txt")
# b, alphas = svmMLiA.smoSimple(dataArr, labelArr, 0.6, 0.001, 40)

b, alphas = svmMLiA.smoP(dataArr, labelArr, 0.6, 0.001, 40)

ws = svmMLiA.calcWs(alphas, dataArr, labelArr)

num, n = shape(alphas[alphas > 0])

print("num = %d, n = %d" % (num, n))

for i in list(range(100)):
    if alphas[i] > 0.0:
        print(dataArr[i], labelArr[i])

print("ws : ")
print(ws)

datMat = mat(dataArr)
caculLabel = datMat[0]*mat(ws) + b
print("caculLabel : ")
print(caculLabel)

print(labelArr[0])


svmMLiA.testRbf()

svmMLiA.testDigits(("rbf", 20))

print("over!")

svmMLiA.py

'''
Created on Nov 4, 2010
Chapter 5 source file for Machine Learing in Action
@author: Peter
'''
import matplotlib.pyplot as plt
from numpy import *
from time import sleep

def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split('	')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat

def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    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
        for i in range(m):
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[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)):
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[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
                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
                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])#update i by the same amount as j
                                                                        #the update is in the oppostie direction
                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): iter += 1
        else: iter = 0
        print("iteration number: %d" % iter)
    return b,alphas

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
        K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
    else: raise NameError('Houston We Have a Problem -- 
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
        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
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

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

def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            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

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

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)):
        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
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        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

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        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]
            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)
    return oS.b,oS.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 testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelMat[svInd];
    print("there are %d Support Vectors" % shape(sVs)[0])
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m, n = shape(datMat)

    # 添加绘图部分
    xcord0 = []
    ycord0 = []
    xcord1 = []
    ycord1 = []

    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
        if sign(predict) == 1:
            xcord0.extend(datMat[i, :].tolist())
        else:
            xcord1.extend(datMat[i, :].tolist())
    print("the test error rate is: %f" % (float(errorCount) / m))

    fig = plt.figure()
    ax = fig.add_subplot(111)
    s0 = shape(xcord0)[0]
    for i in range(s0):
        ax.scatter(xcord0[i][0], xcord0[i][1], marker='o', s=20, c='blue')
    s1 = shape(xcord1)[0]
    for i in range(s1):
        ax.scatter(xcord1[i][0], xcord1[i][1], marker='v', s=35, c='red')
    plt.show()
    # 上绘图部分


def img2vector(filename):
    returnVect = zeros((1,1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0,32*i+j] = int(lineStr[j])
    return returnVect

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)           #load the training set
    m = len(trainingFileList)
    trainingMat = zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     #take off .txt
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels

def testDigits(kTup=('rbf', 10)):
    dataArr,labelArr = loadImages('trainingDigits')
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd]
    labelSV = labelMat[svInd];
    print("there are %d Support Vectors" % shape(sVs)[0])
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr,labelArr = loadImages('testDigits')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount)/m))


'''#######********************************
Non-Kernel VErsions below
'''#######********************************

class optStructK:
    def __init__(self,dataMatIn, classLabels, C, toler):  # Initialize the structure with the parameters
        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

def calcEkK(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJK(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            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

def updateEkK(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]

def innerLK(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)):
        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.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T
        if eta >= 0: print("eta>=0"); return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        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

def smoPK(dataMatIn, classLabels, C, toler, maxIter):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        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]
            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)
    return oS.b,oS.alphas

查看全文

相关阅读:
UVa 11988
UVa 442
.MySQL数据库技术
 Mysql数据库技术
 JDBC技术
 JDBC技术
 JavaSE编程基础
 JavaSE编程基础
 JavaSE编程基础
 web安全性测试

原文地址：https://www.cnblogs.com/Vae1990Silence/p/8377893.html

热门文章
A1021 Deepest Root (25分)
A10131013 Battle Over Cities (25分)
UVa 101
UVa 12169
UVa 11582
UVa 1590
UVa 712
UVa 679
UVa 221
UVa 814

机器学习6—SVM学习笔记

机器学习牛人博客

机器学习实战之SVM

三种SVM的对偶问题

拉格朗日乘子法和KKT条件

支持向量机通俗导论（理解SVM的三层境界）

解密SVM系列（一）：关于拉格朗日乘子法和KKT条件

解密SVM系列（二）：SVM的理论基础

解密SVM系列（三）：SMO算法原理与实战求解

（一）关于拉格朗日乘子法

（二）关于KKT条件

【机器学习详解】SMO算法剖析

两条直线相互垂直的条件: