機械学習のサポートベクトルマシンpython実装

23650 ワード

前言:紙の上で得たのは結局浅く感じて、絶対にこの事がお辞儀をすることを知っています.
  • 理論の基礎
  • 関数間隔と幾何学的間隔
  • 最適化目標
  • 導出プロセス
  • サポートベクトル
  • SMOシーケンス最小最適化アルゴリズム
  • 核テクニック

  • 二python実装
  • から

    一.りろんきそ
    1.関数間隔とジオメトリ間隔
    関数間隔:γ^=yi(w∗xi+b)ジオメトリ間隔:γ=yi(w∗xi+b∣∣w∣∣)
    2.最適化目標:
    minα      12∑i=0m∑j=0mαiαjyiyjxixj−∑i=0mαi
    s.t.   ∑i=0mαiyi=0
    s.t.   0≤αi≤C,   (i=1,2,⋯,m)
    3.導出プロセス:
    我々の目標は境界に最も近い点の幾何学的間隔を最大化し,形式化できることである.
    maxw,bγ
    s.t.        yi(w∗xi+b∣∣w∣∣)≥γ   ,   (i=1,2,…,m)
    ジオメトリ間隔と関数間隔の関係に基づいて
    maxw,bγ^∣∣w∣∣
    s.t.        yi(w∗xi+b)≥γ^   ,   (i=1,2,…,m)
    によってγ^ w,bと同期して変化するのでγ のスケールは上記の制約条件に影響しません.γ^=1
    maxw,b1∣∣w∣∣
    s.t.        yi(w∗xi+b)≥1   ,   (i=1,2,…,m)
    最大化1
    minw,b12∣∣w∣∣2
    s.t.        yi(w∗xi+b)≥1   ,   (i=1,2,…,m)
    ソフト間隔を加えると、トレーニングデータの線形不可分状況に適しています.
    minw,b,ξ12∣∣w∣∣2+C∑i=1mξi
    s.t.        yi(w∗xi+b)≥1−ξi   ,   (i=1,2,…,m)
    ξi≥0   ,   (i=1,2,…,m)
    これは,ラグランジュ対偶性(参照)を適用し,対偶問題を解くことによって元の問題の最適解を得る凸最適化問題である.このような利点:*対偶問題は往々にして解きやすい;対偶形式は自然に核関数を導入することができ、さらに非線形分類問題に広がる.
    まずラグランジュ関数を作成します.
    L(w,b,ξ,α,β)=12∣∣w∣∣2+C∑i=1mξi+∑i=0mαi[1−ξi−yi(w∗xi+b)]−∑i=0mβξi
    ラグランジュの対偶性によれば、元の問題の対偶問題は極めて小さな問題である.
    maxα,βminw,b,ξL(w,b,ξ,α,β)
    1.求める
    minw,b,ξL(w,b,ξ,α,β)
    はい
    w,b,ξ バイアス導関数を求めて0に等しくする
    ∇wL(w,b,ξ,α,β)=w−∑i=0mαiyixi=0
    ∇bL(w,b,ξ,α,β)=−∑i=0mαiyi=0
    ∇ξiL(w,b,ξ,α,β)=C−αi−βi=0
    得る
    w=∑i=0mαiyixi
    ∑i=0mαiyi=0
    C=αi+βi
    ラグランジュ関数の代入
    L(w,b,ξ,α,β)=−12∑i=0m∑j=0mαiαjyiyjxixj+∑i=0mαi
  • maxを求めるαL(w,b,ξ,α,β)

  • maxα      −12∑i=0m∑j=0mαiαjyiyjxixj+∑i=0mαi
    s.t.   ∑i=0mαiyi=0
    s.t.   αi≥0,   (i=1,2,⋯,m)
    に等しい
    ∈α      12∑i=0m∑j=0mαiαjyiyjxixj−∑i=0mαi
    s.t.   ∑i=0mαiyi=0
    s.t.   αi≥0,   (i=1,2,⋯,m)
    4.解を求める
    上記最適化問題の解を求め,KKT条件により元の問題の解を得ることができる.セットα∗は対偶問題の解で、KKTの条件によって
    ∇wL(w∗,b∗,ξ∗,α∗,β∗)=w∗−∑i=0mα∗iyixi=0
    ∇bL(w∗,b∗,ξ∗,α∗,β∗)=−∑i=0mα∗iyi=0
    ∇ξiL(w∗,b∗,ξ∗,α∗,β∗)=C−α∗i−β∗i=0
    α∗i(1−ξ∗i−yi(w∗∗xi+b∗))=0
    β∗iξ∗i=0
    ξ∗i≥0
    α∗i≥0
    β∗i≥0
    少なくとも1つは存在する
    α∗j>0(反証法、もし
    α∗=0であれば得られる
    w∗=0であり、
    w∗=0は元の最適化問題の解ではなく矛盾を生じる),これに対して
    j上記式から導出
    yj(w∗∗xj+b∗)−1=0
    なぜなら
    y 2 j=1、上式の両側を同乗する
    yj,則得
    b∗=yj−∑i=0mα∗iyi(xixj)
    同時に、
    w∗=∑i=0mα∗iy∗ix∗i
    すなわち,原問題の最適化解が得られる.
    5.サポートベクトル
    線形不可分の場合,対偶問題の解をα∗にはαj>0のサンプル点の例xiをサポートベクトルと呼ぶ、ソフト間隔のサポートベクトルは、間隔境界上、または間隔境界と分離超平面との間、または分離超平面の誤分側であってもよい.KKTから得られた公式は、0<αi1である、支持ベクトルは分離超平面の誤分側に位置する.
    6.SMO(シーケンス最小最適化アルゴリズム)
    最適化の目標は次のとおりです.
    minα      12∑i=0m∑j=0mαiαjyiyjxixj−∑i=0mαi
    s.t.   ∑i=0mαiyi=0
                  0≤αi≤C,   (i=1,2,⋯,m)
    基本構想:すべての変数の解がこの最適化問題のKKT条件を満たすならば、この最適化問題の解は得られ、KKT条件はこの最適化問題の十分な必要条件であるからである.そうでなければ、2つの変数を選択して、他の変数を固定し、この2つの変数に対して2次計画問題を構築します.この2つの変数に関する2次計画問題の解は、元の2次計画問題の目標関数値をより小さくするため、元の2次計画問題の解に近いはずです.重要なのは,このときサブ問題を解析法により解くことができ,アルゴリズム全体の計算速度を大幅に向上させることができることである.サブ問題には2つの変数があり,1つはKKT条件に違反する最も深刻なものであり,もう1つは制約条件によって自動的に決定される.このように,SMOアルゴリズムは原問題を絶えずサブ問題に分解し,サブ問題を解き,さらに原問題を解く目的を達成する.(参考)
    アルゴリズムステップ:1.初値をとるα(0)=0、k=0、kを反復回数とする.2.最適化変数の選択αk1,αk 2,2つの変数の最適化問題を解析して解き,最適解を求めるαk+11,αk+12,更新α を選択します.αk+1 ; 3.精度でϵ 範囲内で停止条件を満たす
    ∑i=0mαiyi=0
    0≤αi≤C
    yig(xi)=⎧⎩⎨≥1,=1,≤1,{xi|αi=0}{xi|0≤αi≤C}{xi|αi=C}
    ただし、g(xi)=Σmj=0αjyjK(xj,xi)+bは4ステップ目に移行する.そうでなければk=k+1を2ステップ目に回します.4.取るα^=α(k+1) .
    7.核技術
    入力空間における非線形分類問題については,非線形変換によりある高次元特徴空間における線形分類問題に変換し,高次元特徴空間において線形サポートベクトルマシンを学習することができる.線形サポートベクトルマシン学習の対偶問題では,ターゲット関数と分類決定関数はいずれもインスタンスとインスタンス間の内積にのみ関与するため,非線形変換を表示的に指定するのではなく,中の内積をコア関数で置き換える.コア関数は,1つの非線形変換後の2つのインスタンス間の内積を表す.具体的には、K(x,z)は、入力空間から1つのコア関数、または正定コアであり、χ 特徴空間Hへのマッピングϕ(x):χ→H,任意x,z∈χ ,あります
    K(x,z)=ϕ(x)∗ϕ(z)
    たいしょうかんすう
    K(x,z)が正定核である要件は以下の通りである.
    xi∈χ ,任意の正の整数m、対称関数
    K(x,z)に対応するGram行列は半正定である.
    したがって,ベクトルマシン学習を線形にサポートする対偶問題では,コア関数を用いる.
    K(x,z)は内積の代わりに,非線形支持ベクトルマシンを解く.
    f(x)=sign(∑i=0mα∗iyiK(x,z)+b∗)
    二.pythonインプリメンテーション(引用)
    from numpy import *
    import time
    import matplotlib.pyplot as plt
    
    
    # calulate kernel value
    def calcKernelValue(matrix_x, sample_x, kernelOption):
        kernelType = kernelOption[0]
        numSamples = matrix_x.shape[0]
        kernelValue = mat(zeros((numSamples, 1)))
    
        if kernelType == 'linear':
            kernelValue = matrix_x * sample_x.T
        elif kernelType == 'rbf':
            sigma = kernelOption[1]
            if sigma == 0:
                sigma = 1.0
            for i in xrange(numSamples):
                diff = matrix_x[i, :] - sample_x
                kernelValue[i] = exp(diff * diff.T / (-2.0 * sigma**2))
        else:
            raise NameError('Not support kernel type! You can use linear or rbf!')
        return kernelValue
    
    
    # calculate kernel matrix given train set and kernel type
    def calcKernelMatrix(train_x, kernelOption):
        numSamples = train_x.shape[0]
        kernelMatrix = mat(zeros((numSamples, numSamples)))
        for i in xrange(numSamples):
            kernelMatrix[:, i] = calcKernelValue(
                train_x, train_x[i, :], kernelOption)
        return kernelMatrix
    
    
    # define a struct just for storing variables and data
    class SVMStruct:
        def __init__(self, dataSet, labels, C, toler, kernelOption):
            self.train_x = dataSet  # each row stands for a sample
            self.train_y = labels  # corresponding label
            self.C = C             # slack variable
            self.toler = toler     # termination condition for iteration
            self.numSamples = dataSet.shape[0]  # number of samples
            # Lagrange factors for all samples
            self.alphas = mat(zeros((self.numSamples, 1)))
            self.b = 0
            self.errorCache = mat(zeros((self.numSamples, 2)))
            self.kernelOpt = kernelOption
            self.kernelMat = calcKernelMatrix(self.train_x, self.kernelOpt)
    
    
    # calculate the error for alpha k
    def calcError(svm, alpha_k):
        output_k = float(multiply(svm.alphas, svm.train_y).T *
                         svm.kernelMat[:, alpha_k] + svm.b)
        error_k = output_k - float(svm.train_y[alpha_k])
        return error_k
    
    
    # update the error cache for alpha k after optimize alpha k
    def updateError(svm, alpha_k):
        error = calcError(svm, alpha_k)
        svm.errorCache[alpha_k] = [1, error]
    
    
    # select alpha j which has the biggest step
    def selectAlpha_j(svm, alpha_i, error_i):
        svm.errorCache[alpha_i] = [1, error_i]  # mark as valid(has been optimized)
        candidateAlphaList = nonzero(svm.errorCache[:, 0].A)[
            0]  # mat.A return array
        maxStep = 0
        alpha_j = 0
        error_j = 0
    
        # find the alpha with max iterative step
        if len(candidateAlphaList) > 1:
            for alpha_k in candidateAlphaList:
                if alpha_k == alpha_i:
                    continue
                error_k = calcError(svm, alpha_k)
                if abs(error_k - error_i) > maxStep:
                    maxStep = abs(error_k - error_i)
                    alpha_j = alpha_k
                    error_j = error_k
        # if came in this loop first time, we select alpha j randomly
        else:
            alpha_j = alpha_i
            while alpha_j == alpha_i:
                alpha_j = int(random.uniform(0, svm.numSamples))
            error_j = calcError(svm, alpha_j)
    
        return alpha_j, error_j
    
    
    # the inner loop for optimizing alpha i and alpha j
    def innerLoop(svm, alpha_i):
        error_i = calcError(svm, alpha_i)
    
        # check and pick up the alpha who violates the KKT condition
        # satisfy KKT condition
        # 1) yi*f(i) >= 1 and alpha == 0 (outside the boundary)
        # 2) yi*f(i) == 1 and 0
        # 3) yi*f(i) <= 1 and alpha == C (between the boundary)
        # violate KKT condition
        # because y[i]*E_i = y[i]*f(i) - y[i]^2 = y[i]*f(i) - 1, so
        # 1) if y[i]*E_i < 0, so yi*f(i) < 1, if alpha < C, violate!(alpha = C will be correct)
        # 2) if y[i]*E_i > 0, so yi*f(i) > 1, if alpha > 0, violate!(alpha = 0 will be correct)
        # 3) if y[i]*E_i = 0, so yi*f(i) = 1, it is on the boundary, needless optimized
        if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\
                (svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):
    
            # step 1: select alpha j
            alpha_j, error_j = selectAlpha_j(svm, alpha_i, error_i)
            alpha_i_old = svm.alphas[alpha_i].copy()
            alpha_j_old = svm.alphas[alpha_j].copy()
    
            # step 2: calculate the boundary L and H for alpha j
            if svm.train_y[alpha_i] != svm.train_y[alpha_j]:
                L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])
                H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])
            else:
                L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)
                H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])
            if L == H:
                return 0
    
            # step 3: calculate eta (the similarity of sample i and j)
            eta = 2.0 * svm.kernelMat[alpha_i, alpha_j] - svm.kernelMat[alpha_i, alpha_i] \
                - svm.kernelMat[alpha_j, alpha_j]
            if eta >= 0:
                return 0
    
            # step 4: update alpha j
            svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta
    
            # step 5: clip alpha j
            if svm.alphas[alpha_j] > H:
                svm.alphas[alpha_j] = H
            if svm.alphas[alpha_j] < L:
                svm.alphas[alpha_j] = L
    
            # step 6: if alpha j not moving enough, just return
            if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:
                updateError(svm, alpha_j)
                return 0
    
            # step 7: update alpha i after optimizing aipha j
            svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \
                * (alpha_j_old - svm.alphas[alpha_j])
    
            # step 8: update threshold b
            b1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
                * svm.kernelMat[alpha_i, alpha_i] \
                - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
                * svm.kernelMat[alpha_i, alpha_j]
            b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \
                * svm.kernelMat[alpha_i, alpha_j] \
                - svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \
                * svm.kernelMat[alpha_j, alpha_j]
            if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):
                svm.b = b1
            elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):
                svm.b = b2
            else:
                svm.b = (b1 + b2) / 2.0
    
            # step 9: update error cache for alpha i, j after optimize alpha i, j and b
            updateError(svm, alpha_j)
            updateError(svm, alpha_i)
    
            return 1
        else:
            return 0
    
    
    # the main training procedure
    def trainSVM(train_x, train_y, C, toler, maxIter, kernelOption=('rbf', 1.0)):
        # calculate training time
        startTime = time.time()
    
        # init data struct for svm
        svm = SVMStruct(mat(train_x), mat(train_y), C, toler, kernelOption)
    
        # start training
        entireSet = True
        alphaPairsChanged = 0
        iterCount = 0
        # Iteration termination condition:
        #   Condition 1: reach max iteration
        #   Condition 2: no alpha changed after going through all samples,
        #                in other words, all alpha (samples) fit KKT condition
        while (iterCount < maxIter) and ((alphaPairsChanged > 0) or entireSet):
            alphaPairsChanged = 0
    
            # update alphas over all training examples
            if entireSet:
                for i in xrange(svm.numSamples):
                    alphaPairsChanged += innerLoop(svm, i)
                print '---iter:%d entire set, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
                iterCount += 1
            # update alphas over examples where alpha is not 0 & not C (not on boundary)
            else:
                nonBoundAlphasList = nonzero(
                    (svm.alphas.A > 0) * (svm.alphas.A < svm.C))[0]
                for i in nonBoundAlphasList:
                    alphaPairsChanged += innerLoop(svm, i)
                print '---iter:%d non boundary, alpha pairs changed:%d' % (iterCount, alphaPairsChanged)
                iterCount += 1
    
            # alternate loop over all examples and non-boundary examples
            if entireSet:
                entireSet = False
            elif alphaPairsChanged == 0:
                entireSet = True
    
        print 'Congratulations, training complete! Took %fs!' % (time.time() - startTime)
        return svm
    
    
    # testing your trained svm model given test set
    def testSVM(svm, test_x, test_y):
        test_x = mat(test_x)
        test_y = mat(test_y)
        numTestSamples = test_x.shape[0]
        supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]
        supportVectors = svm.train_x[supportVectorsIndex]
        supportVectorLabels = svm.train_y[supportVectorsIndex]
        supportVectorAlphas = svm.alphas[supportVectorsIndex]
        matchCount = 0
        for i in xrange(numTestSamples):
            kernelValue = calcKernelValue(
                supportVectors, test_x[i, :], svm.kernelOpt)
            predict = kernelValue.T * \
                multiply(supportVectorLabels, supportVectorAlphas) + svm.b
            if sign(predict) == sign(test_y[i]):
                matchCount += 1
        accuracy = float(matchCount) / numTestSamples
        return accuracy
    
    
    # show your trained svm model only available with 2-D data
    def showSVM(svm):
        if svm.train_x.shape[1] != 2:
            print "Sorry! I can not draw because the dimension of your data is not 2!"
            return 1
    
        # draw all samples
        for i in xrange(svm.numSamples):
            if svm.train_y[i] == -1:
                plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'or')
            elif svm.train_y[i] == 1:
                plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'ob')
    
        # mark support vectors
        supportVectorsIndex = nonzero(svm.alphas.A > 0)[0]
        for i in supportVectorsIndex:
            plt.plot(svm.train_x[i, 0], svm.train_x[i, 1], 'oy')
    
        # draw the classify line
        w = zeros((2, 1))
        for i in supportVectorsIndex:
            w += multiply(svm.alphas[i] * svm.train_y[i], svm.train_x[i, :].T)
        min_x = min(svm.train_x[:, 0])[0, 0]
        max_x = max(svm.train_x[:, 0])[0, 0]
        y_min_x = float(-svm.b - w[0] * min_x) / w[1]
        y_max_x = float(-svm.b - w[0] * max_x) / w[1]
        plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')
        plt.show()