2010-03-06
OpenOpt使ってSVM書いた
Python |
追記(5/19):ガウスカーネル2乗してなかった。コード書き忘れ訂正--); ついでに画像も変更
SMO法使った前のエントリは、殆どpureにpythonでコード書いてたせいか、結構時間がかかっててイライラ。ということでOpenOptの二次計画のソルバー使って手抜きに疎な解を求めてみたの巻。
結果はテストデータ200個の↓の図だと200倍の差が…。scipy+OpenOptぱない
コーディングもあっと言う間だし…その…何というか…一昨日の努力は…一体…。まぁデータ200個と少なきゃメモリにのるしね…。
以下適当に書いたpythonのコード。相変わらずグラフの描画とかのコードの筋が悪い気がしてもにょいぜ。
#!/usr/bin/python # -*- coding: utf-8 -*- from scipy import * from scipy.linalg import norm from openopt import QP class svm(object): """ SVM """ def __init__(self, c=10000, kernel=lambda x,y:dot(x,y)): """ Constructor Arguments: - `c`: param - `kernel`: kernel func """ self._c = c self._kernel = kernel def _get_S(self): return self._S S = property(_get_S) def learn(self, x, t): """ Learning SVM (using openopt.QP) Arguments: - `x`: inputs - `t`: targets """ # making Gram matrix N = len(t) K = array([[t[i]*t[j]*self._kernel(x[i], x[j]) for j in range(N)] for i in range(N)]) p = QP(H=K, f=-ones(N), lb=zeros(N), ub=ones(N)*self._c, Aeq=t, beq=0) r = p.solve('nlp:ralg') self._a = r.xf self._S = [i for i in range(N) if 0 < self._a[i]] self._M = [i for i in range(N) if 0 < self._a[i] < self._c] b = 0.0 for i in self._M: b += t[i] for j in self._S: b -= self._a[j]*t[j]*self._kernel(x[i], x[j]) self._b = b/len(self._M) self._x = x self._t = t def calc(self, x): """ Arguments: - `x`: input """ return self._b + sum((self._a[i] * self._t[i] * self._kernel(x, self._x[i]) for i in self._S)) if __name__ == '__main__': from scipy.io import read_array import matplotlib.pyplot as plt import psyco psyco.full() s = svm(c=0.5, kernel=lambda x,y:exp(-square(norm(x-y))/0.45)) data = read_array(open("classification.txt")) p = data[:,0:2] t = data[:,2]*2-1.0 s.learn(p, t) for i in range(len(p)): c = 'r' if t[i] > 0 else 'b' plt.scatter([p[i][0]], [p[i][1]], color=c) X, Y = meshgrid(arange(-2.5, 2.5, 00.1), arange(-2.5, 2.5, 00.1)) w, h = X.shape X.resize(X.size) Y.resize(Y.size) Z = array([s.calc([x, y]) for (x, y) in zip(X, Y)]) X.resize((w, h)) Y.resize((w, h)) Z.resize((w, h)) CS = plt.contour(X, Y, Z, [0.0], colors = ('k'), linewidths = (3,), origin = 'lower') plt.xlim(-2.5, 2.5) plt.ylim(-2.5, 2.5) plt.show()
近況
雑記 |
日常
崩壊する生活リズム。今朝寝てさっき起きました。24h営業のスーパーがあると、思う存分生活崩せてやばい。
来週は高校の友人と東京に行く予定。向こうに行った友人に会ったり飲んだりするのかな。割と未定。
サークル
春合宿が近付いてきた為か、勉強会や読書会の勢いが収まってきた。まぁ皆さん講座の準備で忙しいんでしょう。って自分もだが。
最近機械学習のお勉強が多かったので、講座の内容はそっちの方向で話すことにする。がんばるよっ！
合宿が終われば新歓期、今年はゲーム作りのプロジェクトを持って、新入生と一緒に頑張りたい所。KMCをよろしくね！*1
エロゲ
星メモ本編+EH終了。本編は後半の章がちょっと短くて不満だったが、FDを合わせて見れば良い感じ。メアかわいいよメア。次はAiry[F]ariyを崩す予定。
ラノベ
この2週間で1冊も読んでない…、あれ？
その他
TopCoder参戦2回目。Div1で1問しか解けず1404→1416と微妙な結果。。。
撃墜フェーズでぼんやりしてるのを、そろそろどうにかしたい。
アニメの方はレールガンとプリキュアしか見なくなりました。放送されているほぼ全てのアニメを見ている、某アニメ廃人の友人には程遠い。いやなりたく無いけど。
SMO法でSVMの学習してみた
SMO法はPRMLに名前はあるけど詳細が無かったので、別の本*2買ってきて読んで理解した気分になったのでコードを書いてみた。
コードはpythonで書いたのだが、アルゴリズム的に行列の計算が全く無かった。
実際に試した対象は学習データが200程度だったので、scipyとOpenOptで素直に2次計画問題を計算させた方が早かった気がしてもにょい。*3
今はRVMの節を読んでいる。まだ7章終わってない＞ω＜；
#!/usr/bin/python # -*- coding: utf-8 -*- import sys import random from scipy import * from scipy.linalg import norm from scipy.io import read_array import matplotlib.pyplot as plt class svm(object): """ 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 j, lst) else: return reduce(lambda j,k: j if self._e[j] > self._e[k] else j, 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 - self._eps) or (r2 > self._tol and alph2 > self._eps)): alst1 = [i for i in range(len(self._alpha)) if self._eps < self._alpha[i] < self._c - self._eps] 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] <= self._eps or self._alpha[i] >= self._c - self._eps)] random.shuffle(alst2) for i1 in alst2: if self._takeStep(i1, i2): return True return False def _calc_b(self): 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 self._eps < self._alpha[j] < self._c - self._eps): changed |= self._examinEx(i) if examine_all: examine_all = False elif not changed: examine_all = True self._s = [i for i in range(len(self._target)) if self._eps < self._alpha[i]] self._m = [i for i in range(len(self._target)) if self._eps < self._alpha[i] < self._c - self._eps] self._calc_b() def _get_s(self): return self._s s = property(_get_s) def _get_m(self): return self._m m = property(_get_m) def _get_alpha(self): return self._alpha alpha = property(_get_alpha) if __name__ == '__main__': try: import psyco psyco.full() except ImportError: pass s = svm(c=2, kernel=lambda x,y:exp(-norm(x-y)/0.45)) data = read_array(open("classification.txt")) p = data[:,0:2] t = data[:,2]*2-1.0 s.learn(p, t) for i in range(len(p)): c = 'r' if t[i] > 0 else 'b' plt.scatter([p[i][0]], [p[i][1]], color=c) X, Y = meshgrid(arange(-2.5, 2.5, 00.1), arange(-2.5, 2.5, 00.1)) w, h = X.shape X.resize(X.size) Y.resize(Y.size) Z = array([s.calc([x, y]) for (x, y) in zip(X, Y)]) X.resize((w, h)) Y.resize((w, h)) Z.resize((w, h)) CS = plt.contour(X, Y, Z, [0.0], colors = ('k'), linewidths = (3,), origin = 'lower') plt.xlim(-2.5, 2.5) plt.ylim(-2.5, 2.5) plt.show() print s.alpha
takeStepメソッドのeta <= 0の処理をカットされた理由は、本文には説明がなく、必要性がないと判断されたからですか？
また、tolの値としてepsと同じ1e-2を使っていらっしゃいますが、javaで実装された方の例ですと、
http://d.hatena.ne.jp/nowokay/20080730/1217371769
1 - ε を使っています、この辺もどのような値にすればよいかお教え下さい。
最後に、「サポートベクターマシン入門」のアルゴリズムではbの更新をtakeStepメソッド内で行っていますが、learnメソッドに変更された理由をお教え下さい。
よろしくお願いします。
1. eta<=0 をカット: 覚えてないです．見逃してるだけかも……
2. メソッド名: 多分その時の気分です……．今ならそのままにすると思います．
そこで気づいた点をご報告します。
・_searchメソッドのelse jの部分はelse kの間違いでだと思います。
・kernelの定義で0.45で割ってありますが、図7.4の注釈では*0.45となっています。
（この場合、c=1.0でコンターとサポートベクタの分布が似た形になりました）
以上、貴重な情報をありがとうございました。