非線形SVM
今回は、非線形サポートベクトルマシンを試してみます。線形SVM(2010/5/1)は、カーネル関数に線形カーネル(ただの内積)を使いましたが、これを多項式カーネル(A)やガウスカーネル(B)に変更します。
カーネル関数は元のベクトルxを非線形写像によって高次元空間に写像した特徴ベクトルφ(x)の内積(C)で定義されます。
一般に特徴ベクトルφ(x)は高次元空間(無限次元空間でもOK)になるので普通にやってたら内積の計算量が非常に大きくなります。そこで、特徴ベクトルφ(x)の内積を計算せずに多項式カーネル(A)やガウスカーネル(B)の計算で置き換えるテクニックをカーネルトリックと呼ぶとのこと。多項式カーネルやガウスカーネルを使うとφ(x)を陽に計算する必要がなくなります。ただ、元の空間xでの内積は必要なんですよね・・・最初は、カーネルトリックのありがたみがよくわからなかったのですが、「入力空間xの次元 << 特徴空間φ(x)の次元」という前提があったようです。そういえば、非線形SVMの基本的な考え方は、入力空間xでは線形分離不可能でも高次元空間φ(x)に写してやれば線形分離可能になることを期待した手法でした。カーネルの理論は実のところまだよくわかってません・・・もっと深耕する必要があるなぁ。
では、プログラムを書いて試してみます。Pythonのグラフ描画ライブラリmatplotlibで等高線を描く方法がよくわからなかったので下のサイトを参考にさせていただきました。下のサイトでは二次計画法のライブラリとしてcvxoptの代わりにOpenOptという最適化ライブラリを使ってます。
#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 # データ数 P = 3 # 多項式カーネルのパラメータ SIGMA = 5.0 # ガウスカーネルのパラメータ # 多項式カーネル 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 # Sを渡してサポートベクトルだけで計算した方が早い # サポートベクトルはa[n]=0なのでsumに足す必要ない 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__": # 訓練データを作成 cls1 = [] cls2 = [] mean1 = [-1, 2] mean2 = [1, -1] mean3 = [4, -4] mean4 = [-4, 4] cov = [[1.0,0.8], [0.8, 1.0]] cls1.extend(np.random.multivariate_normal(mean1, cov, N/4)) cls1.extend(np.random.multivariate_normal(mean3, cov, N/4)) cls2.extend(np.random.multivariate_normal(mean2, cov, N/4)) cls2.extend(np.random.multivariate_normal(mean4, cov, N/4)) # データ行列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)) G = cvxopt.matrix(np.diag([-1.0]*N)) h = cvxopt.matrix(np.zeros(N)) 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 = [] for n in range(len(a)): if a[n] < 1e-5: continue S.append(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='g', marker='o') # 識別境界を描画 X1, X2 = meshgrid(linspace(-6,6,50), linspace(-6,6,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(-6, 6) ylim(-6, 6) show()
データは、赤がクラス1、青がクラス2です。見ればわかるように元の空間xでは、直線「1本」では分離できません(線形分離不可能)。まず、多項式カーネルでパラメータはp=3としてみました。
高次元の特徴空間における超平面の分類境界を入力空間に戻すと上みたいなうねうねした曲線になるんですねー。次に、ガウスカーネルを使ってみます。パラメータはσ=5.0です。
こっちは少し違いますね。面白いです。パラメータをいろいろ変えて試してみるとよいかも。一般に、カーネル関数の最良のパラメータは実験で求めるしかないようです。有名なSVMライブラリlibsvmのA practical guide to SVM classification(PDF)を読むと最適パラメータはグリッドサーチ(パラメータの組み合わせを網羅的に試す)で見つけようと書いてあります。パラメータを変えるとどのように分類境界が変わるかは次回試して見ます。
最後に、PRMLの図7.4のデータを使って試してみましょう。下のプログラムを実行するには、PRMLの原書サポートページにあるclassification.txtというデータを同じディレクトリにおいてください。データはファイルから読み込みます。
#coding:utf-8 # 非線形SVM # cvxoptのQuadratic Programmingを解く関数を使用 # 図7.4をハードマージンSVMで解いた場合 # 解けないで当たり前?ソフトマージンにしないとダメ!!!! # 参考:http://d.hatena.ne.jp/se-kichi/20100306/1267871768 import numpy as np from scipy.linalg import norm import cvxopt import cvxopt.solvers from pylab import * N = 200 # データ数 P = 3 # 多項式カーネルのパラメータ SIGMA = 5.0 # ガウスカーネルのパラメータ # 多項式カーネル 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 # Sを渡してサポートベクトルだけで計算した方が早い # サポートベクトルはa[n]=0なのでsumに足す必要ない 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)) G = cvxopt.matrix(np.diag([-1.0]*N)) h = cvxopt.matrix(np.zeros(N)) 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 = [] for n in range(len(a)): if a[n] < 1e-5: continue S.append(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 # 訓練データを描画 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()
結果は、下のようになります。
あれ?うまく分類できてないですねー。でも分類できなくて当たり前。PRMLの7.1.1 重なりのあるクラス分布に書いてありますが、今までやってきたのはハードマージンSVMといって入力空間xでデータを完全に分離できることを仮定してきた手法だからです。上のデータを見ると2つのクラスのデータが入り組んでいて曲線でも完全に分離できません。少しくらいデータが重なっていて誤分類も許すようなSVMは、ソフトマージンSVM(C-SVM)と呼ぶそうです。この拡張は次回試してみますが、修正は驚くほど簡単です。