ソフトマージンSVM - 人工知能に関する断創録

前回（2010/5/2）のハードマージンSVMでは、データに重なりがある場合、下のようにちゃんと分類境界を求められませんでした。今回は、重なりのあるクラス分布に対応できるように拡張してみます。このようなSVMはハードマージンSVMに対してソフトマージンSVMと呼ばれます。別名としてC-SVMとも呼ばれるようです。

f:id:aidiary:20100502213243p:plain

PRMLの7.1.1にあるように、データの誤分類を許すようにSVMを修正します。ハードマージンSVMでは、データ点がマージン内（-1 < y < 1）に絶対に入らないことを前提にしていましたが、ソフトマージンSVMでは「入ってしまったものは仕方ない、だがペナルティを与える！」と少し条件を緩めます。

まず、スラック変数ζ（ゼータ）をデータごとに導入します。スラック変数は、データが正しく分類されかつマージン境界上または外側にある場合は0、正しく分類されているがマージン内に侵入してしまっているデータは0 < ζ <= 1、分類境界を越えて誤分類されてしまったデータはζ > 1となります。ハードマージンSVMでは、後の2つのようにマージン内に侵入したり、誤分類があることは仮定してなかったわけですね！で、マージン最大化するわけですが、上記のスラック変数をペナルティとして与えます。

$\displaystyle \zeta_n = |t_n - y(\bf{x_n})| = |t_n - (\bf{w}^T \phi (\bf{x}) + b)|$
$(7.21) \hspace{50pt} \displaystyle C \sum_{n=1}^N \zeta_n + \frac{1}{2} ||\bf{w}||^2 \rightarrow \min$

マージンは、1/|w|より上の式は最大化ではなく最小化を目指すので注意です。ここで、C > 0はスラック変数を用いて表されるペナルティと、マージンの大きさの間のトレードオフを制御するパラメータです。このパラメータのCからとってソフトマージンSVMがC-SVMと呼ばれるのだと思います。ν-SVM（ニューSVM）というのもあるんですが、このνもパラメータなので。誤分類が多いほど（7.21）の第1項スラック変数の合計は大きくなるので最小化するのに妨げとなります。なので、最小化しようとするとなるべく誤分類が少なくなるようにパラメータwを調整しようとするはずですね！（ζの定義にもwが含まれてます）正則化もそうですが誤差関数に項を付け加えていいようにしようってのは自分にとって新鮮な考え方でした。最適化恐るべし。ちなみにハードマージンSVMはC→∞のときです。ペナルティがものすごく大きいので誤りは絶対に許せない！ってことですね。これは後で試してみましょう。

ハードマージンSVMのときと同様にラグランジュ関数をたててラグランジュ乗数aを用いた双対表現に変換するわけですが、双対表現が何とハードマージンSVMのときとまったく同じになってしまいます。

$(7.32) \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})$

実際に導出してみますと、スラック変数ζが入った項は+と-できれいに相殺されて全部消えてしまいます。じゃ、ハードマージンSVMと同じじゃないかと思ったのですが、制約条件が違うんですねー。ソフトマージンSVMのラグランジュ乗数aの制約条件は、

$(7.33) \hspace{50pt} \displaystyle 0 \le a_n \le C$
$(7.34) \hspace{50pt} \displaystyle \sum_{n=1}^N a_n t_n = 0$

です。ハードマージンSVMの制約条件は（7.11）（7.12）でした。

$(7.11) \hspace{50pt} \displaystyle a_n \ge 0$
$(7.12) \hspace{50pt} \displaystyle \sum_{n=1}^N a_n t_n = 0$

何が違うかというとラグランジュ乗数にC以下という条件が加わっただけですね。あれだけスラック変数云々と導入してきたのに結局パラメータCしか残らないなんて。何かすごいです。では、早速このパラメータCのすごさをプログラミングして試してみます。

#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 = 200         # データ数
C = 0.5         # スラック変数を用いて表されるペナルティとマージンのトレードオフパラメータ
SIGMA = 0.3     # ガウスカーネルのパラメータ

# 多項式カーネル
def polynomial_kernel(x, y):
    return (1 + np.dot(x, y)) ** P

# ガウスカーネル
def gaussian_kernel(x, y):
    return np.exp(-norm(x-y)**2 / (2 * (SIGMA ** 2)))

# どちらのカーネル関数を使うかここで指定
kernel = gaussian_kernel

def f(x, a, t, X, b):
    sum = 0.0
    for n in range(N):
        sum += a[n] * t[n] * kernel(x, X[n])
    return sum + b

if __name__ == "__main__":
    # 訓練データをロード
    data = np.genfromtxt("classification.txt")
    X = data[:,0:2]
    t = data[:,2] * 2 - 1.0  # 教師信号を-1 or 1に変換
    
    # ラグランジュ乗数を二次計画法（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))
    temp1 = np.diag([-1.0]*N)
    temp2 = np.identity(N)
    G = cvxopt.matrix(np.vstack((temp1, temp2)))
    temp1 = np.zeros(N)
    temp2 = np.ones(N) * C
    h = cvxopt.matrix(np.hstack((temp1, temp2)))
    A = cvxopt.matrix(t, (1,N))
    b = cvxopt.matrix(0.0)
    sol = cvxopt.solvers.qp(Q, p, G, h, A, b)
    a = array(sol['x']).reshape(N)
    print a
    
    # サポートベクトルのインデックスを抽出
    S = []
    M = []
    for n in range(len(a)):
        if 0 < a[n]:
            S.append(n)
        if 0 < a[n] < C:
            M.append(n)
    
    # bを計算
    sum = 0
    for n in M:
        temp = 0
        for m in S:
            temp += a[m] * t[m] * kernel(X[n], X[m])
        sum += (t[n] - temp)
    b = sum / len(M)
    
    print S, b
    
    # 訓練データを描画
    for n in range(N):
        if t[n] > 0:
            scatter(X[n,0], X[n,1], c='b', marker='o')
        else:
            scatter(X[n,0], X[n,1], c='r', marker='o')
    
    # サポートベクトルを描画
#    for n in S:
#        scatter(X[n,0], X[n,1], s=80, c='c', marker='o')
    
    # 識別境界を描画
    X1, X2 = meshgrid(linspace(-2,2,50), linspace(-2,2,50))
    w, h = X1.shape
    X1.resize(X1.size)
    X2.resize(X2.size)
    Z = array([f(array([x1, x2]), a, t, X, b) for (x1, x2) in zip(X1, X2)])
    X1.resize((w, h))
    X2.resize((w, h))
    Z.resize((w, h))
    CS = contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower')
    
    for n in S:
        print f(X[n], a, t, X, b)
    
    xlim(-2, 2)
    ylim(-2, 2)
    show()

サポートベクトルとbの求め方もハードマージンSVMと少し違います。詳しくはPRMLを参照してください。また、ちょっとした実装上のテクニックですが、最適化ライブラリcvxoptの二次計画法の制約条件はGa <= hのように1つの行列の形式でまとめる必要があります。0 < a < Cを1つの行列でどうやって書けばよいか少し悩んだのですが、かなり強引に2つの行列をくっつける形で上のように書きました。たとえば、データ数N=2の場合には、下のようにGとhを決めるとほしかった4つの不等式が得られます。上のプログラムではtemp1とtemp2という上半分の行列と下半分の行列を積み重ねてGとhを作ってます。これでちゃんと制約されているのかcvxoptの中身を見てないのでわからないのですが結果が正しいようなのでOKなのでしょう。

f:id:aidiary:20100503203247p:plain

では、結果です。ガウシアンカーネルを用い、パラメータはC=0.5です。ガウシアンカーネルのパラメータσは上から0.1, 0.3, 1.0です。

f:id:aidiary:20100503203248p:plain
f:id:aidiary:20100503203249p:plain
f:id:aidiary:20100503203250p:plain

どうやらガウシアンカーネルのパラメータσを大きくすると分類境界が滑らかになるようですね。これは、ガウス分布の分散が大きいほど滑らかになるイメージと直感的にも合っているようです。では、次にガウシアンカーネルのパラメータはσ=0.3で固定してCを20, 1000, 1e+10としてみます。

f:id:aidiary:20100503203344p:plain
f:id:aidiary:20100503203345p:plain
f:id:aidiary:20100503203346p:plain

うーむ。Cを大きくすると誤分類のペナルティが大きくなるせいか、「分類境界がどんなに細かくなってもいいから何とか分類するんだ！」という意地が伝わってきますね（笑）実際は、細かすぎると汎化性能が落ちるのでペナルティを大きくするのも考えものだと思いますが・・・これは、「ペナルティが大きいために目先の成果にとらわれると大局を見失ってしまう」というよい教訓になりそうです。まあ、冗談はさておき、C→∞のときハードマージンSVMになるそうですが、完全に分離しようとするという意味ではそうなってそうですね。実際、C=1e+10では完全に分離できています。

調整しなくちゃいけないパラメータは（C, σ）です。非線形SVM（2010/5/2）でも書きましたが、最適なパラメータは、問題にあわせて実験によって求めるのが普通だそうです。グリッドサーチといってCとσの組み合わせのグリッドをつくり網羅的に分類精度を評価するのがよいとのこと。

一応、サポートベクトルマシンは今回が最後です。二次計画法以外は自分で実装してみましたが、普通は車輪の再発明はせずに既存のライブラリを使います（もちろん、SVMのアルゴリズム自体を研究する場合は別ですが）。朱鷺の杜にいろいろなライブラリへのリンクがあるので参考にしてください。今度は、テキストや画像など具体的な問題に応用してみたいと思います。乞うご期待！