スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

GWに学んだことシリーズその1 サポートベクトル回帰

サポートベクトル回帰について学んだことを自分用にまとめる. $\bm x$を入力(説明変数),$\bm\phi(\bm x)$を特徴ベクトルとしたとき, 回帰に用いるモデルを \begin{align} \hat y(\bm x) = \bm w^\top \bm \phi(\bm x) + b \end{align} とする.また,N組の説明変数と目的変数 $\bm x_n, $,$y_n$($n = 1,\dots N$)が得られているとする.このとき,評価関数 \begin{align} C\sum^N_{n=1}E_\varepsilon\left(\hat y(\bm x_n)-y_n\right ) + \frac{1}{2}\|\bm w\|^2 \end{align} を最小化するパラメータ$\bm w$,$b$を求めるのがサポートベクトル回帰である.ただし, \begin{align} E_\varepsilon(e) = \begin{cases} 0 & (|e| < \varepsilon) \\ |e| - \varepsilon & (\mathrm{otherwise}) \end{cases} \end{align} である.この最適化問題をいろいろ頑張って書き換えていくと,カーネル関数 \begin{align} k(\bm x, \bm x') = \bm \phi^\top(\bm x)\bm \phi(\bm x') \end{align}を用いてモデルを \begin{align} \hat y(\bm x) = -\sum^N_{n=1}\left(a_n^+ - a_n^-\right)k(\bm x,\bm x_n) + b \end{align}のように書き換えることができて,このパラメータ推定問題は, \begin{align} 0 \leq a_n^+ \leq C\\ 0 \leq a_n^- \leq C \end{align}のもとで \begin{align} L &= \bm a^\top \bm P \bm a + \bm a^\top \bm q\\ \end{align} を最小化する二次計画問題に落ちる. ただし, \begin{align} \bm a &= \begin{bmatrix} a^+_1 & \dots & a^+_N & a_1^- & \dots & a_N^- \end{bmatrix}^\top\\ \bm P &= \begin{bmatrix}\bm K & -\bm K\\ -\bm K& \bm K\end{bmatrix}\\ \bm q &= \begin{bmatrix} y_1 + \varepsilon & \dots & y_N + \varepsilon & -y_1 + \varepsilon & \dots & -y_N + \varepsilon \end{bmatrix}^\top \end{align} である.ここで,$\bm K$はグラム行列と呼ばれて,$(n,m)$成分が$k(\bm x_n, \bm x_m)$で与えられる$N \times N$行列である.
この最適化問題を解くと多くの$n$について$a_n^+ = a_n^- = 0$となり,スパースなモデルが得られる.
このとき, $a_n^+ \neq 0$または$a_n^- \neq 0$が成り立つサンプル$\bm x_n$をサポートベクトルと呼ぶ.

最後に,$b$は$a_n^+ \neq 0$または$a_n^- \neq 0$が成り立つ$n$について \begin{align} \begin{cases} y_n + \varepsilon + \sum^N_{m=1}\left(a_m^+ - a_m^-\right)k(\bm x_n,\bm x_m) & (a_n^+\neq 0)\\ y_n - \varepsilon + \sum^N_{m=1}\left(a_m^+ - a_m^-\right)k(\bm x_n,\bm x_m) & (a_n^-\neq 0) \end{cases} \end{align}を計算して,平均すれば求められる.


以上をpython(numpy, cvxopt, matplotlib)を用いて実装したのが以下のコード.
import numpy as np
import matplotlib.pylab as plt
import cvxopt
import cvxopt.solvers
class SVR:
    sv = None
    coeff = None
    b = None
    def __init__(self, kernel, e, C):
        self.kernel = kernel
        self.epsilon = e
        self.C = C
    def train(self, x, y):
        N = x.shape[0]
        K = self._makeKmat(x)
        P0 = np.vstack([np.hstack([K, -K]), np.hstack([-K, K])])
        q0 = np.hstack([y, -y]) + self.epsilon
        G0 = np.vstack([np.eye(2 * N), -np.eye(2 * N)])
        h0 = np.hstack([np.ones(2 * N) * self.C, np.zeros(2 * N)])
        A0 = np.hstack([np.ones(N), -np.ones(N)])[np.newaxis, :]
        b0 = 0.
        P = cvxopt.matrix(P0)
        A = cvxopt.matrix(A0)
        q = cvxopt.matrix(q0)
        G = cvxopt.matrix(G0)
        h = cvxopt.matrix(h0)
        b = cvxopt.matrix(b0)
        sol = cvxopt.solvers.qp(P, q, G, h, A, b)
        a = np.array(sol['x']).flatten()
        a[a < 1e-5] = 0
        a_tilde = -a[0:N] + a[N:2*N]
        self.sv = x[a_tilde != 0, :]
        self.coeff = a_tilde[a_tilde != 0]
        self.b = 0
        yhat = np.dot(K, a_tilde)
        bhat = y[a[N:2*N]!=0] - self.epsilon - yhat[a[N:2*N]!=0]
        bhat2 = y[a[0:N]!=0] + self.epsilon - yhat[a[0:N]!=0]
        self.b = np.mean(np.hstack([bhat, bhat2]))
        print "number of support vector is", self.coeff.shape[0]
    def _makeKmat(self, x):
        N = x.shape[0]
        K = np.zeros((N, N))
        for i in range(N):
            K[:,i] = self.kernel.get_value(x, x[i, :])
        return K
    def get_value(self, x):
        N = x.shape[0]
        k = np.zeros((N, self.coeff.shape[0]))
        for i in range(N):
            k[i,:] = self.kernel.get_value(x[i, :], self.sv)
        yhat = np.dot(k, self.coeff) + self.b
        return yhat
    
class gaussianKernel:
    def __init__(self, sigma = None):
        if not sigma:
            self.sigma = 1
        else:
            self. sigma = sigma
    def get_value(self, x1, x2):
       x =  x1 - x2
       if x.ndim == 1:
           e2 = np. inner(x, x)
       else:
           e2 = np.sum(x * x, 1)
       val = np.exp(-e2/(2 * self.sigma**2))
       return val

if __name__ == '__main__':
    kernel = gaussianKernel()
    svr = SVR(kernel, 0.1, 100)
    N = 200
    X = np.random.normal(size = (N, 2))
    y = np.sin(X[:,0]) + np.sin(2 * X[:, 1]) + 5 
    svr.train(X, y + 1e-3 * np.random.normal(size = N))
    yhat = svr.get_value(X)
    plt.plot(y)
    plt.plot(yhat)
    plt.show()
スポンサーサイト
プロフィール

kawa0616

Author:kawa0616
某企業で研究開発してます.専門は制御工学と機械学習を少々.

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
検索フォーム
RSSリンクの表示
リンク
QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。