Python 3:『機械学習実戦』のサポートベクトルマシン(3)完全版SMO


Python 3:『機械学習実戦』のサポートベクトルマシン(3)完全版SMO
  • 転載作者と出典を明記してください:http://blog.csdn.net/u011475210
  • コードアドレス:https://github.com/WordZzzz/ML/tree/master/Ch06
  • オペレーティングシステム:WINDOWS 10
  • ソフトウェアバージョン:python-3.6.2-amd64
  • 編  者:WordZzzz
  • Python 3マシン学習実戦のサポートベクトルマシン3完全版SMO
  • 前言
  • サポート関数
  • 最適化ルーチン
  • 外部循環コード
  • 分類テスト

  • 前言
    小規模なデータセットでは、前の記事の簡略化版SMOは問題ありませんが、より大きなデータセットでは実行速度が遅くなります.
    完全版SMOと簡略版SMOは,alphaの変更個代数演算の最適化の一環とそっくりである.最適化の過程で、唯一の違いはalphaを選択する方法です.完全版のSMOアルゴリズムはいくつかのスピードアップできる啓発方法を応用した.
    IXプラットSMOアルゴリズムは、1つの外部サイクルによって最初のalphaを選択し、その選択プロセスは、すべてのデータセット上で単一パススキャンを行う2つの方法で切り替えられます.もう1つは、非境界alpha(境界0またはCのalpha値に等しくない)で単一パススキャンを実現することです.データセット全体のスキャンは簡単で、前述したように実現されていますが、非境界alpha値のスキャンを実現する場合は、これらのalpha値のリストを作成してから、このテーブルを遍歴する必要があります.同時に、このステップは既知の変更されないalpha値をスキップします.
      最初のalpha値を選択すると、アルゴリズムは1つの内ループで2番目のalphaを選択します.最適化の過程で、ステップ長を最大化することによって2番目のalpha値が得られます.簡略化版SMOアルゴリズムでは,jを選択した後に誤り率Ejを計算する.しかし、ここでは、誤差値を保存するためのグローバルキャッシュを作成し、ステップ長またはEi-Ejを最大にするalpha値を選択します.
    サポート関数
      は簡略化版と同様に、完全版にもいくつかのサポート関数が必要です.
  • の最も重要なことは、すべての重要な値を保存するためにデータ構造を構築することであり、このプロセスは1つのオブジェクトによって完了することができる.
  • 与えられたalpha値に対して、第1の補助関数calcEk()は、E値を計算して戻ることができる(呼び出しが頻繁であるため、単独で持ち出さなければならない).
  • selectJ()は、2番目のalphaまたはイントラサイクルのalpha値を選択し、最適化のたびに最大ステップ長を採用することを保証するために適切な値を選択するために使用される.
  • updateEk()は、誤差値を計算してキャッシュに格納するために使用されます.
  • '''#######********************************
    Non-Kernel VErsions below
    '''#######********************************
    
    class optStruct:
        """
        Function:            
    
        Input:      dataMatIn:   
                    classLabels:    
                    C:  C
                    toler:   
    
        Output: X:   
                    labelMat:    
                    C:  C
                    tol:   
                    m:     
                    b:   
                    alphas:alphas  
                    eCache:    
        """ 
        def __init__(self, dataMatIn, classLabels, C, toler):
            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)))
    
    def calcEk(oS, k):
        """
        Function:        E
    
        Input:      oS:    
                    k:  
    
        Output: Ek:   E 
        """ 
        #  fXk,        f(x)=w`x + b
        fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k,:].T)) + oS.b   
        #  E 
        Ek = fXk - float(oS.labelMat[k])
        #        E
        return Ek
    
    def selectJ(i, oS, Ei):
        """
        Function:        alpha  
    
        Input:      i:   alpha   
                    oS:    
                    Ei:       alpha    
    
        Output: j:   alpha   
                    Ej:       alpha    
        """ 
        #      
        maxK = -1; maxDeltaE = 0; Ej = 0
        #      
        oS.eCache[i] = [1, Ei]
        #        ,         E    alpha ,   E  
        validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
        #        1,         
        if (len(validEcacheList)) > 1:
            #         
            for k in validEcacheList:
                #      alpha   ,       
                if k == i: continue
                #  k        
                Ek = calcEk(oS, k)
                #   alpha          
                deltaE = abs(Ei - Ek)
                #     
                if (deltaE > maxDeltaE):
                    maxK = k; maxDeltaE = deltaE; Ej = Ek
            #         maxK    Ej
            return maxK, Ej
        #        ,     alpha,      
        else:
            j = selectJrand(i, oS.m)
            Ej = calcEk(oS, j)
        #    j       Ej
        return j, Ej
    
    def updateEk(oS, k):
        """
        Function:         
    
        Input:      oS:    
                    j:alpha   
    
        Output:  
        """ 
        #     k      
        Ek = calcEk(oS, k)
        #       
        oS.eCache[k] = [1, Ek]

    最適化ルーチン
    次に、決定境界を探すための最適化ルーチンについて簡単に説明する.
      大部分のコードは、以前のsmoSimple()と同じであり、違いは以下のとおりである.
  • は、oSで伝達される独自のデータ構造を使用する.
  • は、selectJrand()ではなくselectJ()を使用して2番目のalphaの値を選択する.
  • alpha値の変更時にEcacheを更新します.
  • def innerL(i, oS):
        """
        Function:     SMO        
    
        Input:      oS:    
                    i:alpha   
    
        Output:  
        """ 
        #    
        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)):
            #        alpha 
            j, Ej = selectJ(i, oS, Ei)
            #  copy        ,      
            alphaIold = oS.alphas[i].copy(); alpahJold = oS.alphas[j].copy();
            #  alpha 0 C  
            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
            #         0,             ,     SMO      
            if eta >= 0: print("eta>=0"); return 0
            #    alphas[j]  
            oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
            #   alphas[j]      
            oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
            #      
            updateEk(oS, j)
            #        ,           
            if (abs(oS.alphas[j] - alpahJold) < 0.00001): print("j not moving enough"); return 0
            # i    ,     ,      
            oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alpahJold - oS.alphas[j])
            #      
            updateEk(oS, i)
            #     
            b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :]*oS.X[i, :].T - oS.labelMat[j] * (oS.alphas[j] - alpahJold) * 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] - alpahJold) * oS.X[j, :]*oS.X[j, :].T
            #  0 C  ,    ,       
            if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
            elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[i]): oS.b = b2
            else: oS.b = (b1 + b2) / 2.0
            #    1
            return 1
        #    0
        else: return 0

    がいぶサイクルコード
    外部ループコードの入力は、関数smoSimple()と全く同じです.コード全体のボディはwhileループであり、終了条件:反復回数が指定した最大値を超えた場合、またはセット全体が任意のalphaペアを変更していない場合、ループを終了します.whileループの内部はsmoSimple()とは異なり、最初のforループはデータセット上で任意の可能なalphaを遍歴する.innerl()によって2番目のalphaを選択し、可能な場合に最適化処理を行います.いずれかのペアのalpha値が変化すると、1が返されます.2番目のforループは、すべての非境界alpha値、すなわち境界0またはCにない値を遍歴する.
    def smoP(dataMatIn, classLabels, C, toler, maxIter):
        """
        Function:     SMO  
    
        Input:      dataMatIn:   
                    classLabels:    
                    C:  C
                    toler:   
                    maxIter:       
    
        Output: b:   
                    alphas:    
        """ 
        #        
        oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
        #       
        iter = 0
        #      
        entireSet = True; alphaPairsChanged = 0
        #    :      、         alpha    
        while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
            alphaPairsChanged = 0
            #              
            if entireSet:
                #       alpha 
                for i in range(oS.m):
                    #     alpha ,             
                    alphaPairsChanged += innerL(i, oS)
                    print("fullSet, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
                #      
                iter += 1
            else:
                #        alpha 
                nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
                #        alpha 
                for i in nonBoundIs:
                    #     alpha ,             
                    alphaPairsChanged += innerL(i, oS)
                    print("non-bound, iter: %d i: %d, pairs changed %d" % (iter, i, alphaPairsChanged))
                #      
                iter += 1
            #                 
            if entireSet: entireSet = False
            elif (alphaPairsChanged == 0): entireSet =True
            print("iteration number: %d" % iter)
        #          
        return oS.b, oS.alphas

    テストコードは、興味があれば何度も実行して実行時間の平均値を計算し、簡略化版に比べてどれだけ速くなったかを見ることができます.
    >>> reload(svmMLiA)
    'svmMLiA' from 'E:\\      \\mycode\\Ch06\\svmMLiA.py'>
    >>> dataArr, labelArr = svmMLiA.loadDataSet('testSet.txt')
    >>> b, alphas = svmMLiA.smoP(dataArr, labelArr, 0.6, 0.001, 40)
    L==H
    fullSet, iter: 0 i: 0, pairs changed 0
    L==H
    fullSet, iter: 0 i: 1, pairs changed 0
    fullSet, iter: 0 i: 2, pairs changed 1
    L==H
    ···
    j not moving enough
    fullSet, iter: 2 i: 97, pairs changed 0
    fullSet, iter: 2 i: 98, pairs changed 0
    fullSet, iter: 2 i: 99, pairs changed 0
    iteration number: 3

    以前のようにbとalphaを印刷し、得られたデータを図を描くために使用します.
    >>> b
    matrix([[-2.89901748]])
    >>> alphas[alphas > 0]
    matrix([[ 0.06961952,  0.0169055 ,  0.0169055 ,  0.0272699 ,  0.04522972,
              0.0272699 ,  0.0243898 ,  0.06140181,  0.06140181]])
    >>> from numpy import *
    >>> shape(alphas[alphas > 0])
    (1, 9)
    >>> for i in range(100):
    ...     if alphas[i] > 0.0: print(dataArr[i], labelArr[i])
    ... 
    [3.542485, 1.977398] -1.0
    [2.114999, -0.004466] -1.0
    [8.127113, 1.274372] 1.0
    [4.658191, 3.507396] -1.0
    [8.197181, 1.545132] 1.0
    [7.40786, -0.121961] 1.0
    [6.960661, -0.245353] 1.0
    [6.080573, 0.418886] 1.0
    [3.107511, 0.758367] -1.0

      定数Cは、すべてのサンプルの間隔が1.0以上であることを保証する一方で、分類間隔をできるだけ大きくし、この2つの面でバランスをとる必要がある.Cが大きい場合、分類器は、多くのサンプルに対して超平面を区切ることによって、力図を正しく分類します.この最適化結果は下図のように,支持ベクトルが多くなることが明らかになった.データセットが非線形に分割可能である場合,サポートベクトルは超平面付近で凝集することが分かった.
    Python3:《机器学习实战》之支持向量机(3)完整版SMO_第1张图片
    分類テスト
    さあ、やっと私たちが計算したalpha値を分類することができます.まずalpha値に基づいて超平面に行かなければならない.これにはwの計算も含まれている.
    def calcWs(alphas, dataArr, classLabels):
        """
        Function:     W
    
        Input:      alphas:    
                    dataArr:   
                    classLabels:    
    
        Output: w:w*x+b  w
        """ 
        #     
        X = mat(dataArr); labelMat = mat(classLabels).transpose()
        #       
        m,n = shape(X)
        #   w
        w = zeros((n,1))
        #  alpha,  w
        for i in range(m):
            w += multiply(alphas[i]*labelMat[i],X[i,:].T)
        #  w 
        return w

      上記のコードの中で最も重要なのはforループであり,複数の数の積を実現する.forループはデータセットのすべてのデータを巡回するが,最終的に機能するのはサポートベクトルのみである.
    >>> reload(svmMLiA)
    'svmMLiA' from 'E:\\      \\mycode\\Ch06\\svmMLiA.py'>
    >>> from numpy import *
    >>> ws = svmMLiA.calcWs(alphas, dataArr, labelArr)
    >>> ws
    array([[ 0.65307162],
           [-0.17196128]])
    >>> datMat = mat(dataArr)
    >>> datMat[0]* mat(ws)+b
    matrix([[-0.92555695]])
    >>> labelArr[0]
    -1.0
    >>> datMat[1]* mat(ws)+b
    matrix([[-1.36706674]])
    >>> labelArr[1]
    -1.0
    >>> datMat[2]* mat(ws)+b
    matrix([[ 2.30436336]])
    >>> labelArr[2]
    1.0

      上記のテストでは、計算値が0より大きいと1クラスに属し、0より小さいと-1クラスに属する.
    ここで,線形分類器の紹介は終わり,データセットが非線形に分割可能であれば,コア関数の概念を導入する必要があるので,次の編で紹介する.
    シリーズチュートリアル継続リリース中、購読、注目、コレクション、コメント、いいね~( ̄▽ ̄)~
    完の汪(′▽`▽`)。。。zz