FrontPage

2011/06/28からのアクセス回数 3011

ここで紹介した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、

Figure7.4.jpg

で使用された(赤と青の点が入り混ざっている)人工データを使います。

カーネルには、\(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()

sage0.png

SMO

「きちめも」で実装されている、「サポートベクターマシン入門」の付録に掲載されているSMO(Sequenctial Minimal Optimization) は、J. Platt(1988)で発表されたもので、2つのデータ点に対する最適化問題を解析的に求め、ヒューリスティクスな選択で反復的な処理を 高速に収束させる手法です。

以下では、「きちめも」の実装に若干の修正を加えています。

  • ガウスカーネルを図7.4の説明に合わせた
  • 「サポートベクターマシン入門」の注釈7.14からcalc_bを_examinExの最後にも追加
  • KKTのチェックは、「サポートベクターマシン入門」の付録Aに合わせた。

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()

sage0-1.png

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()

sage0-2.png

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} $$ とすることで解決されたそうです。

「人工知能に関する断想録」からの変更点は

  • 「人工知能に関する断想録」のコメントからサポートベクタの抽出条件を \( \epsilon \le a_n\)に変更しました。

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()

sage0-3.png

コメント

選択肢 投票
おもしろかった 4  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。

  • 最後のプログラムで、bを計算した後、下記の最小チェックを行ってみたのですが、この計算が悪いせいか最小を示せなくて悩んでいます. print np.sum([max(0, 1 - y * f(x, a, t, X, b)) for (x, y) in zip(X, t)]) <br/> -- zyx? 2012-11-13 (火) 20:01:26
  • 確認していませんがf(x, a, t, X, b)は、f([x, y], a, t, X, b)の間違いではないかと思いますが、いかがでしょう。 -- 竹本 浩? 2012-11-14 (水) 12:49:11
  • ありがとうございます。b の計算は、if eps < a[n] < C - eps: M.append(n) などと書き換えて計算すれば、np.sum([max(0, 1 - y * f(x, a, t, X, b)) for (x, y) in zip(X, t)]) も最小になり、この M に属するサポートベクトルも f(x, a, t, X, b) = ±1 の等高線上に乗りました。 -- zyx? 2012-11-14 (水) 19:53:34
  • 結局、そのままでは収束不足で、cvxopt.solvers で options['abstol'] = 1e-5 options['reltol'] = 1e-10 などと指定して計算する必要があるようです。サポートベクトルかどうかの判定や、±1の等高線も正しく計算できるようになりました。 -- zyx? 2013-06-16 (日) 22:05:47
  • zyxさま、貴重な情報ありがとうございます。 -- 竹本 浩? 2013-06-17 (月) 06:27:41

(Input image string)


添付ファイル: filesage0-3.png 454件 [詳細] filesage0-2.png 458件 [詳細] filesage0-1.png 475件 [詳細] filesage0.png 433件 [詳細]

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-10-10 (金) 15:26:11 (1108d)
SmartDoc