- 追加された行はこの色です。
- 削除された行はこの色です。
#freeze
[[FrontPage]]
#contents
2011/06/28からのアクセス回数 &counter;
ここで紹介した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で試す [#a31f14d0]
今回は、SVMに関する他のサイトの情報を元にSageでSVMを試してみます。
参考(ほとんどコピー)にしたのは以下のサイトです。
-
se-kichiさんの
[[「きちめも」>http://d.hatena.ne.jp/se-kichi/20100306/1267858745#20100306f1]]
-
[[「人工知能に関する断想録」>http://d.hatena.ne.jp/aidiary/20100503/1272889097]]
上記で説明されたSVMを計算する2通りの方法が7章を理解する上で助けになりました。
** テストデータ [#x95a1b9a]
SVMのテストには、PRMLの図7.4、
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure7.4.jpg,center,40%);
で使用された(赤と青の点が入り混ざっている)人工データを使います。
カーネルには、\(exp(-\gamma||x - x'||^2)\)のガウスカーネルを使い、\(\gamma=0.45\)とします。
sageへの入力:
#pre{{
# 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への入力:
#pre{{
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()
}}
&ref(sage0.png);
** SMO [#pd9704b1]
「きちめも」で実装されている、「サポートベクターマシン入門」の付録に掲載されているSMO(Sequenctial Minimal Optimization)
は、J. Platt(1988)で発表されたもので、2つのデータ点に対する最適化問題を解析的に求め、ヒューリスティクスな選択で反復的な処理を
高速に収束させる手法です。
以下では、「きちめも」の実装に若干の修正を加えています。
- ガウスカーネルを図7.4の説明に合わせた
- 「サポートベクターマシン入門」の注釈7.14からcalc_bを_examinExの最後にも追加
- KKTのチェックは、「サポートベクターマシン入門」の付録Aに合わせた。
sageへの入力:
#pre{{
# 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への入力:
#pre{{
# 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の計算結果 [#iacaff7e]
SMOの計算結果は、図7.4と同様のパターンを得ました。
- 最初の図は、境界線(黒線)とサポートベクタ(縁取りした大きな点)
- コンター図で分布パターンを示したもの
sageへの入力:
#pre{{
# コンター図を作成
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()
}}
&ref(sage0-1.png);
sageへの入力:
#pre{{
# コンター図を作成
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()
}}
&ref(sage0-2.png);
** CVXOPTを使ったSVM [#p668efbe]
「人工知能に関する断想録」では、適化ライブラリ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}
$$
とすることで解決されたそうです。
「人工知能に関する断想録」からの変更点は
- 「人工知能に関する断想録」のコメントからサポートベクタの抽出条件を
\( \epsilon \le a_n\)に変更しました。
sageへの入力:
#pre{{
# 同様の計算をCVXOPTを使って計算する例
# 「人工知能に関する断想録」から引用
# http://d.hatena.ne.jp/aidiary/20100503/1272889097
X = data[:,0:2]
t = data[:,2]*2-1.0
}}
sageへの入力:
#pre{{
from cvxopt import solvers
from cvxopt.base import matrix as _matrix
}}
sageへの入力:
#pre{{
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への入力:
#pre{{
# ラグランジュ乗数を二次計画法(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からの出力:
#pre{{
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の計算結果 [#z35612d6]
CVXOPTの計算は、SMOよりも高速で解も安定しています。
CVXOPTの計算結果も、図7.4と同様のパターンを得ました。
- 最初の図は、境界線(黒線)とサポートベクタ(縁取りした大きな点)
- コンター図で分布パターンを示したもの
sageへの入力:
#pre{{
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への入力:
#pre{{
# サポートベクターのポイントを強調して表示
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への入力:
#pre{{
(qp_plt + data_plt + sv2_plt).show()
}}
&ref(sage0-3.png);
** コメント [#l3c0f1ed]
#vote(おもしろかった,そうでもない,わかりずらい)
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha