[[FrontPage]]

#contents

2011/06/06からのアクセス回数 &counter;

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://sage.math.canterbury.ac.nz/home/pub/100/

また、Sageのサーバを公開しているサイト(http://sage.math.canterbury.ac.nz/ 、http://www.sagenb.org/)にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。


* Sageでスケール共役勾配法を試す [#lf27bb4c]

PRMLの第5章のニューラルネットワークの例題説明に良く出てくる「スケール共役勾配法」
Scaled Conjugate Gradientsという名前が出てきますが、Googleで検索しても
日本語のアルゴリズム説明が見当たらないので、オリジナルの論文を基にsageを使って
実装してみました。		

スケール共役勾配法は、ヘッセ行列を計算することなく、共役勾配を求めることができる
貴重な手法です。


** 論文 [#xaae54e8]

スケール共役勾配法の論文は、Moller, M (1993),
"A Scaled Conjugate Gradient Algorithm for Fast Supervised Learning"
に紹介されています。		

論文の終わりのアルゴリズム説明を書くと、

+ 重みベクトル\(w_1、スカラー \sigma > 0, \lambda_1 > 0, \bar{\lambda_1} = 0\)
を選択
+ success = trueなら、2次情報を計算
$$
\sigma_k = \frac{\sigma}{|p_k|},
$$
$$
s_k = \frac{E'(w_k + \sigma_k p_k) - E'(w_k)}{\sigma_k},
$$
$$
\delta_k = p_k^T s_k.
$$
+ スケール\(s_k\)は、
$$
s_k = s_k + (\lambda_k -\bar{\lambda_k}) p_k,
$$
$$
\delta_k = delta_k + (\lambda_k -\bar{\lambda_k}) |p_k|^2
$$
+ \(\delta_k \le 0\) なら、ヘッセマトリックスを正定値にする、
$$
s_k = s_k + (\lambda_k - 2 \frac{\delta_k}{|p_k|^2} p_k,
$$
$$
\bar{\lambda_k} = 2 (\lambda_k - \frac{\delta_k}{|p_k|^2}),
$$
$$
\delta_k = - \delta_k \lambda_k|p_k|^2, \lambda_k = \bar{\lambda_k}.
$$
+ 刻み幅を計算、
$$
\mu_k = p_k^T r_k, \alpha_k = \frac{\mu_k}{\delta_k}.
$$
+ 比較パラメータを計算、
$$
\Delta_k = \frac{2\delta_k[E(w_k) - E(w_k + \alpha_k p_k)}{\mu_k^w}.
$$
+ \(\Delta_k \ge 0\)なら、後退が正常に処理されており、
$$
w_{k+1} = w_k + \alpha_k p_k,
$$
$$
r_{k+1} = -E'(w_{k+1}),
$$
$$
\bar{\lambda_k} = 0, success = true.
$$
-- k mod N = 0なら、アルゴリズムをリスタート:
$$
p_{k+1} = r_{k+1}
$$
そうでなければ、
$$
\beta_k = \frac{|r_{k+1}|^2 - r_{k+1} r_k}{\mu_k},
$$
$$
p_{k+1} = r_{k+1} + \beta_k p_k
$$
-- \(\Delta_k \ge 0.75\) なら、スケールパラメータを減らす:
$$
\lambda_k = \frac{1}{2}\lambda_k
$$											
そうでなければ、後退は正ではなく、\(\bar{\lambda_k} = \frac{1}{2}\lambda_k, success = false\)とする
+ \(\Delta_k \lt 0.25\)なら、スケールを増やす:
$$
\lambda_k = 4 \lambda_k
$$
+ 勾配方向 \( r_k \ne 0 なら、  k = k + 1 \)とし、2に進む、
そうでなければ、終了し\(w_{k+1}\)を最小値として返す


*** テストデータの準備 [#r0d55e30]

スケール共役勾配法のテストは、図5.9

&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.9a.jpg,,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.9b.jpg,,30%);
&ref(http://research.microsoft.com/en-us/um/people/cmbishop/PRML/prmlfigs-jpg/Figure5.9c.jpg,,30%);

と同じsin曲線のデータを使用します。


sageへの入力:
#pre{{
# Scaled Conjugate Gradients(スケール共役勾配法)
# 論文 A Scaled Conjugate Gradient Algorithm for Fast Supervised Leaning
# テストデータ
data = matrix([
       [0.000000, 0.349486],
       [0.111111, 0.830839],
       [0.222222, 1.007332],
       [0.333333, 0.971507],
       [0.444444, 0.133066],
       [0.555556, 0.166823],
       [0.666667, -0.848307],
       [0.777778, -0.445686],
       [0.888889, -0.563567],
       [1.000000, 0.261502],
       ]);
X = data.column(0)
T = data.column(1)
}}


*** 設定 [#f99820a5]

ニューラルネットワークの設定は、逐次的勾配降下法とスケール共役勾配法の比較では

- 入力ユニット数1個(バイアス項を除く)
- 隠れ層ユニット数4個(バイアス項を除く)
- 出力ユニット数1個

図5.9の再現では隠れ層のユニット数をそれぞれ、1, 3, 10としました。

sageへの入力:
#pre{{
# 入力1点、隠れ層3個、出力1点
N_in = 1   
N_h = 4      
N_out = 1
# 隠れ層の活性化関数
h_1(x) = tanh(x)
# 出力層の活性化関数
h_2(x) = x
}}


** ニューラルネットワークのクラス化 [#qa8f46d8]

スケール共役勾配法では、何度もwの設定を変えたり、勾配を求めたりするため、
ニューラルネットワークをクラス化する方が便利です。

以下の様にニューラルネットワーク用のクラスを定義します。


sageへの入力:
#pre{{
class NeuralNetwork:
    """ ニューラルネットワーク用のクラス(バイアス項を含む) """
    def __init__(self, N_in, N_h, N_out, h_1, h_2, beta):
        # N_in : 入力ユニットの数(バイアス項を除く)
        # N_h  : 隠れ層の数(バイアス項を除く)
        # N_out: 出力ユニットの数
        # h_1  : 隠れ層の活性化関数
        # h_2  : 出力層の活性化関数
        # beta : wの初期化分布β
        self.N_in = N_in+1; self.N_h = N_h+1
        self.N_out = N_out; self.h_1 = h_1
        self.dif_h_1 = diff(h_1, x); self.h_2 = h_2
        self.w_1 = matrix(RDF, [[gauss(0,beta) for i in range(self.N_in)] for j in range(self.N_h)])
        self.w_2 = matrix(RDF, [[gauss(0,beta) for j in range(self.N_h)] for k in range(self.N_out)])
        self.x_i = vector(RDF, self.N_in); self.y_k = vector(RDF, self.N_out)
        self.z_j = vector(RDF, self.N_h); self.z_j[0] = 1
        self.d_j = vector(RDF, self.N_h); self.d_k = vector(RDF, self.N_out)
    # フィードフォワード
    def feedForward(self, x):
        # x : 入力値
        self.x_i = vector(RDF, flatten([1, x]))
        # フィードフォワード処理
        for j in range(1, self.N_h):           
            self.z_j[j] = self.h_1(self.w_1.row(j)*self.x_i).n()
        for k in range(self.N_out):
            self.y_k[k] = self.h_2(self.w_2.row(k)*self.z_j).n()
        return self.y_k
    # バックプロパゲート処理
    def backPropagate(self, t):
        # t : 観測値
        self.t = vector(RDF, flatten([t]))
        self.d_k = self.y_k - self.t
        for j in range(self.N_h):
            self.d_j[j] = self.dif_h_1(self.z_j[j])*self.w_2.column(j)*self.d_k
        return self.d_k
    # 1/2 En(w)を返す
    def En(self):
        return 1/2*self.d_k*self.d_k
    # ▽Enを返す
    def gradEn(self):
        dE_1 = [[self.d_j[j]*self.x_i[i] for i in range(self.N_in)] for j in range(self.N_h)]
        dE_2 = [[self.d_k[k]*self.z_j[j] for j in range(self.N_h)] for k in range(self.N_out)]
        return vector(RDF, flatten(dE_1 + dE_2))
    # wを返す
    def getW(self):
        return vector(RDF, self.w_1.list()+self.w_2.list())
    # wをセットする
    def setW(self, w):
       self.w_1 = matrix(RDF, self.N_h, self.N_in, w.list()[0:self.N_in*self.N_h])
       self.w_2 = matrix(RDF, self.N_out, self.N_h, w.list()[self.N_in*self.N_h:])
}}


** 逐次的勾配降下法による学習 [#nc7b5012]

比較のために、逐次的勾配降下法をニューラルネットワークのクラスを使って
実装します。


sageへの入力:
#pre{{
# 逐次勾配降下法による学習
def _learn_sg(net, X, T, eta, LEARN_COUNT):
    for m in range(LEARN_COUNT):
        for n in range(len(X)):
            net.feedForward(X[n])
            net.backPropagate(T[n])
            # print net.En()
            net.setW(net.getW() - eta*net.gradEn())
}}


** スケール共役勾配法による学習 [#h7a10d67]

スケール共役勾配法の学習アルゴリズムを実装します。


sageへの入力:
#pre{{
# Scaled Conjugate Gradients(スケール共役勾配法)
# 論文 A Scaled Conjugate Gradient Algorithm for Fast Supervised Leaning
def _learn_csg(net, X, T, LEARN_COUNT, RESET_COUNT):
    lam_min = 1e-15
    lam_max = 1e100
    w_len = len(net.getW())
    # ▽E(w), 1/2E(w)を返す
    def gradE(w):
        dE = vector(RDF, w_len)
        net.setW(w)
        E = 0
        for n in range(len(X)):
            net.feedForward(X[n])
            net.backPropagate(T[n])
            E += net.En()
            dE = dE + net.gradEn()
        return (dE, E)
    # 初期値設定
    sig0 = 1e-4
    w = net.getW()
    (dE, E) = gradE(w)
    o_E = E
    o_dE = dE
    p = r = -dE
    lam_bar = 0
    lam = 1
    success = true
    k = 1
    for k in range(1, LEARN_COUNT):
        if success :
            # 2次の計算が可能
            sig = sig0/sqrt(p*p)
            o_E = E; o_dE = dE
            (dE, E) = gradE(w + sig*p)
            s = (dE - o_dE)/sig
            d = p*s
        s = s + (lam - lam_bar)*p
        d = d + (lam - lam_bar)*(p*p)
        if d <= 0:
            s = s + (lam - 2*d/(p*p))*p
            lam_bar = 2*(lam - d /(p*p))
            d = -d + lam*(p*p)
            lam = lam_bar
        mu = p*r
        alpha = mu/d
        (dE, E) = gradE(w + alpha*p)
        delta = 2*d*(o_E - E)/(mu*mu)
        if delta >= 0:
            #print "update"
            o_E = E; o_dE = dE; o_r = r
            w = w + alpha*p
            (dE, E) = gradE(w)
            r = -dE
            lam_bar = 0; success = true
            k = k + 1
            if k % RESET_COUNT == 0:
                p = r
            else:
                beta = (r*r - r*o_r)/mu
                p = r + beta*p           
            if delta >= 0.75:
                lam = max(lam/2, lam_min)
            else:
                lam_bar =lam; success = false
        if delta < 0.25:
            lam = min(lam*4, lam_max)
        if r.norm() == 0:
            return
        #print k, delta, E
    net.setW(w)
}}


sageへの入力:
#pre{{
# 出力関数
def _f(x, net):
    y = net.feedForward(x)
    return y[0]
}}


sageへの入力:
#pre{{
# トレーニングデータのプロット
data_plt = list_plot(zip(X, T))
}}


** 逐次的勾配降下法の結果 [#m712691c]

比較のために、逐次的勾配降下法での学習結果をプロットします。

500回の学習回数では、sin曲線に近づけません。


sageへの入力:
#pre{{
# 比較のために、逐次勾配降下法で解く
# 定数設定
eta = 0.1    # η(学習率)
LEARN_COUNT = 500
RESET_COUNT = 50
beta = 0.01
# ニューラルネットワークの生成
net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# wの初期値を保存
saved_w = net.getW()
# 逐次勾配降下法による学習
_learn_sg(net, X, T, eta, LEARN_COUNT)
# 結果のプロット
var('x')
f_plt = plot(lambda x : _f(x, net), [x, 0, 1], color='red')
(data_plt + f_plt).show()
}}

&ref(sage0.png);


** スケール共役勾配法の結果 [#ub889c57]

スケール共役勾配法の学習結果をプロットします。

少ない学習回数でもかなり良くフィットしています。


sageへの入力:
#pre{{
# ニューラルネットワークの生成
net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# 逐次勾配降下法と同じwをセット
net.setW(saved_w)
# スケール共役勾配法による学習
_learn_csg(net, X, T, LEARN_COUNT, RESET_COUNT)
# 結果のプロット
var('x')
f_plt = plot(lambda x : _f(x, net), [x, 0, 1], color='red')
(data_plt + f_plt).show()
}}

&ref(sage0-1.png);


** Mを変えたスケール共役勾配法の結果 [#n25f49b7]

図5.9と同じようにMを1, 3, 10に変えた結果をプロットします。

M=10のような過学習の結果にはなりません、この辺は重みwの初期値の
分布を変えると変化することが分かりましたが、処理時間の関係で詳しくは調べていません。


sageへの入力:
#pre{{
# PRML fig. 5.9
# M=1
N_h = 1; net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# スケール共役勾配法による学習
_learn_csg(net, X, T, LEARN_COUNT, RESET_COUNT)
# 結果のプロット
f_plt = plot(lambda x : _f(x, net), [x, 0, 1], color='red')
(data_plt + f_plt).show()
}}

&ref(sage0-2.png);


sageへの入力:
#pre{{
# M=3
N_h = 3; net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# スケール共役勾配法による学習
_learn_csg(net, X, T, LEARN_COUNT, RESET_COUNT)
# 結果のプロット
f_plt = plot(lambda x : _f(x, net), [x, 0, 1], color='red')
(data_plt + f_plt).show()
}}

&ref(sage0-3.png);


sageへの入力:
#pre{{
# M=10
N_h = 10; net = NeuralNetwork(N_in, N_h, N_out, h_1, h_2, beta)
# スケール共役勾配法による学習
_learn_csg(net, X, T, LEARN_COUNT, RESET_COUNT)
# 結果のプロット
f_plt = plot(lambda x : _f(x, net), [x, 0, 1], color='red')
(data_plt + f_plt).show()
# 結果は、図5.9のM=10とは異なり、Over fittingは見えない
}}

&ref(sage0-4.png);



** コメント [#s79edfdd]
#vote(おもしろかった[4],そうでもない[0],わかりずらい[1])
#vote(おもしろかった[5],そうでもない[0],わかりずらい[1])

皆様のご意見、ご希望をお待ちしております。
- 自分でアップして何ですが、論文のステップ3の\( s_k = s_k + (\lambda_k -\bar{\lambda_k}) p_k,\)は、λのバー付きの位置がその前の説明と逆になっているのが、気になります。 -- [[竹本 浩]] &new{2011-06-07 (火) 01:48:42};

#comment_kcaptcha

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
SmartDoc