![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
|
2011/06/28からのアクセス回数 5005 ここで紹介したSageワークシートは、以下のURLからダウンロードできます。 http://sage.math.canterbury.ac.nz/home/pub/105/ また、Sageのサーバを公開しているサイト（http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。 第7章-SVM（サポートベクトルマシン）をSageで試す †今回は、SVMに関する他のサイトの情報を元にSageでSVMを試してみます。 参考（ほとんどコピー）にしたのは以下のサイトです。
上記で説明されたSVMを計算する2通りの方法が7章を理解する上で助けになりました。 テストデータ †SVMのテストには、PRMLの図7.4、
で使用された（赤と青の点が入り混ざっている）人工データを使います。 カーネルには、\(exp(-\gamma||x - x'||^2)\)のガウスカーネルを使い、\(\gamma=0.45\)とします。 sageへの入力: # SMOで識別を行う from numpy import loadtxt data = loadtxt(DATA+"classification.txt") # 1,2カラムはx, y、3カラムには、0:red cosses×, 1:blue circle ○がセットされている # これをx, yと識別子(-1, 1)に変換する X = [[x[0], x[1]] for x in data[:,0:2]] T = data[:,2]*2-1 sageへの入力: var('x y')
data_plt = Graphics()
for i in range(len(T)):
if T[i] == -1:
data_plt += point(X[i], rgbcolor='red')
else:
data_plt += point(X[i], rgbcolor='blue')
data_plt.show()
SMO †「きちめも」で実装されている、「サポートベクターマシン入門」の付録に掲載されているSMO(Sequenctial Minimal Optimization) は、J. Platt(1988)で発表されたもので、2つのデータ点に対する最適化問題を解析的に求め、ヒューリスティクスな選択で反復的な処理を 高速に収束させる手法です。 以下では、「きちめも」の実装に若干の修正を加えています。
sageへの入力: # se-kichiさんの「きちめも」から引用
# http://d.hatena.ne.jp/se-kichi/20100306/1267858745#20100306f1
# 修正点： _searchのelse jをelse kに変更
# 「サポートベクターマシン入門」の注釈7.14からcalc_bを_examinExの最後にも追加
# KKTのチェックにself._epsが含まれているが、「サポートベクターマシン入門」の付録Aに合わせた
import random
from scipy import zeros, array
from scipy.linalg import norm as _norm
class SMO:
"""
Support Vector Machine
using SMO Algorithm.
"""
def __init__(self,
kernel=lambda x,y:dot(x,y),
c=10000,
tol=1e-2,
eps=1e-2,
loop=float('inf')):
"""
Arguments:
- `kernel`: カーネル関数
- `c`: パラメータ
- `tol`: KKT条件の許容する誤差
- `eps`: αの許容する誤差
- `loop`: ループの上限
"""
self._kernel = kernel
self._c = c
self._tol = tol
self._eps = eps
self._loop = loop
def _takeStep(self, i1, i2):
if i1 == i2:
return False
alph1 = self._alpha[i1]
alph2 = self._alpha[i2]
y1 = self._target[i1]
y2 = self._target[i2]
e1 = self._e[i1]
e2 = self._e[i2]
s = y1 * y2
if y1 != y2:
L = max(0, alph2 - alph1)
H = min(self._c, self._c-alph1+alph2)
else:
L = max(0, alph2 + alph1 - self._c)
H = min(self._c, alph1+alph2)
if L == H:
return False
k11 = self._kernel(self._point[i1], self._point[i1])
k12 = self._kernel(self._point[i1], self._point[i2])
k22 = self._kernel(self._point[i2], self._point[i2])
eta = 2 * k12 - k11 - k22
if eta > 0:
return False
a2 = alph2 - y2 * (e1 - e2) / eta
a2 = min(H, max(a2, L))
if abs(a2 - alph2) < self._eps * (a2 + alph2 + self._eps):
return False
a1 = alph1 + s * (alph2 - a2)
# update
da1 = a1 - alph1
da2 = a2 - alph2
self._e += array([(da1 * self._target[i1] * self._kernel(self._point[i1], p) +
da2 * self._target[i2] * self._kernel(self._point[i2], p))
for p in self._point])
self._alpha[i1] = a1
self._alpha[i2] = a2
return True
def _search(self, i, lst):
if self._e[i] >= 0:
return reduce(lambda j,k: j if self._e[j] < self._e[k] else k, lst)
else:
return reduce(lambda j,k: j if self._e[j] > self._e[k] else k, lst)
def _examinEx(self, i2):
y2 = self._target[i2]
alph2 = self._alpha[i2]
e2 = self._e[i2]
r2 = e2*y2
if ((r2 < -self._tol and alph2 < self._c) or
(r2 > self._tol and alph2 > 0)):
alst1 = [i for i in range(len(self._alpha))
if 0 < self._alpha[i] < self._c]
if alst1:
i1 = self._search(i2, alst1)
if self._takeStep(i1, i2):
return True
random.shuffle(alst1)
for i1 in alst1:
if self._takeStep(i1, i2):
return True
alst2 = [i for i in range(len(self._alpha))
if (self._alpha[i] <= 0 or
self._alpha[i] >= self._c)]
random.shuffle(alst2)
for i1 in alst2:
if self._takeStep(i1, i2):
return True
# 注釈7.14にtakeStepでbの更新は必要ないが、停止基準を評価する場合には必要とのコメントから挿入
self._calc_b()
return False
def _calc_b(self):
self._s = [i for i in range(len(self._target))
if 0 < self._alpha[i]]
self._m = [i for i in range(len(self._target))
if 0 < self._alpha[i] < self._c]
self._b = 0.0
for i in self._m:
self._b += self._target[i]
for j in self._s:
self._b -= (self._alpha[j]*self._target[j]*
self._kernel(self._point[i], self._point[j]))
self._b /= len(self._m)
def calc(self, x):
ret = self._b
for i in self._s:
ret += (self._alpha[i]*self._target[i]*
self._kernel(x, self._point[i]))
return ret
def learn(self, point, target):
self._target = target
self._point = point
self._alpha = zeros(len(target), dtype=float)
self._b = 0
self._e = -1*array(target, dtype=float)
changed = False
examine_all = True
count = 0
while changed or examine_all:
count += 1
if count > self._loop:
break
changed = False
if examine_all:
for i in range(len(self._target)):
changed |= self._examinEx(i)
else:
for i in (j for j in range(len(self._target))
if 0 < self._alpha[j] < self._c):
changed |= self._examinEx(i)
if examine_all:
examine_all = False
elif not changed:
examine_all = True
self._calc_b()
def _get_s(self):
return self._s
def _get_m(self):
return self._m
def _get_alpha(self):
return self._alpha
sageへの入力: # numpyのnormは、ベクトルの絶対値を返すから自乗する
import numpy as np
N = len(T)
C = 30
def kernel(x, y):
return np.exp(-0.45*_norm(x-y)**2)
s = SMO(c=C, kernel=kernel, loop=20)
p = data[:,0:2]
t = data[:,2]*2-1.0
s.learn(p, t)
SMOの計算結果 †SMOの計算結果は、図7.4と同様のパターンを得ました。
sageへの入力: # コンター図を作成
cnt_plt = contour_plot(lambda x, y : s.calc([x, y]), [x, -2.5, 2.5], [y, -3, 3], contours=(0.0,), fill=False)
# サポートベクターのポイントを強調して表示
sv_plt = Graphics()
for i in s._get_s():
if T[i] == -1:
sv_plt += point(X[i], rgbcolor='red', pointsize=15, faceted=True)
else:
sv_plt += point(X[i], rgbcolor='blue', pointsize=15, faceted=True)
(cnt_plt + data_plt + sv_plt).show()
sageへの入力: # コンター図を作成 cnt_plt = contour_plot(lambda x, y : s.calc([x, y]), [x, -2.5, 2.5], [y, -3, 3], fill=False, cmap='hsv') (cnt_plt + data_plt + sv_plt).show()
CVXOPTを使ったSVM †「人工知能に関する断想録」では、適化ライブラリCVXOPTを使ってSVMを計算しています。 「人工知能に関する断想録」のポイントは、CVXOPTの入力形式に $$ \begin{array}{ll} \mbox{minimize} & (1/2)x^TPx + q^Tx \\ \mbox{subject to} & G x \leq h \\ & A x = b \end{array} $$ 制約条件\(0 \le a_n \le C\)をどうやって組み入れるかですが、 $$ G x \leq h $$ を $$ \begin{pmatrix} -1 & 0 \\ 0 & -1 \\ 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} a_1\\ a_2 \end{pmatrix} \le \begin{pmatrix} 0 \\ 0 \\ C \\ C \end{pmatrix} $$ とすることで解決されたそうです。 「人工知能に関する断想録」からの変更点は
sageへの入力: # 同様の計算をCVXOPTを使って計算する例 # 「人工知能に関する断想録」から引用 # http://d.hatena.ne.jp/aidiary/20100503/1272889097 X = data[:,0:2] t = data[:,2]*2-1.0 sageへの入力: from cvxopt import solvers from cvxopt.base import matrix as _matrix sageへの入力: 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
sageへの入力: # ラグランジュ乗数を二次計画法（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 = _matrix(K)
p = _matrix(-np.ones(N))
temp1 = np.diag([-1.0]*N)
temp2 = np.identity(N)
G = _matrix(np.vstack((temp1, temp2)))
temp1 = np.zeros(N)
temp2 = np.ones(N) * C
h = _matrix(np.hstack((temp1, temp2)))
A = _matrix(t, (1,N))
b = _matrix(np.array([[1.0]]))
sol = solvers.qp(Q, p, G, h, A, b)
a = array(sol['x']).reshape(N)
#print a
# サポートベクトルのインデックスを抽出
S = []
M = []
# aidiaryさんのコメントを反映して、αの誤差としてepsに1e-3を使用するように変更
eps = 1e-3
for n in range(len(a)):
if eps < a[n]:
S.append(n)
if eps < a[n] < C:
M.append(n)
# bを計算
sum = 0.0
for n in M:
sum = t[n]
for m in S:
sum -= a[m] * t[m] * kernel(p[n], p[m])
b = sum / len(M)
#print S, b
sageからの出力: pcost dcost gap pres dres 0: 0.0000e+00 -6.0000e+03 4e+02 1e+00 1e+00 1: -6.4529e+01 -2.3215e+03 3e+02 1e+00 1e+00 2: -1.8484e+02 -1.0747e+03 2e+02 9e-01 9e-01 3: -3.7697e+02 -1.0604e+03 2e+02 8e-01 8e-01 4: -6.3559e+02 -1.3224e+03 1e+02 7e-01 7e-01 5: -8.8941e+02 -1.6107e+03 1e+02 7e-01 7e-01 6: -1.1981e+03 -1.9385e+03 1e+02 6e-01 6e-01 7: -1.4595e+03 -2.1869e+03 9e+01 5e-01 5e-01 8: -1.7532e+03 -2.4195e+03 7e+01 4e-01 4e-01 9: -1.9720e+03 -2.5556e+03 6e+01 3e-01 3e-01 10: -2.2906e+03 -2.6936e+03 4e+01 2e-01 2e-01 11: -2.5358e+03 -2.7486e+03 3e+01 8e-02 8e-02 12: -2.6860e+03 -2.7580e+03 1e+01 2e-02 2e-02 13: -2.7456e+03 -2.7557e+03 4e+00 3e-03 3e-03 14: -2.7530e+03 -2.7541e+03 6e-01 2e-04 2e-04 15: -2.7537e+03 -2.7539e+03 1e-01 2e-05 2e-05 16: -2.7538e+03 -2.7538e+03 5e-03 1e-06 4e-14 17: -2.7538e+03 -2.7538e+03 1e-04 5e-08 4e-14 CVXOPTの計算結果 †CVXOPTの計算は、SMOよりも高速で解も安定しています。 CVXOPTの計算結果も、図7.4と同様のパターンを得ました。
sageへの入力: qp_plt = contour_plot(lambda x, y : f([x, y], a, t, X, b), [x, -3, 3], [y, -3, 3], fill=False, cmap='hsv') sageへの入力: # サポートベクターのポイントを強調して表示
sv2_plt = Graphics()
for i in S:
if T[i] == -1:
sv2_plt += point(X[i], rgbcolor='red', pointsize=15, faceted=True)
else:
sv2_plt += point(X[i], rgbcolor='blue', pointsize=15, faceted=True)
sageへの入力: (qp_plt + data_plt + sv2_plt).show()
コメント †皆様のご意見、ご希望をお待ちしております。
Tweet |
sage0-3.png
sage0-2.png
sage0-1.png
sage0.png
