サポートベクターマシンを勾配法で解く。

最も代表的な教師あり分類手法の一つであるサポートベクターマシン。
訓練サンプルから、各データ点との距離が最大となるマージン最大化超平面を求めるという明確な分類基準がある、かつ、一見求めづらそうなこの超平面を一意に求められるということから
数学的に美しい分類手法と言えるでしょう。
サポートベクターマシンの詳しい説明はこちらになります。
サポートベクターマシン - Wikipedia

さて、上のWikipediaの記述にもある通り、
マージン最大化超平面は最適化問題の凸二次計画問題となります。
今回は、この凸二次計画問題を勾配法で解くことで、
ハードマージンのサポートベクターマシンを実装しました。

使用したデータ

線形分離可能なデータです。

また、こちらの座標関係をテキストデータにしたものがこちらです。 各行がデータを表しており、スペースで区切られています。
左から、X座標、Y座標、クラス(1or-1)となります。

1 2 1
0 1 1
-1 2 1
2 1 1
3 5 1
0 -1 -1
-1 -2 -1
-2 0 -1
-4 0 -1
-3 -3 -1
(data/svm_hard_data.txtの中身。プログラム中で使用。)

実装したコード

本記事のプログラムは
線形SVM - 人工知能に関する断創録
を大いに参考にしました。
言語はpythonです。

作成したプログラム

#coding:utf-8

import numpy as np
from pylab import *

LR = 0.05    # 学習率
C = sys.maxint    # Cを無限大に設定
CountMax = 1000   # 1000回更新

def f(x1, w, b):
    return - (w[0] / w[1]) * x1 - (b / w[1])

def kernel(x, y):
    return np.dot(x, y)  # 線形カーネル

def dL(i):
    ans = 0
    for j in range(0,N):
        ans += L[j] * t[i] * t[j] * kernel(X[i], X[j])
    return (1 - ans)

if __name__ == "__main__":
    data = np.loadtxt("data/svm_hard_data.txt")
    N = data.size / 3                   # 教師データ数
    X = np.delete(data,2,1)             # 入力の座標
    X = np.c_[X,np.ones(X.shape[0])]    # 入力の空間を1次元拡張
    t = np.array(data[:,2])             # 入力のラベル
    L = np.zeros((N,1))                 # データの個数分のラグランジュ乗数

    count = 0
    while (count < CountMax):
        for i in range(N):
            L[i] = L[i] + LR * dL(i)    # ラグランジュ乗数の更新
            if (L[i] < 0):
                L[i] = 0
            elif (L[i] > C):
                L[i] = C
        count += 1

    for i in range(N):
        print L[i]

    # サポートベクトルのインデックスを抽出
    S = []
    for i in range(len(L)):
        if L[i] < 0.00001: continue
        S.append(i)

    # wを計算
    w = np.zeros(3)
    for n in S:
        w += L[n] * t[n] * X[n]

    # wの3次元目は拡張次元のbとなる。
    b = w[2]
    np.delete(w, 2, 0)

    # 訓練データを描画
    for i in range(0,N):
        if(t[i] > 0):
            plot(X[i][0],X[i][1], 'rx')
        else:
            plot(X[i][0],X[i][1], 'bx')

    # サポートベクトルを描画
    for n in S:
        scatter(X[n,0], X[n,1], s=80, c='y', marker='o')

    # 識別境界を描画
    x1 = np.linspace(-6, 6, 1000)
    x2 = [f(x, w, b) for x in x1]
    plot(x1, x2, 'g-')

    xlim(-6, 6)
    ylim(-6, 6)
    show()

出力結果

黄色の○で囲まれている座標がサポートベクトルです。 また、このときのラグランジュ乗数は以下になります。

[ 0.]
[ 0.625]
[ 0.]
[ 0.]
[ 0.]
[ 0.375]
[ 0.]
[ 0.25]
[ 0.]
[ 0.]

このように、サポートベクトル以外のラグランジュ乗数が0になっていることが確認できます。

理解に必要な数学

• 勾配法
• ラグランジュ乗数
• クーン・タッカー条件
• 特徴空間を1次元拡張（バイアスを固定）

以上です。
これらが理解できていれば、今回のコードの理解もできるでしょう。

ただ、余談ですが、これらの数学的理論についてのわかりやすい説明が身近にないように感じています。
どれも幾何学的に説明すれば、理解が容易だと思うのですが……。
が、図を書くのって正直面倒ですよね、わかります。
本当はこれら3つの数学的理論についても説明したいのですが、
私も泣く泣く割愛致します。

勾配法によるSVM

まずは、参考文献を掲載します。

サポートベクターマシン入門

こちらの本に勾配法よるSVMの実装の詳細が書かれています。

大枠はこちらの本による説明に任せるとして、
今回は勾配法によるSVMの肝となる処理であるたった1行だけを軽く触れます。

L[i] = L[i] + LR * dL(i)    # ラグランジュ乗数の更新

• L[i]は（現在の）ラグランジュ乗数。
• LRは学習率。
• dL(i)はラグランジュ関数を微分。

となっています。

参考までに、式を載せておきます。

あとはこの式に従って、ひたすら（サンプルだと1000回）回すだけで、
マージンが最大化となる超平面を求めることができます。

数学的に美しいですよね。素晴らしいです。