線形SVM - 人工知能に関する断創録

下巻に入って7章のサポートベクトルマシン（Support Vector Machine: SVM）を実装してみます。SVMに関しては、有名なSVMのライブラリ（libsvm）を使ったことがあるだけで、アルゴリズム詳細はPRMLで初めて学習しました。なので変なことを書いていたらコメント欄で指摘してもらえると助かります。

まずは、一番簡単な線形SVMを実装してみます。今までと同様に直線（超平面）でデータが完全に分離できる場合です。PRMLの7章には特に説明がありませんが、カーネル関数に下の線形カーネル（データのただの内積）を用いた場合に相当するようです。このカーネル関数を多項カーネルやガウシアンカーネルに変更すると線形分離不可能なデータも分類できるようになるとのこと。非線形SVMは次回ためしてみます。

$k(\bf{x_n}, \bf{x_m}) = \bf{x_n}^T \bf{x_m}$

まず、SVMの識別関数は、式（7.1）で表せます。

$(7.1) \displaystyle \hspace{50pt} y(\bf{x}) = \bf{w}^T \phi(\bf{x}) + b$

今までと違ってバイアスパラメータをまとめないのが流儀なんでしょうか？今までと同様に訓練データと教師信号（クラス1のとき+1、クラス2のとき-1）の組からパラメータwとbを学習するのが目的です。SVMの特徴は、マージン（分類境界と訓練データ間の最短距離）を最大化するようなwとbを求めることです。このとき、wに式（7.5）のような制約を付けます。制約は訓練データの1つにつき1つなのでN個データがあったらN個制約の式があります。

$(7.5) \hspace{50pt} \displaystyle t_n (\bf{w}^T \phi(\bf{x_n}) + b) \ge 1 \hspace{20pt} (n = 1, \cdots, N)$

この制約は、すべてのデータがy>=1かy<=-1にあることを言っています。つまり、-1 < y < 1の間にデータが存在しない空白地帯（マージン）ができます。この空白地帯が最大になるようにしたいってことですね。線形分離可能などんなデータでも式（7.5）の制約を満たすようにwを調節できます（PRMLの図7.1の右側がマージン最大化したイメージ図だと思いますが、左側とy=1、y=-1が入れ替わっているのはなぜなんでしょうね？誤植かな？）。そして、マージン1/|w|を最大化するのですが、これは|w|^2を最小化する問題と等価になり、式（7.5）の制約のもとで式（7.6）を満たすwとbを求めればよいとのこと。

$(7.6) \hspace{50pt} \displaystyle \arg \min_{\bf{w}, b} \quad \frac{1}{2} ||w||^2$

制約付きの極値問題はラグランジュの未定乗数法で解くことができ、ラグランジュ乗数a_nを導入すると、ラグランジュ関数は式（7.7）になります。

$(7.7) \hspace{50pt} \displaystyle L(\bf{w}, b, \bf{a}) = \frac{1}{2} ||w||^2 - \sum_{n=1}^N a_n \{t_n (\bf{w}^T \phi(\bf{x_n}) + b) - 1\}$

これを、wとbで偏微分して0とおく（極値条件）と式（7.8）と式（7.9）の2つの条件が出てきました。

$(7.8) \hspace{50pt} \displaystyle \bf{w} = \sum_{n=1}^N a_n t_n \phi(\bf{x_n})$
$(7.9) \hspace{50pt} \displaystyle 0 = \sum_{n=1}^N a_n t_n$

でここからが面白いところでこの条件を「元の」ラグランジュ関数に代入し、wとbを消去してaの関数にしてしまいます。ここで出てきた式（7.10）を双対表現と言うとのこと。双対表現にするとカーネル関数が出てきます。

$(7.10) \hspace{50pt} \displaystyle \tilde{L}(\bf{a}) = \sum_{n=1}^N a_n - \frac{1}{2} \sum_{n=1}^N \sum_{m=1}^N a_n a_m t_n t_m k(\bf{x_n}, \bf{x_m})$

aが満たすべき制約は、式（7.9）と式（7.11）です。

$(7.11) \hspace{50pt} \displaystyle a_n \ge 0 \hspace{20pt} (n = 1, \cdots, N)$

ここまでやると、wとbに関する最適化問題がaに関する最適化問題に置き換わってしまいます。そして、aに関する最適化問題は、二次計画法（quadratic programming）で解ける形式になっている！そして、aさえ求まればwも決まるので識別関数（7.1）はすぐに求められます。式（7.1）に式（7.8）を代入すると式（7.13）になり、wではなくaを使った識別関数の式に変換できます。

$(7.13) \hspace{50pt} \displaystyle y(\bf{x}) = \sum_{n=1}^N a_n t_n k(\bf{x}, \bf{x_n}) + b$

また、bはサポートベクトルの条件から式（7.18）で求まります。

$(7.18) \hspace{50pt} \displaystyle b = \frac{1}{N_S} \sum_{n \in S} (t_n - \sum_{m \in S} a_m t_m k(\bf{x_n}, \bf{x_m}))$

Sはサポートベクトルの集合、N_Sはサポートベクトルの数です。何はともあれaが求まればwもbも求まって最終的に識別関数も求まるってことですね。じゃ、次はラグランジュ乗数aを二次計画法で解こうって話になります。二次計画法は線形不等式のもとで二次関数を最小化する手法とのこと。解き方はいろいろあるそうで、PRMLでは逐次最小問題最適化法（SMO）ってのがよくつかわれると紹介があります。ただ、説明がほとんどなくこれでは実装できないので、ずるをして（今回はSVMが目的なんで）cvxoptというPythonの最適化ライブラリを使ってみます。SMOかはわかりませんが、二次計画法も実装されています。Quadratic Programmingのドキュメントと例題を読むと、

$\displaystyle \frac{1}{2} \bf{x}^T Q \bf{x} + p^T \bf{x} \rightarrow \min \hspace{20pt} (Gx \le h, Ax = b)$

という二次計画問題を解くには、

     cvxopt.solvers.qp(Q, p, G, h, A, b)

とすればよいとあります（ドキュメントと例題はQとpの使い方が逆になっていますが、例題側にあわせています）。式（7.10）と全然形が違うと思ったのですが、実は下のように置くと同じ形になります。ラグランジュ乗数L(a)は最大化が目的ですが、下のようにおくと-L(a)とマイナスがつくので最小化を目指すことになります。Kはグラム行列ってやつですが、要素に教師信号をかけたt[i]*t[j]がつくので注意です。Nは訓練データ数です。実際にN=2くらいで行列を書き下してみると、下の設定による二次計画法がラグランジュ関数-L(a)と等しいことが確かめられます。

    # ラグランジュ乗数を二次計画法（Quadratic Programming）で求める
    K = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            K[i, j] = t[i] * t[j] * kernel(X[i], X[j])
    Q = cvxopt.matrix(K)
    p = cvxopt.matrix(-np.ones(N))             # -1がN個の列ベクトル
    G = cvxopt.matrix(np.diag([-1.0]*N))       # 対角成分が-1のNxN行列
    h = cvxopt.matrix(np.zeros(N))             # 0がN個の列ベクトル
    A = cvxopt.matrix(t, (1,N))                # N個の教師信号が要素の行ベクトル（1xN）
    b = cvxopt.matrix(0.0)                     # 定数0.0
    sol = cvxopt.solvers.qp(Q, p, G, h, A, b)  # 二次計画法でラグランジュ乗数aを求める
    a = array(sol['x']).reshape(N)             # 'x'がaに対応する
    print a

numpyのarrayではなく、cvxoptのmatrixではないと受け付けません。ただ、arrayとmatrixは相互変換できるように作られているようです。では、ためしてみます。cvxoptを事前にインストールする必要があります。

#coding:utf-8

# 線形SVM
# cvxoptのQuadratic Programmingを解く関数を使用

import numpy as np
from scipy.linalg import norm
import cvxopt
import cvxopt.solvers
from pylab import *

N = 100  # データ数

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

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

if __name__ == "__main__":
    # 訓練データを作成
    cls1 = []
    cls2 = []

    mean1 = [-1, 2]
    mean2 = [1, -1]
    cov = [[1.0,0.8], [0.8, 1.0]]
    
    cls1.extend(np.random.multivariate_normal(mean1, cov, N/2))
    cls2.extend(np.random.multivariate_normal(mean2, cov, N/2))
    
    # データ行列Xを作成
    X = vstack((cls1, cls2))
    
    # ラベルtを作成
    t = []
    for i in range(N/2):
        t.append(1.0)   # クラス1
    for i in range(N/2):
        t.append(-1.0)  # クラス2
    t = array(t)
    
    # ラグランジュ乗数を二次計画法（Quadratic Programming）で求める
    K = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            K[i, j] = t[i] * t[j] * kernel(X[i], X[j])
    Q = cvxopt.matrix(K)
    p = cvxopt.matrix(-np.ones(N))             # -1がN個の列ベクトル
    G = cvxopt.matrix(np.diag([-1.0]*N))       # 対角成分が-1のNxN行列
    h = cvxopt.matrix(np.zeros(N))             # 0がN個の列ベクトル
    A = cvxopt.matrix(t, (1,N))                # N個の教師信号が要素の行ベクトル（1xN）
    b = cvxopt.matrix(0.0)                     # 定数0.0
    sol = cvxopt.solvers.qp(Q, p, G, h, A, b)  # 二次計画法でラグランジュ乗数aを求める
    a = array(sol['x']).reshape(N)             # 'x'がaに対応する
    print a

    
    # サポートベクトルのインデックスを抽出
    S = []
    for i in range(len(a)):
        if a[i] < 0.00001: continue
        S.append(i)
    
    # wを計算
    w = np.zeros(2)
    for n in S:
        w += a[n] * t[n] * X[n]
    
    # bを計算
    sum = 0
    for n in S:
        temp = 0
        for m in S:
            temp += a[m] * t[m] * kernel(X[n], X[m])
        sum += (t[n] - temp)
    b = sum / len(S)
    
    print S, b
    
    # 訓練データを描画
    x1, x2 = np.array(cls1).transpose()
    plot(x1, x2, 'rx')
    
    x1, x2 = np.array(cls2).transpose()
    plot(x1, x2, 'bx')
    
    # サポートベクトルを描画
    for n in S:
        scatter(X[n,0], X[n,1], s=80, c='c', 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()

ラグランジュ乗数aの出力結果は、下のようになります。

Optimal solution found.
[  9.97273590e-10   7.98981786e-10   1.03506919e-09   6.17912793e-10
   9.10675053e-10   9.01468645e-10   5.56458118e-10   1.16859402e-09
   8.43724806e-10   8.30459486e-10   4.41787152e-10   1.13298134e-09
   5.91637823e-10   1.34747211e-09   9.30731884e-10   9.76136081e-10
   9.37096128e-10   1.16101659e-09   1.10479480e-09   8.54622020e-10
   5.88324059e-10   6.99580873e-01   8.99946843e-10   3.56721637e-10
   1.24235147e-09   8.22830144e-10   1.04218873e-09   5.93072343e-10
   1.24675759e-09   1.48237157e-09   1.44087609e-09   4.96985056e-10
   1.01088309e-09   5.02370334e-10   4.16597956e-10   5.86217425e-10
   2.41784674e-09   1.17027497e-09   3.57412003e-09   1.16853725e-09
   1.06772274e-09   7.02814767e-10   4.88276504e-10   6.05880284e-10
   3.44003509e-10   6.74570426e-10   5.84429191e-09   1.82715181e-09
   5.45153517e-10   5.05061091e-10   3.87154248e-10   5.34313746e-09
   4.65807585e-10   4.84783794e-10   8.23285140e-10   5.38684427e-09
   7.53777723e-10   7.99663665e-10   1.65324805e-09   8.29040400e-10
   5.58985700e-10   4.57652358e-10   9.50314068e-10   3.99388053e-10
   8.20185718e-10   9.44976171e-10   5.32702107e-10   3.85934835e-10
   6.85196286e-10   6.13486583e-10   8.27424160e-10   7.09635348e-10
   1.01215519e-09   5.76159099e-10   8.40025594e-10   4.12163546e-10
   1.28867545e-09   8.08776675e-10   1.08172075e-09   7.45955517e-09
   1.52212095e-01   2.82647704e-10   7.49962266e-10   8.29774662e-10
   5.47368766e-01   6.96348345e-10   6.09702674e-10   4.80598085e-09
   5.60964351e-10   9.46656577e-10   6.88250137e-09   3.92380165e-10
   1.20997099e-09   8.90907369e-10   5.75479302e-10   5.38351970e-10
   4.39364577e-10   8.01302244e-10   7.76893463e-10   3.06443105e-09]

ほとんどのデータに対応するラグランジュ乗数はa_n=0になっています（実際はすごく小さな数です）。a_n>0となるラグランジュ乗数に対応するデータx_nをサポートベクトルと呼びます。上の結果では、3つだけ0でないラグランジュ乗数があります。なので、サポートベクトルは3つですね！分類境界は下のようになります。

f:id:aidiary:20100501195851p:plain

上の結果で水色の丸がa_n>0のサポートベクトルです。式（7.13）の識別関数からわかりますが、サポートベクトル以外のデータ（a_n=0）は識別関数にまったく影響しません。サポートベクトルのみで分類境界が決まります（これが疎な解の意味？）。疎な解だと何がうれしいのかPRMLには書かれてないようなんですが、ノイズに対して頑強で、サポートベクトルだけ予測できるので使用メモリも少なくてよいとかでしょうか？ためしに最小二乗法（2010/4/24）でうまくいかなかった例題を試してみます。

f:id:aidiary:20100501200612p:plain