ラベル Chainer の投稿を表示しています。 すべての投稿を表示
ラベル Chainer の投稿を表示しています。 すべての投稿を表示

2018年7月1日日曜日

Adversarial Variational Bayes 〜理論と実装〜

はじめに


 論文Adversarial Variational Bayesの前半部分の解説とChainerによる実装を行う。

論文概要


 深層学習において生成モデルを求める代表的な手法は、Generative Adversarial Network(GAN)とVariational AutoEncoder(VAE)である。最初にGANが最適化する目的関数を示す。 \begin{equation} \min_{\theta}\max_{\phi} \left\{ {\rm E}_{p_{d}(x)}\left[\log{D_{\phi}(x)}\right]+{\rm E}_{p(z)}\left[1-\log{D_{\phi}\left(G_{\theta}(z)\right)}\right] \right\} \end{equation} ここで、$p_d(x)$と$p(z)$はそれぞれ観測値$x$と潜在変数$z$の分布を表す。 ${\rm E}_{p}\left[\cdot\right]$は分布$p$による期待値である。$D_{\phi}$は識別器、$G_{\theta}$は生成器であり、それぞれパラメータ$\phi$と$\theta$で象徴される重みを持つネットワークである。上の式は次の2つの最適化を行うことを表す。
  1. $\theta$を固定して考える。観測値$x$と生成器$G_{\theta}(z)$の出力値を識別器$D_{\phi}$に与えたとき、$D_{\phi}(x)$の値を最大にし、$D_{\phi}(G_{\theta}(z))$の値を最小にするような$\phi$を見つける。
  2. $\phi$を固定して考える。観測値$x$と生成器$G_{\theta}(z)$の出力値を識別器$D_{\phi}$に与えたとき、$D_{\phi}(x)$の値を最小にし、$D_{\phi}(G_{\theta}(z))$の値を最大にするような$\theta$を見つける。
訓練の後、生成モデル$G_{\theta}(z)$を得ることになる。一方、VAEが最適化すべき式は次式である。 \begin{equation} \max_{\theta} \max_{\phi} \left\{ -{\rm KL} \left( q_{\phi}(z|x), p(z) \right) + {\rm E}_{q_{\phi}(z|x)} \left[\ln{p_{\theta}(x|z)}\right] \right\} \end{equation} ${\rm KL}$はKullback Leibler divergence、$q_{\phi}(z|x)$は観測値$x$が与えられた下での潜在変数$z$の事後確率を近似する分布である。上式第1項は正則化項、第2項は尤度を表す。すなわち、正則化項を考慮して尤度を最大化することになる。訓練の後、$q_{\phi}(z|x)$が推論モデル(inference model)、$p_{\theta}(x|z)$が生成モデル(generative model)になる。VAEの処理の流れを以下に示す。
GANとVAEを比べると、後者からは生成モデル以外に推論モデルも求めることができる。また、GANで決まる生成器は決定論的な関数であるが、VAEから求まる生成器は確率モデルである。さらに、GANには尤度に相当する量が存在せず、VAEには識別器に相当する仕組みが存在しない。

 GANを自然画像に適用すると、シャープな画像を生み出すが、VAEではボケた画像になることが知られている。これは、VAEの推論モデル$q_{\phi}(z|x)$が真の事後確率$p(z|x)$を再現できていないためである。今回紹介する文献は、任意の形状を持つ推論モデルを扱うことのできる手法(Adversarial Variational Bayes:AVB)を提案している。本文献の主張をまとめると次のようになる。
  1. 任意の複雑な推論モデルを扱うことができる。
  2. VAEの目的関数に数学的な変換を施すことによりGAN的な目的関数を導出する。すなわち、VAEに識別器に相当する量が導入される。
  3. 上の目的関数は、ある極限の下で、VAEの元の目的関数と厳密に一致する。

論文詳細


 観測値$x$の分布$p(x)$を考え、潜在変数$z$を導入する。 \begin{eqnarray} \ln{p(x)} &=&\ln\int dz\;p(x|z)p(z) \\ &=&\ln\int dz\;q(z|x) \frac{p(x|z)p(z)}{q(z|x)} \end{eqnarray} Jensenの不等式を用いて \begin{eqnarray} \ln{p(x)} &\geqq&\int dz\;q(z|x) \ln\frac{p(x|z)p(z)}{q(z|x)} \\ &=&\int dz\;q(z|x)\ln{p(x|z)}-\int dz\;q(z|x)\ln{\frac{q(z|x)}{p(z)}} \\ &=&{\rm E}_{q(z|x)}\left[\ln{p(x|z)}\right]-{\rm KL}\left[q(z|x),p(z)\right] \\ &\equiv&\mathcal{L} \end{eqnarray} を得る。上式右辺はEvidence Lower Bound(ELOB)と呼ばれる量である。ここで、次の量を考える。 \begin{eqnarray} {\rm KL}\left(q(z|x),p(z|x)\right) &=&\int dz\;q(z|x)\ln{\frac{q(z|x)}{p(z|x)}} \\ &=&\int dz\;q(z|x)\ln{\frac{p(x)q(z|x)}{p(x,z)}} \\ &=&\int dz\;q(z|x)\ln{\frac{p(x)q(z|x)}{p(x|z)p(z)}} \\ &=&\int dz\;q(z|x)\ln{p(x)}-\int dz\;q(z|x)\ln{p(x|z)}+{\rm KL}(q(z|x),p(z)) \\ &=&\ln{p(x)}-\int dz\;q(z|x)\ln{p(x|z)}+{\rm KL}(q(z|x),p(z)) \\ \end{eqnarray} 上式は先に定義したELBO($\mathcal{L}$)を用いると \begin{equation} {\rm KL}\left(q(z|x),p(z|x)\right)=\ln p(x)-\mathcal{L} \end{equation} と書くことができる。すなわち次式が成り立つ。 \begin{equation} \ln p(x)=\mathcal{L}+{\rm KL}\left(q(z|x),p(z|x)\right) \end{equation} $q(z|x)$は分布$p(z|x)$を近似するために導入された分布であり、パラメータ$\phi$を持つとする。一方、分布$p(x|z)$を表現するモデルはパラメータ$\theta$を持つとする。これらパラメータを顕に書くと \begin{eqnarray} \ln p(x)&=&\mathcal{L}+{\rm KL}\left(q_{\phi}(z|x),p(z|x)\right) \label{eq2}\\ \ln p(x)&\geqq&\mathcal{L} \label{eq1}\\ \mathcal{L}&=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-{\rm KL}\left[q_{\phi}(z|x),p(z)\right] \end{eqnarray} となる。VAEでは式(\ref{eq1})の右辺を$\phi$と$\theta$について最大化する。式(\ref{eq2})から、$\ln p(x)=\mathcal{L}$が成り立つのは$q_{\phi}(z|x)=p(z|x)$となる時であることが分かるが、一般にこの等式が成り立つことはない。通常のVAEでは実際に計算を進める際、$q_{\phi}(z|x)$を正規分布で近似する。本文献では$q_{\phi}(z|x)$に対してそのような近似を行わない。

 $\mathcal{L}$は以下のように書き換えることができる。 \begin{eqnarray} \mathcal{L} &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-{\rm KL}\left[q_{\phi}(z|x),p(z)\right] \\ &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-\int dz\;q_{\phi}(z|x)\ln{\frac{q_{\phi}(z|x)}{p(z)}} \\ &=&{\rm E}_{q_{\phi}(z|x)}\left[\ln{p_{\theta}(x|z)}\right]-\int dz\;q_{\phi}(z|x) \left[ \ln{q_{\phi}(z|x)}-\ln{p(z)} \right] \\ &=&{\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-\ln{q_{\phi}(z|x)}+\ln{p(z)} \right]\\ \end{eqnarray} 従って、VAEの目的関数は \begin{equation} \max_{\theta}\max_{\phi} {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-\ln{q_{\phi}(z|x)}+\ln{p(z)} \right] \label{eq4} \end{equation} となる。ここまでは通常のVAEである。本文献では、ここで、関数$T^*(x,z)$を次式で導入する。 \begin{equation} T^*(x,z)=\arg\max_{T} \left\{ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{\sigma\left(T(x, z)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq3} \end{equation} この式は、$q_{\phi}(z|x)$からサンプリングされた$z$のとき$T(x,z)$を大きくし、$p(z)$からサンプリングされた$z$のとき$T(x,z)$を小さくすることを意味する。すなわち、$T(x,z)$は$q_{\phi}(x|z)$と$p(z)$を識別する識別器である。実際に$T^*(x,z)$を求めるため、式(\ref{eq3})の期待値を積分に書き換える。 \begin{equation} \int dz \left[ q_{\phi}(z|x) \ln{\sigma\left(T(x, z)\right)} + p(z) \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \end{equation} 上式の被積分関数の最大値を求めるため、$a=q_{\phi}(z|x)$、$b=p(z)$、$t=\sigma\left(T(x,z)\right)$として次式を考える。 \begin{equation} y=a\ln{t}+b\ln{(1-t)} \end{equation} 両辺$t$で微分して \begin{equation} \frac{dy}{dt}=\frac{a}{t}-\frac{b}{1-t} \end{equation} 右辺を0とおいて \begin{equation} t=\frac{a}{a+b} \end{equation} を得る。各変数を元に戻して計算すると \begin{equation} T(x,z)=\ln{q_{\phi}(z|x)}-\ln{p(z)}\equiv T^{*}(x,z) \end{equation} を得る。$T^{*}$を用いると、式(\ref{eq4})は以下のように書くことができる。 \begin{equation} \max_{\theta}\max_{\phi} {\rm E}_{q_{\phi}(z|x)} \left[ \ln{p_{\theta}(x|z)}-T^{*}(x,z) \right] \label{eq5} \end{equation} ただし \begin{equation} T^{*}(x,z)=\arg{\max_{T}} \left\{ {\rm E}_{q_{\phi}(z|x)} \left[ \ln{\sigma\left(T(x, z)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq6} \end{equation} である。式(\ref{eq5})と(\ref{eq6})に再パレメータ化トリックを用いると \begin{eqnarray} T^{*}(x,z)&=&\arg{\max_{T}} \left\{ {\rm E}_{p(\epsilon)} \left[ \ln{\sigma\left(T(x, z_{\phi}(x,\epsilon)\right)} \right] + {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \\ (\theta^*,\phi^*)&=&\arg\max_{\theta,\phi} {\rm E}_{p(\epsilon)} \left[ \ln{p_{\theta}\left(x|z_{\phi}\left(x,\epsilon\right)\right)}-T^{*}\left(x,z_{\phi}\left(x,\epsilon\right)\right) \right] \label{eq7} \end{eqnarray} を得る。$p(\epsilon)$は標準正規分布とすれば良い。観測値の分布$p_d(x)$による期待値も考慮して最終的に次の2つの目的関数を考えることになる。 \begin{eqnarray} T^{*}(x,z)&=&\arg{\max_{T}} \left\{ {\rm E}_{p_d(x)} {\rm E}_{p(\epsilon)} \left[ \ln{\sigma\left(T(x, z_{\phi}(x,\epsilon)\right)} \right] + {\rm E}_{p_d(x)} {\rm E}_{p(z)} \left[ \ln{\left(1-\sigma\left(T(x, z)\right)\right)} \right] \right\} \label{eq9}\\ (\theta^*,\phi^*)&=&\arg\max_{\theta,\phi} {\rm E}_{p_d(x)}{\rm E}_{p(\epsilon)} \left[ \ln{p_{\theta}\left(x|z_{\phi}\left(x,\epsilon\right)\right)}-T^{*}\left(x,z_{\phi}\left(x,\epsilon\right)\right) \right] \label{eq8} \end{eqnarray} 上の2式が本文献の最初の手法である。$T$を最適化したあと$(\theta,\phi)$を最適化する。アルゴリズムが文献に掲載されている。
ニューラルネットワークで表現される量は次の3つである。
  1. $z_{\phi}(x,\epsilon)$
  2. $p_{\theta}(x|z)$のパラメータ
  3. $T_{\psi}(x,z)$
文献の図5は、本アルゴリズムを評価するために用いられた実験データである。4次元データであり1つの成分だけが1、残りが0のデータである。
図6は$z_{\phi}(z|x)$から求まる潜在変数$z$の分布である。実験データは4種類のベクトルからなる。それぞれのデータごとに対応する$z$が色分けされている。同じデータに対する従来のVAEの結果も示されている。AVBの方が密に分布していることが分かる。 VAEの方はそれぞれの分布がガウス分布しているように見える。
下の表1は、各種評価値である。
表中のlog-likelihoodは$\ln{p_{\theta}(x|z)}$を、reconstruction errorは観測値$x$と再構築した$x$の間のクロスエントロピーを表す。いずれの値も従来法より精度が上がっていることが分かる。ここまでが本文献の前半の内容である。

Chainerによる実装


 ここからは、Chainerを用いて実際に実装を行い、上の図5に示された実験データから図6(b)の結果を再現するまでの様子を記載する。ソースコードはここにある。既存の実装として参考にしたのは以下の3つである。
  1. https://gist.github.com/poolio/b71eb943d6537d01f46e7b20e9225149
  2. https://github.com/gdikov/adversarial-variational-bayes
  3. https://github.com/LMescheder/AdversarialVariationalBayes
1はtensorflow、2はkeras、3は本文献の著者によるtensorflowによる実装である。

Encoder


 まず最初にEncoder($z_{\phi}(x,\epsilon)$)の実装を示す(encoder.py内のコードである)。
class Encoder_2(chainer.Chain):

    def __init__(self, x_dim, eps_dim, h_dim=512):
        super(Encoder_2, self).__init__()
        with self.init_scope():
            self.l1 = L.Linear(x_dim + eps_dim, h_dim, initialW=xavier.Xavier(x_dim + eps_dim, h_dim))
            self.l2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l4 = L.Linear(h_dim, eps_dim, initialW=xavier.Xavier(h_dim, eps_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, xs, es, activation=F.relu):
        xs = 2 * xs - 1 
        h = F.concat((xs, es), axis=1)

        h = self.l1(h)
        h = activation(h)

        h = self.l2(h)
        h = activation(h)

        h = self.l3(h)
        h = activation(h)

        h = self.l4(h)
        return h
  1. 入力値は観測値xsと標準正規分布からのサンプル値esである。
  2. 15行目:xsは0と1である。これを-1と1に置き換える。
  3. 16行目:xsとesを連結する。
  4. あとは、全結合層と活性化関数の繰り返しである。
  5. 出力層には活性化関数を適用しない。
  6. 活性化関数としてreluを採用した。

Decoder


 次はDecoderである(decoder.py内のコードである)。実験データは0と1から構成されるので$p_{\theta}(x|z)$としてBernoulli分布を用いる。 \begin{eqnarray} p_{\theta}(x|z) &=&{\rm Bern}(x|\mu(z))\\ &=&\mu(z)^x(1-\mu(z))^{1-x} \end{eqnarray} Decoderはパラメータ$\mu(z)$を計算する。
class Decoder_1(chainer.Chain):

    def __init__(self, z_dim, x_dim=1, h_dim=512):
        super(Decoder_1, self).__init__()
        with self.init_scope():
            self.l1 = L.Linear(z_dim, h_dim, initialW=xavier.Xavier(z_dim, h_dim))
            self.l2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.l4 = L.Linear(h_dim, x_dim, initialW=xavier.Xavier(h_dim, x_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, zs, activation=F.tanh, is_sigmoid=False):
        h = self.l1(zs)
        h = activation(h)

        h = self.l2(h)
        h = activation(h)

        h = self.l3(h)
        h = activation(h)

        h = self.l4(h)
        if is_sigmoid:
            h = F.sigmoid(h)
        return h
  1. 入力値は潜在変数zsである。
  2. 最終層にだけ仕掛けがしてある。訓練時はis_sigmoid=False、テスト時にはis_sigmoid=Trueとする。

Discriminator


 次はDiscriminator($T_{\psi}(x,z)$)である(discriminator.py)。
class Discriminator_1(chainer.Chain):

    def __init__(self, x_dim, z_dim, h_dim=512):
        super(Discriminator_1, self).__init__()
        self.h_dim = h_dim
        with self.init_scope():
            self.xl1 = L.Linear(x_dim, h_dim, initialW=xavier.Xavier(x_dim, h_dim))
            self.xl2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.xl3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.zl1 = L.Linear(z_dim, h_dim, initialW=xavier.Xavier(z_dim, h_dim))
            self.zl2 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))
            self.zl3 = L.Linear(h_dim, h_dim, initialW=xavier.Xavier(h_dim, h_dim))

    def update(self, updates):
        update_links(self, updates)

    def __call__(self, xs, zs, activation=F.relu):
        xs = 2 * xs - 1 
        hx = self.xl1(xs)
        hx = activation(hx)
        hx = self.xl2(hx)
        hx = activation(hx)
        hx = self.xl3(hx)
        hx = activation(hx)

        hz = self.zl1(zs)
        hz = activation(hz)
        hz = self.zl2(hz)
        hz = activation(hz)
        hz = self.zl3(hz)
        hz = activation(hz)
        h = F.sum(hx * hz, axis=1) / self.h_dim
        return h
  1. 入力値は観測値xsと潜在変数zsである。
  2. 18行目:xsは0と1である。これを-1と1に置き換える。
  3. 19行目から24行目:xsについての計算である。
  4. 26行目から31行目:zsについての計算である。
  5. 32行目:それぞれの結果の内積を取る。

$(\theta,\phi)$についての目的関数(式(\ref{eq8}))


 式(\ref{eq8})を実装したものが以下である(phi_loss_calculator.py)。
class PhiLossCalculator_2(chainer.Chain):
    
    def __init__(self, encoder, decoder, discriminator):
        super(PhiLossCalculator_2, self).__init__()
        with self.init_scope():
            self.encoder = encoder
            self.decoder = decoder
            self.discriminator = discriminator

    def __call__(self, xs, zs, es):
        batch_size = xs.shape[0]
        encoded_zs = self.encoder(xs, es)
        ys = self.decoder(encoded_zs)
        d_loss = F.bernoulli_nll(xs, ys) / batch_size
        t_loss = F.sum(self.discriminator(xs, encoded_zs)) / batch_size
        return t_loss + d_loss, encoded_zs
  1. 式(\ref{eq8})では最大化しているが、実装では符号を反転させて最小化させる。

$\psi$についての目的関数(式(\ref{eq9}))


 式(\ref{eq9})を実装したものが以下である(psi_loss_calculator.py)。
class PsiLossCalculator_3(chainer.Chain):

    def __init__(self, encoder, discriminator):
        super(PsiLossCalculator_3, self).__init__()
        with self.init_scope():
            self.encoder = encoder
            self.discriminator = discriminator

    def __call__(self, xs, zs, es):
        # batch_size = xs.shape[0]
        encoded_zs = self.encoder(xs, es) 
        posterior = self.discriminator(xs, encoded_zs)
        prior = self.discriminator(xs, zs) 
        a = F.sigmoid_cross_entropy(posterior, np.ones_like(posterior).astype(np.int32))
        b = F.sigmoid_cross_entropy(prior, np.zeros_like(prior).astype(np.int32))
        c = F.sum(a + b)
        return c
  1. こちらの目的関数も符号を反転させて最小化させる。

訓練コード


 訓練コードの一部を掲載する。
    for epoch in range(args.epochs):
        with chainer.using_config('train', True):
            # shuffle dataset
            sampler.shuffle_xs()

            epoch_phi_loss = 0 
            epoch_psi_loss = 0 

            for i in range(batches):
                xs = sampler.sample_xs()
                zs = sampler.sample_zs()
                es = sampler.sample_es()

                # compute psi-gradient(eq.3.3)
                update_switch.update_models(enc_updates=False, dec_updates=False,
                                            dis_updates=True)
                psi_loss = psi_loss_calculator(xs, zs, es) 
                update(psi_loss, psi_loss_calculator, psi_optimizer)
                epoch_psi_loss += psi_loss

                # compute phi-gradient(eq.3.7)
                update_switch.update_models(enc_updates=True, dec_updates=True,
                                            dis_updates=False)
                phi_loss, _ = phi_loss_calculator(xs, zs, es) 
                update(phi_loss, phi_loss_calculator, phi_optimizer)
                epoch_phi_loss += phi_loss

            # end for ...
            # see loss per epoch
            epoch_phi_loss /= batches
            epoch_psi_loss /= batches

        # end with ...
        print('epoch:{}, phi_loss:{}, psi_loss:{}'.format(epoch, epoch_phi_loss.data,
                                                          epoch_psi_loss.data))
        epoch_phi_losses.append(epoch_phi_loss.data)
        epoch_psi_losses.append(epoch_psi_loss.data)

    # end for ...
  1. 文献に掲載されているアルゴリズム図をほぼ踏襲している。
  2. 15行目から19行目:目的関数(\ref{eq9})を更新する。その際、EncoderとDecoderの重みを固定する(15行目)。
  3. 22行目から26行目:目的関数(\ref{eq8})を更新する。その際、Discriminatorの重みを固定する(22行目)。

予測のためのコード


 訓練後に実行するコートはpredict.pyである。掲載は略。

実行


 次を実行する。
$> ./train.py
$> ./predict.py


結果


 結果の表示はvisualize.ipynbで行った。最初に式(\ref{eq9})と(\ref{eq8})の変化を示す。
 以下が学習後の$z_{\phi}(x,\epsilon)$の結果である。
Lossの振る舞いは良くわからないが、潜在変数$z$の分布は$(0,0)$を中心に綺麗に4分割されている。文献掲載の結果をそれなりに再現できているように見えるが。。。

まとめ


 論文Adversarial Variational Bayesの前半部分の解説と実装を示した。正直なところ、上のような形に到達するまでかなりの時間を費やして試行錯誤を繰り返した。2つの目的関数を交互に最適化するのは大変難しい。ご批判いただければ幸いである。

2018年5月2日水曜日

Variational Auto Encoder 〜その3〜

はじめに


 先のページで、Variational Auto Encoder(VAE)を実装したChainerのサンプルコードをそのまま動かし、生成される画像を見た。符号化・復号化して得られる画像は、入力画像をそれなりに再現していたが、乱数から生成される画像は精度が良くないことを示した。今回は、サンプルコードに手を加えて実験を行う。

ソースコードの場所


 今回のソースコードはここにある。

評価基準


 このページで示したように、勾配降下法で最適化すべき式は次式である。 \begin{equation} \min_{\phi} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = \min_{\phi} {\left[ D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right] \right] } \end{equation} この右辺第1項は次のように書けた。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \label{eq1} \end{equation} また、潜在変数$\vec{z}$の成分は次式で与えられた。 \begin{equation} z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d \end{equation} 式(\ref{eq1})で理想的な最適化を実現できれば、$\mu_{\phi,d}(X)\rightarrow 0,\sigma_{\phi,d}(X)\rightarrow 1$とできる。このとき、$\vec{z}$は標準正規分布から生成される値$\vec{\epsilon}$に置き換えることができる。net.py内にある関数decodeは本来ならば$\vec{z}$を与えて画像を生成するコードであるが、実際に使われる際は標準正規分布から作られる値を与えている。従って、上記の収束が不十分であると、意図した振る舞いをしないことになる。最適化を行う際に、$\mu_{\phi,d}(X)$と$\sigma_{\phi,d}(X)$の値の変化も追跡すべきである。以上の考察に基づきnet.pyに実装されている関数get_loss_funcを以下のように変更した。
    def get_loss_func(self, C=1.0, k=1):
        """Get loss function of VAE.
        The loss value is equal to ELBO (Evidence Lower Bound)
        multiplied by -1.
        Args:
            C (int): Usually this is 1.0. Can be changed to control the
                second term of ELBO bound, which works as regularization.
            k (int): Number of Monte Carlo samples used in encoded vector.
        """
        def lf(x):
            mu, ln_var = self.encode(x)
            mean_mu, mean_sigma = calculate_means(mu, ln_var)
            batchsize = len(mu.data)
            # reconstruction loss
            rec_loss = 0
            for l in six.moves.range(k):
                z = F.gaussian(mu, ln_var)
                rec_loss += F.bernoulli_nll(x, self.decode(z, sigmoid=False)) \
                    / (k * batchsize)
            self.rec_loss = rec_loss
            kl = gaussian_kl_divergence(mu, ln_var) / batchsize
            self.loss = self.rec_loss + C * kl
            chainer.report(
                {
                    'rec_loss': rec_loss,
                    'loss': self.loss,
                    'kl': kl,
                    'mu': mean_mu,
                    'sigma': mean_sigma,
                },
                observer=self)
            return self.loss
        return lf
12行目でmuとsigmaの平均値を計算している。これに伴い、実行中にコマンドラインに表示する項目を増やしている( 25行目から29行目)。関数calculate_meansの中身は以下の通りである。
def calculate_means(mu, ln_var):
    xp = chainer.cuda.get_array_module(mu)
    mean_mu = xp.mean(mu.data)
    sigma = xp.exp(ln_var.data / 2)
    mean_sigma = xp.mean(sigma)
    return mean_mu, mean_sigma

実験-1


 先のページと同じパラメータ(epoch=100,dimz=20)で実行したときのmuとsigmaの変化は以下の通りである。
mu

sigma

訓練データに対するmuの値は0近傍を推移しているがepochの後半部分で振動しており、テストデータに対するそれは終始振動していることが分かる。sigmaの方は1からは程遠い値である。このときのトータルのlossは以下のようになる。
loss


実験-2


 次に、こちらで紹介されている実装で実験を行った。ソースコード内にnet_2.pyとして保存してある。epoch=100、dimz=100とした結果を以下に示す。
mu

sigma

muもsigmaもかなり改善されていることが分かる。muは0に、sigmaは1に近づいている。このときのトータルのlossは以下のようになる。

こちらも少しだけ改善されている。

生成画像の比較


 実験-1と2で、標準正規分布から生成した同じ乱数値から復号された画像を示す。
実験-1

実験-2

後者の方が数字であることを判別できるので精度は良くなっている(と思う)。

まとめ


 今回は、Chainerのサンプルコードに手を加えて実験を行った。muは0に、sigmaは1に近い方が復号化される画像の精度も上がることを示した。VAEを使用する場合、lossだけでなく、muやsigmaの値の推移も確認した方が良い。

2018年4月22日日曜日

Variational Auto Encoder 〜その2〜

はじめに


 先のページで、Variational Auto Encoder(VAE)をBayes推論の枠組みで解説し、Chainerのサンプルコードをみた。今回は、サンプルコードを実際に動かし、その動作を確認する。

前回の補足


 前回の式(16)は、観測値$X$の下での未知変数$\vec{x}$の実現確率であった。 \begin{equation} p(\vec{x}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(\vec{x}|\vec{z}) \end{equation} ここで、$z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d$である。$\vec{x}$の各成分が独立同分布で生成されると仮定すると、以下のように書き換えることができる。 \begin{equation} \prod_{m=1}^M p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\prod_{m=1}^M p(x_m|\vec{z}) \end{equation} すなわち、要素$x_m$ごとに次式が成り立つ。 \begin{equation} p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;p(x_m|\vec{z}) \end{equation} $p(x_m|\vec{z})$としてBernoulli分布を仮定したから \begin{equation} p(x_{m}|X)\approx \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\mathrm{Bern}(x_m|\eta_{\theta,m}(\vec{z})) \end{equation} となる。確率分布$p(x_{m}|X)$の下での$x_m$の期待値は \begin{eqnarray} <x_m>&=&\sum_{x_m=0,1} x_m\;p(x_{m}|X) \\ &\approx& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\sum_{x_m=0,1} x_m \mathrm{Bern}(x_m|\eta_{\theta,m}(\vec{z})) \\ &=& \int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\;\eta_{\theta,m}(\vec{z}) \end{eqnarray} となる。$\vec{z}$は$\vec{\epsilon}$に依存する項であることに注意する。ここまでの議論で言えることは、復号化した結果を得るには、$\eta_{\theta,m}(\vec{z})$を標準正規分布に従ってサンプリングすれば良いということである。$\eta_{\theta,m}(\vec{z})$はDecoderの出力である。

Chainerの実装の確認


 前回は、net.pyを見たので、今回はtrain_vae.pyを見る。trainerを用いた実装部分は特に指摘することはないので、結果を描画している箇所を解説する。最初は訓練データに関わる描画部分である。 適当に画像を9枚選択し、これを関数__call__で符号化・復号化している。前回指摘したように、復号化の際、$\sigma_{\phi,d}(X)$は無視されている。計算のあと元画像と復号化画像を保存している。epochを100とした結果は以下の通りである。

訓練画像


復号化した訓練画像

次はテスト画像に関わる部分である。 ここでも適当に9枚の画像を選択し、元画像と復号化画像を保存している。epochを100とした結果は以下の通りである。

テスト画像


復号化したテスト画像

最後に、標準正規分布に従う値$\vec{z}$から復号化する部分である。 9個の乱数$\vec{z}$を作り、関数decodeを呼び出している。epochを100とした結果は以下の通りである。
訓練・テスト画像を符号化・復号化した結果と比べるとかなり精度の悪い結果である。関数decodeは$\eta_{\theta,m}(\vec{z})$を出力する。上で見たように本来の$\vec{z}$は、$z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d$として与えらるべきものである。精度が悪いのは、$\vec{z}$を標準正規分布に置き換えたことが原因である。ところで、勾配降下法で最小にすべき式の1つが次式であった(前回の式(19))。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \end{equation} これを十分小さくできれば、すなわち、$\sigma_{\phi,d}\rightarrow 1, \mu_{\phi,d}\rightarrow 0$とできれば、標準正規分布による置き換えは意味があるものになる。残念ながらepochを1000としても大して精度は良くならない。Chainerのサンプル実装を変更する必要があるかもしれない。

まとめ


 今回は、$<x_m>$が$\eta_{\theta,m}(\vec{z})$を求めることに帰着すること示し、Chainerのサンプルコードの計算結果を考察した。 さらに、標準正規分布による$\vec{z}$から計算した画像の精度が悪い理由についても述べた。改善するには、epoch数を増やすだけではなくコードの見直し(多層化、初期化関数の変更)も必要であろう。次回はこの辺りのことについてまとめたい。

2018年4月15日日曜日

Variational Auto Encoder

はじめに


 Variational Auto Encoder(VAE)をBayes推論の枠組みで解説し、Chainerのサンプルコードを読解する。

問題設定


 観測値$X=\{\vec{x}_1,\cdots,\vec{x}_N\}$が与えられたとき、未知の値$\vec{x}_*$を生成する確率分布$p(\vec{x}_*|X)$を求めたい。潜在変数$\vec{z}$を導入し、$X$と$\vec{z}$の同時確率分布$p(X,\vec{z})$を考え、Bayesの定理を適用すると次式を得る。 \begin{equation} p(\vec{z}|X) = \frac{p(X|\vec{z})p(\vec{z})}{p(X)} \label{eq9} \end{equation} 事後確率$p(\vec{z}|X)$が求まれば、次式により$\vec{x}_*$を生成する確率分布を求めることができる。 \begin{equation} p(\vec{x}_*|X)=\int d\vec{z}\;p(\vec{x}_*|\vec{z})p(\vec{z}|X) \label{eq3} \end{equation} 事後確率$p(\vec{z}|X)$を求めることが目的である。

最適化すべき量


 $p(\vec{z}|X)$を直接求めることはせず、パラメータ$\phi$を持つ関数$q_{\phi}(\vec{z}|X)$を導入し、次のKullback Leibler divergenceを最小にすることを考える。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X) \ln{ \frac{ q_{\phi}(\vec{z}|X) } { p(\vec{z}|X) } } \end{equation} これを変形すると次式を得る。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]+\ln{p(X)} \label{eq1} \end{equation} ただし、式変形の途中で式(\ref{eq9})を用いた。式(\ref{eq1})右辺にある$\ln{p(X)}$は$\phi$に依存せず、観測値だけから決まる定数である。従って、次式が成り立つ。 \begin{equation} \min_{\phi} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}|X) \right] = \min_{\phi} {\left[ D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]-E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right] \right] } \label{eq2} \end{equation} 式(\ref{eq2})の右辺第1項を小さく、第2項の期待値を大きくすれば良い。第1項は$q_{\phi}(\vec{z}|X)$をできるだけ$p(\vec{z})$に近い形の分布にすることを要求し、この分布の下で対数尤度$\ln{p(X|\vec{z})}$の期待値を大きくすることを第2項は要求する。第1項は正則化項に相当する。

KL divergenceの計算


 式(\ref{eq2})の右辺第1項を考える。いま次の仮定をおく。 \begin{eqnarray} q_{\phi}(\vec{z}|X)&=&\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X),\Sigma_{\phi}(X)) \\ p(\vec{z})&=&\mathcal{N}(\vec{z}|\vec{0},I_D) \end{eqnarray} ここで、$\vec{z}$の次元を$D$とした。$I_D$は$D\times D$の単位行列である。どちらの分布も正規分布とし、$q_{\phi}(\vec{z}|X)$の平均と共分散行列は$\phi$と$X$から決まる量とする。これらは、入力$X$、パラメータ$\phi$のニューラルネットワークを用いて計算される。一方、$p(\vec{z})$の方は平均0、分散1の標準正規分布である。このとき、$D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]$は解析的に計算することができる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]=\frac{1}{2}\left[ -\ln{|\Sigma_{\phi}(X)|} -D +\mathrm{Tr}\left(\Sigma_{\phi}(X)\right)+\vec{\mu}_{\phi}(X)^T\vec{\mu}_{\phi}(X) \right] \label{eq4} \end{equation}

ここまでの処理の流れ


 式(\ref{eq2})の最適化を行う際に勾配降下法を用いる。処理の流れは以下のようになる(下図参照)。
分布$q_{\phi}(\vec{z}|X)$は$X$から$\vec{z}$を生成するEncoder、$p(X|\vec{z})$は$\vec{z}$から$X$を生成するDecoderとみなすことができる。青色で示した部分は最小化すべき量であり、赤字はサンプリングするステップである。勾配降下法を実現するには、誤差逆伝播ができなければならない。$q_{\phi}(\vec{z}|X)$はその$\phi$依存性のため誤差逆伝播時の微分鎖の中に組み込まれるが、サンプリグという処理の勾配を定義することができない。対数尤度の期待値の計算に工夫が必要である。

対数尤度の期待値の計算


 計算したい式は次式である。 \begin{equation} E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]=\int d\vec{z}\;q_{\phi}(\vec{z}|X)\ln{p(X|\vec{z})} \end{equation} この式に再パラメータ化トリック(re-parametrization trick)を適用する。すなわち \begin{equation} \vec{z}\sim\mathcal{N}(\vec{z}|\vec{\mu}_{\phi}(X),\Sigma_{\phi}(X)) \end{equation} の代わりに \begin{eqnarray} \vec{\epsilon}&\sim&\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\\ \vec{z}&=& \vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon} \label{eq7} \end{eqnarray} を用いてサンプリングを行う。これを用いて期待値を書き直すと次式を得る。 \begin{equation} E_{q_{\phi}(\vec{z}|X)}\left[\ln{p(X|\vec{z})}\right]=\int d\vec{\epsilon}\;\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)\ln{p(X| \vec{z}=\vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon})} \end{equation} このときの処理の流れは以下のようになる(下図参照)。
上図であれば誤差逆伝播が可能となる。

未知変数の生成


 未知変数を生成する確率分布は次式で与えられた。 \begin{equation} p(\vec{x}_*|X)=\int d\vec{z}\;p(\vec{x}_*|\vec{z})p(\vec{z}|X) \end{equation} 事後確率$p(\vec{z}|X)$の近似解$q_{\phi}(\vec{z}|X)$を用いると \begin{equation} p(\vec{x}_*|X)\approx\int d\vec{z}q_{\phi}(\vec{z}|X)p(\vec{x}_*|\vec{z}) \end{equation} を得る。先と同様に再パラメータ化トリックを適用すると \begin{equation} p(\vec{x}_*|X)\approx\int d\vec{\epsilon}\mathcal{N}(\vec{\epsilon}|\vec{0},I_D)p(\vec{x}_*| \vec{z}=\vec{\mu}_{\phi}(X)+\Sigma_{\phi}^{1/2}(X)\vec{\epsilon}) \end{equation} となる。

Chainer実装の確認


 ChainerのサンプルコードにVAEがある。これはMNISTデータセットにVAEを適用したものである。MNISTは2値画像であるから$\vec{x}_n$として0と1が784(=28$\times$28)個並んだベクトルを考えることになる。実際にコードを見て行く前に先に導出した式をもう少し詳細に計算しておく。

 最初に$\vec{\mu}_{\phi}(X)$と$\Sigma_{\phi}(X)$を次のように置く。 \begin{eqnarray} \vec{\mu}_{\phi}(X)&=&(\mu_{\phi,1}(X),\cdots,\mu_{\phi,D}(X))^T \label{eq5}\\ \Sigma_{\phi}(X)&=&\mathrm{diag}(\sigma^2_{\phi,1}(X),\cdots,\sigma^2_{\phi,D}(X)) \label{eq6} \end{eqnarray} このとき式(\ref{eq4})は次式となる。 \begin{equation} D_{KL} \left[ q_{\phi}(\vec{z}|X)||p(\vec{z}) \right]= \frac{1}{2} \sum_{d=1}^{D}\left\{ -\ln{\sigma^2_{\phi,d}(X)}-1+\sigma^2_{\phi,d}(X)+\mu_{\phi,d}^2(X) \right\} \label{eq8} \end{equation} また、$\vec{z}$の成分は次式で与えられる。 \begin{equation} z_d=\mu_{\phi,d}(X)+\sigma_{\phi,d}(X)\epsilon_d \end{equation} 観測値が独立同分布に従うと仮定すると対数尤度は以下のように変形される。 \begin{equation} \ln{p(X|\vec{z})}= \sum^{N}_{n=1}\ln{p(\vec{x}_n|\vec{z})} \end{equation} $\vec{x}_n$の次元数を$M$(=784)とすると \begin{equation} \ln{p(\vec{x}_n|\vec{z})}=\sum_{m=1}^{M}\ln{p(x_{n,m}|\vec{z})} \end{equation} となる。いま考える画像は0と1から構成される。従って、$p(x_{n,m}|\vec{z})$として0と1を生成する確率分布であるBernoulli分布を仮定する。 \begin{eqnarray} p(x_{n,m}|\vec{z})&=&\mathrm{Bern}\left(x_{n,m}|\eta_{\theta,m}\left(\vec{z}\right)\right) \\ \mathrm{Bern}(x|\eta)&=&\eta^{x}(1-\eta)^{1-x} \end{eqnarray} $\eta_{\theta,m}\left(\vec{z}\right)$は、入力$\vec{z}$、パラメータ$\theta$のネットワークで学習される量である。 以上を踏まえて処理の流れを書き直すと下図となる。

以上で準備が整ったので順にコードを見ていく。まず最初にネットワークを定義したクラスVAEを見る。コンストラクタは以下の通り。 Encoderとして$\vec{\mu}_{\phi}(X)$と$\Sigma_{\phi}(X)$を生成する層がそれぞれ1層ずつ定義されている。Decoderとして$\vec{\eta}_{\theta}(\vec{z})$を生成する2層が定義されている。次に関数encodeを見る。 $\mu_{\phi,d}$と$\ln{\sigma_{\phi,d}^2}$を生成する処理が記述されている。次は関数decodeである。 $\vec{z}$を受け取り$\vec{\eta}_{\theta}(\vec{z})$を返す処理が記述されている。次は__call__である。 Encoderで計算した平均を入力としてDecoderを呼び出している。分散を無視して復号化している。次はget_loss_funcである。
  • 13行目:Encoderで平均と対数分散を計算する。
  • 17行目:ここから始まるループは$\vec{z}$のサンプリングのためのものである。
  • 18行目:正規分布からサンプリングする。関数F.gaussianの中で再パラメータ化トリックが実行されていることに注意する。
  • 19行目:Bernoulli分布の対数にマイナスを付けたものが計算される。
  • 23行目:式(\ref{eq8})が計算される。

  • まとめ


     今回は、VAEをBayes推論の枠組みで解説し、Chainerのサンプルコードを見た。ニューラルネットワークで計算されるのは確率分布のパラメータであることを明確に示した。言い換えると他の手法でパラメータを計算できるのであればそれでも構わないということである。計算の仮定で確率分布にいくつかの仮定(ガウス分布やBernoulli分布)をしていることに注意しなければならない。今回Bernoulli分布を使用したのはターゲットとした観測値が0と1から構成される2値画像であるためである。2値でない観測値を対象とするならそれに見合った確率分布を導入する必要がある。次回はChainerのコードを動かして得られる結果を考察したい。

    参考文献


  • Tutorial on Variational Autoencoders
  • Pattern Recognition and Machine Learning
  • ベイズ推論による機械学習入門
  • Variational Autoencoder徹底解説
  • 2016年7月24日日曜日

    LSTMによる正弦波の予測 〜 Chainerによる実装 〜


    はじめに


     「RNNにsin波を学習させて予測してみた」ではTensorflowを使って、「深層学習ライブラリKerasでRNNを使ってsin波予測」ではKerasを使って、RNNによる正弦波の学習・予測が行われている。ここでは同じことをChainerを使って実装する。

    ネットワークの構造


     実装は以下の通りである。

    -- lstm.py -- コンストラクタの引数の意味は以下の通り。

    引数名 意味 デフォルト値
    in_units 入力層のユニット数 1
    hidden_units 隠れ層のユニット数 2
    out_units 出力層のユニット数 1

    隠れ層の各ユニットはLSTM(Long Short Term Memory)、損失関数は2乗平均誤差である。in_unitsout_unitsは1に固定し(実数値を1つ受け取り実数値を1つ返す)、hidden_unitsの値を変えた時の精度の変化を見る(後述)。

    訓練データの作成


     実装は以下の通りである。

    -- make_data.py -- コンストラクタの引数の意味は以下の通り。

    引数名 意味
    steps_per_cycle 正弦波1周期の分割数
    number_of_cycles 生成する周期の数

    steps_per_cycleを50に、number_of_cyclesを100としてデータを作成した。関数makeは全データを1次元の配列(np.ndarray)に入れて返す。関数make_mini_batchは、長さlength_of_sequenceのデータをバッチサイズ分(mini_batch_size)だけ切り出す。

    訓練


     実装は以下の通りである。

    -- train.py -- プログラム冒頭の定数の意味は以下の通り。

    引数名 意味
    IN_UNITS 入力層のユニット数
    HIDDEN_UNITS 隠れ層のユニット数
    OUT_UNITS 出力層のユニット数
    TRAINING_EPOCHS 訓練時のエポック数
    DISPLAY_EPOCH ログを出力するタイミング
    MINI_BATCH_SIZE ミニバッチの数
    LENGTH_OF_SEQUENCE 学習する実数値のシーケンスの長さ
    STEPS_PER_CYCLE 正弦波1周期の分割数
    NUMBER_OF_CYCLES 周期の数

    関数compute_lossでは、長さLENGTH_OF_SEQUENCEの配列をMINI_BATCH_SIZE個一括して処理を行っている(下図参照)。
    また、LENGTH_OF_SEQUENCE$\times$MINI_BATCH_SIZE個のデータの塊を1epochと定義した。以下を実行する。 DISPLAY_EPOCH(上の場合は10)epoch ごとに訓練時の損失値が表示される。上記のパラメータのとき計算時間は27分ほどであった。

    学習曲線


     下図はHIDDEN_UNITSの値を2,3,4,5に変えた時の学習曲線である。縦軸は損失(2乗平均誤差)の対数($\log_{10}$)を取った値である。10epoch毎に1回ずつ表示していので横軸の値を10倍したものが実際のepoch数である。
    隠れ層のユニット数を増やすと損失は小さくなるが、3以上になるとそれほど差はないことが分る。

    予測


     予測時のスクリプトは以下の通りである。

    -- predict.py -- 上記のスクリプトでしていることを以下に図示する。
    ここで、input_seqは最初に与える波形、pre_lengthはその波形をもとに予測する波形の長さを表す。例えば、input_seqの長さを4、pre_lengthを3とすると予測手順は以下のようになる。
    1. 4つの実数値の並びから最後の値を予測する(これを4とする)。
    2. 最初の値を捨てて4を最後尾に追加し、これを使って最後の値を予測する(これを5とする)。
    3. 最初の値を捨てて5を最後尾に追加し、これを使って最後の値を予測する(これを6とする)。
    4. (4,5,6)が求める予測値である。
    input_seqの長さを50として、pre_lengthを50とした時の結果を以下に示す。
    input_seqの長さを25として、pre_lengthを75とした時の結果を以下に示す。
    隠れ層のユニットの数が多いほど精度は良い。しかし、いずれのユニット数の場合も、最初に与える波形から遠ざかるにつれて精度は悪くなる。

    2016年7月18日月曜日

    Fully Convolutional Networks 〜 Chainerによる実装 再考2 〜


    はじめに


     先のページで、fcn32sを Chainer を使って実装した。今回は、fcn32sの訓練済みモデルを使って、fcn16sの学習を行う。

    学習曲線


     最初に学習曲線を示す。

    Accuracy
    Loss
    今回はEpoch60回で打ち切った(かろうじて収束していないように見えるが課金量の関係である)。テスト画像の正解率は94%程度となり前回の92%より向上した。このあとfcn8sを行うことになる。

    計算機環境


     前回と同じく、Amazon EC2 にある g2.2xlarge を利用した。GPU を搭載したインスタンスである。今回もミニバッチ処理ができないので1枚ずつ学習する。

    データセット


     訓練画像、テスト画像ともに前回と全く同じである。

    ネットワークの構造


     ネットワークの構造は以下の通りである。
    前回と同じである。上記の表では224$\times$224の正方形画像を想定してその1辺の長さだけをinputに記してあるが、縦横の長さは任意で構わない。表の各項目の意味は以下の通りである。
    1. name: 層の名前
    2. input: 入力featureマップの1辺のサイズ
    3. in_channels: 入力featureマップ数
    4. out_channels: 出力featureマップ数
    5. ksize: カーネルのサイズ
    6. stride: ストライドのサイズ
    7. pad: paddingのサイズ
    8. output: 出力featureマップの1辺のサイズ
    今回はfcn16sを実装するのでpool4、pool5 の出力に接続する層(score-pool4、score-pool5)が必要である。これを通ったあとの出力データ名をそれぞれp4、p5とする。p4とp5のfeatureマップの数はクラスの数に等しくなる。
    この表も、正方形画像を想定してその1辺の長さをinputに記してあるが縦横任意で構わない。出力p5はbilinear補間で2倍され、p4に加算される。そのあと、元画像と同じサイズに拡大され、ラベル画像と比較される。


    ネットワークの実装


     上記の構造をそのままChainerで記述する。

    -- myfcn_16s_with_any_size.py -- 物体の境界線には-1を配置し、softmax_cross_entropyの計算時に境界の寄与を無視するようにした。また、関数F.accuracyの引数 ignore_labelを使えば境界線上の画素を除いてaccuracyを計算することができる。112行目から113行目にかけて記述した関数でp5を2倍に拡大し、116行目でp4に加算している。122行目から123行目にかけて記述した関数で加算結果を16倍し、126行目でラベル画像との間のクロスエントロピーを計算している。

    訓練


     訓練時のスクリプトは以下の通りである。

    -- train_16s.py -- 関数 copy_model は前回と同じである。 train_16s.pyでしていることは、
    1. オブジェクトmini_batch_loader_with_any_sizeを作る。
    2. fcn32sで学習したモデルを読み込む。
    3. MyFcn16sWithAnySize.pyで定義されたネットワークのオブジェクトを作る。
    4. 学習済みモデルのパラメータをこれにコピーする。
    5. 最適化アルゴリズムとして MomentumSGD を選択する。
    6. あとは、通常の訓練過程である。
    mini_batch_loader_with_any_size.pyは前回と同じである。

    結果画像


     以下にテスト画像に適用した結果を示す。左の列はfcn32sの結果、中央の列はfcn16sの結果、右の列はGround Truthである。下に記した数値はaccuracyである。ここでaccuracyとは、1枚の画像の中で何%の画素が正解したかを計算したものである(境界線上の画素は除く)。

    一番最後の画像以外は、精度は向上している。

    ダウンロード


     ここの、タグが2016-07-17のものです。

    2016年7月16日土曜日

    Fully Convolutional Networks by Chainer 〜 revisit 1 〜

    in Japanese

    Introduction


     In the previous post, a simplified Fully Convolutional Network (FCN) was implemented by means of Chainer. Its accuracy is lower than that of the original FCN. In this post, I have remarkably improved the accuracy by the following modifications:
    1. The original implementation of the FCN is written with Caffe. Caffe supports a "CropLayer" which the original code uses. Unfortunately Chainer does not support layers corresponding to the "CropLayer." I have noticed, however, that the similar outcome is achieved by utilizing an argument "outsize" which is passed to a function "chainer.functions.deconvolution_2d." The use of the function allows my code to accept inputs of any size. In the previous post, an input was made square with 224$\times$224 in advance.
    2. In the previous post, I didn't understand the fact that an image with information on labels in VOC dataset is not RGB mode, but index one. Therefore, I spent time to implement a function by which RGB is converted to an integer value. In this post, an index image is read using index mode directly.
    3. The procedures of the construction of the original FCN is fcn32s $\rightarrow$ fcn16s $\rightarrow$ fcn8s. The accuracy of the FCN depends mostly on that of the fcn32s, and the contributions from the fcn16s and the fcn8s are at most 1%. Though I showed the implementation of the fcn8s in the previous post, in this post the code of the fcn32s is demonstrated.
    4. As a mini batch learning was employed in the previous post, "chainer.links.BatchNormalization" was able to be introduced. In this post, there is no room to introduce the normalization because of the use of a one-by-one learning.
    5. In the previous post, weights in deconvolutional layers were learned using "chainer.links.Deconvolution2D." In this post, the deconvolutional layers are fixed to the simple bilinear interpolation.

    Learning Curve


     The resultant learning curves are as follows:

    Accuracy
    Loss
    Both curves have the ideal behaviour. The accuracy for test images achieves about 92%, which is much better than the previous accuracy. (If the trained model of fcn32s is succeeded by fcn16s and fcn8s, the accuracy will be about 93%.) More detail is provided below.

    Computation Environment


     The same instance as the previous post, g2.2xlarge in the Amazon EC2, is used. As input images have various sizes, there is no choice but to run the one-by-one learning. The mini batch learning is not able to be used. It takes 3.5 days for the above learning curves to converge.

    Dataset


     Just as before, the FCN is trained on a dataset VOC2012 which includes ground truth segmentations. The number of images is 2913. I divided them by the split ratio 4:1. The former set of images corresponds to a training dataset, and the latter a testing one.

    number of train number of test
    2330 583

    As described above, all images are not resized in advance. The number of the categories is 20+1(background).

    label category
    0 background
    1 aeroplane
    2 bicycle
    3 bird
    4 boat
    5 bottle
    6 bus
    7 car
    8 cat
    9 chair
    10 cow
    11 diningtable
    12 dog
    13 horse
    14 motorbike
    15 person
    16 potted plant
    17 sheep
    18 sofa
    19 train
    20 tv/monitor


    Network Structure


     The detailed structure of a network is as follows:
    In the original FCN, there are some layers that follow the layer pool5. For simplicity, these layers are again removed. Under the assumption that an image is a square with 224 $\times$ 224, one side length of an image is written to a column "input" of the table shown above. It should be noticed that the FCN implemented in this post is actually permitted to accept inputs of any size. Each column of the table indicates:
    1. name: a layer name
    2. input: a size of an input feature map
    3. in_channels: the number of input feature maps
    4. out_channels: the number of output feature maps
    5. ksize: a kernel size
    6. stride: a stride size
    7. pad: a padding size
    8. output: a size of an output feature map
    Because the fcn32s is implemented, only score-pool5 following pool5 is required. I refer to an output of score-pool5 as p5. The number of feature maps of p5 is equal to that of categories (21).
    Though a column "input" of that figure also shows one side length of an image under the same assumption as shown above, inputs of any size are permitted. After upsampling p5 back to input size, it is compared with a label image.


    Implementation of Network


     The network can be written in Chainer like this:

    -- myfcn_32s_with_any_size.py -- By assigning a label of -1 to pixels on the borderline between objects, we can ignore contribution from those pixels when calculating softmax_cross_entropy (see the Chainer's specification for details). Using an argument "ignore_label" passed to a function "F.accuracy" allows us to calculate the accuracy without contribution from pixels on the borderlines. The function "F.deconvolution_2d" shown from the 99th line to the 100th line upsamples p5 back to the input size.

    Training


     The script for training is as follows:

    -- train_32s_with_any_size.py -- I used copy_model and VGGNet.py described in the previous post. The procedures shown in the script train_32s_with_any_size.py are as follows:
    1. make an instance of the type MiniBatchLoaderWithAnySize
    2. make an instance of the type VGGNet
    3. make an instance of the type MyFcn32sWithAnySize
    4. copy parameters of the instance of VGGNet to those of that of MyFcn32sWithAnySize
    5. select MomentumSGD as an optimization algorithm
    6. run a loop to train the net
    The script mini_batch_loader_with_any_size.py is written like this:

    -- mini_batch_loader_with_any_size.py --
    The function load_voc_label shown from the 91th line to the 96th line loads a label image using the index mode and replaces 255 with -1. To load all the training data on the GPU memory at a time causes the error "cudaErrorMemoryAllocation: out of memory" to occur. Therefore, only one pair of an image and a label is loaded every time a training procedure requires it.

    Prediction Results


     The left column shows prediction images overlaid on testing images. A value under each image indicates an accuracy which is defined as the percentage of the pixels which have labels classified correctly that are not on the borderlines between objects. The right column shows the ground truths.

    I think that if fcn16s and fcn8s follow fcn32s, outlines of the objects will be much finer.

    Download


     You can download tag 2016-07-09

    2016年7月11日月曜日

    Fully Convolutional Networks 〜 Chainerによる実装 再考1 〜

    in English

    はじめに


     先のページで、簡易化したFully Convolutional Networks(FCN)を Chainer を使って実装した。残念ながら、その精度は文献のものより低かった。今回以下の改良を行ったところ、格段に精度は向上した。
    1. FCNのソースはCaffeである。CaffeにはCropLayerなる層が実装されており、FCNはこれを利用している。Chainerにはこれに相当するものがないが、chainer.functions.deconvolution_2dのoutsizeという引数を利用すれば同じような効果が得られることが分った。この関数の利用により、任意サイズの画像を受け付けることができるようになった。従って、入力画像に対しては何も手を加えていない。前回は$224\times 224$にリサイズしたのであった。
    2. ラベルの情報を持つ画像(ラベル画像)はRGB画像でなく、インデックス画像である。前回はこれを理解しておらず、わざわざRGBから整数値に変換する関数を介在させていた。今回はインデックス画像としてラベル画像を読み込むようにした。
    3. FCNの正式な構築手順は、fcn32s $\rightarrow$ fcn16s $\rightarrow$ fcn8sである。FCNの精度は、ほぼfcn32sの精度で決まっており、fcn16s/fcn8sからの寄与は1%程度である。前回はfcn8sを実装したが、今回はfcn32sの実装を示す。
    4. 前回はミニバッチ処理を行ったので関数chainer.links.BatchNormalizationを導入したが、今回は1枚ずつ処理するのでこの層を削除した。
    5. 前回はchainer.links.Deconvolution2Dを利用して拡大時のパラメータも学習していた。今回はchainer.functions.deconvolution_2dを利用し、単純なbilinear補間による拡大を採用した。

    学習曲線


     最初に学習曲線を示す。

    Accuracy
    Loss
    どちらも理想的な曲線となった。またテスト画像の正解率は92%弱となり前回よりも格段に向上した。このあとfcn16s/fcn8sを行えば93%程度にはなりそうである。以下詳細を記す。

    計算機環境


     前回と同じく、Amazon EC2 にある g2.2xlarge を利用した。GPU を搭載したインスタンスである。 ただし、今回はミニバッチ処理ができないので1枚ずつ学習する。上の結果を得るのに3.5日ほどかかっている。課金量も相当である。

    データセット


     前回と同じく、データセットは VOC2012 である。領域分割の教師データの数は2913枚、これを4:1に分割し、前者を訓練データ、後者をテストデータとした。

    number of train number of test
    2330 583

    前回は、訓練データ数は10で、テストデータ数は5で割り切れるように端数を切り捨てたが、今回はミニバッチ処理を行わないので、そのような操作はしていない。上述したように画像のサイズもそのままである。 カテゴリ数は前回と同じく20+1(背景)である。

    label category
    0 background
    1 aeroplane
    2 bicycle
    3 bird
    4 boat
    5 bottle
    6 bus
    7 car
    8 cat
    9 chair
    10 cow
    11 diningtable
    12 dog
    13 horse
    14 motorbike
    15 person
    16 potted plant
    17 sheep
    18 sofa
    19 train
    20 tv/monitor

    ネットワークの構造


     ネットワークの構造は以下の通りである。
    これは前回と同じである。文献では、pool5 の後ろにfc6とfc7があるが、これらを入れると拡大時の処理が煩雑となる(というかよく分らない)ので、これら2層は残念ながら今回も削除した。上記の表では224$\times$224の正方形画像を想定してその1辺の長さだけをinputに記してあるが、縦横の長さは任意で構わない。表の各項目の意味は以下の通りである。
    1. name: 層の名前
    2. input: 入力featureマップの1辺のサイズ
    3. in_channels: 入力featureマップ数
    4. out_channels: 出力featureマップ数
    5. ksize: カーネルのサイズ
    6. stride: ストライドのサイズ
    7. pad: paddingのサイズ
    8. output: 出力featureマップの1辺のサイズ
    今回はfcn32sを実装するのでpool5 の出力に接続する層(score-pool5)だけが必要である。これを通ったあとの出力データ名をp5とする。このデータのfeatureマップの数はクラスの数に等しくなる。
    この表も、正方形画像を想定してその1辺の長さをinputに記してあるが縦横任意で構わない。出力p5はbilinear補間で元画像と同じサイズに拡大され、ラベル画像と比較される。


    ネットワークの実装


     上記の構造をそのままChainerで記述する。

    -- myfcn_32s_with_any_size.py -- 物体の境界線には-1を配置し、softmax_cross_entropyの計算時に境界の寄与を無視するようにした(Chainerの仕様ではラベルが-1の画素は評価されない)。また、関数F.accuracyの引数 ignore_labelを使えば境界線上の画素を除いてaccuracyを計算することができる(前回はわざわざ実装していた)。99行目から100行目にかけて記載した関数F.deconvolution_2dでp5を32倍している。さらに、引数outsizeにtのサイズを渡すことで結果がtと同じになるようにリサイズを行う。引数padには辻褄合わせのpadding量を渡している。

    訓練


     訓練時のスクリプトは以下の通りである。

    -- train_32s_with_any_size.py -- 関数 copy_model とVGGNet.pyは前回と同じである。 train_32s_with_any_size.pyでしていることは、
    1. オブジェクトmini_batch_loader_with_any_sizeを作る。
    2. VGGNet.pyで定義されたネットワークのオブジェクトを作る。
    3. MyFcn32sWithAnySize.pyで定義されたネットワークのオブジェクトを作る。
    4. VGGNetのパラメータをこれにコピーする。
    5. 最適化アルゴリズムとして MomentumSGD を選択する。
    6. あとは、通常の訓練過程である。
    mini_batch_loader_with_any_size.pyの内容は以下の通りである。

    -- mini_batch_loader_with_any_size.py --
    91行目から96行目までの関数load_voc_labelでは、pngのインデックス画像をインデックスのまま読み込み、値が255の画素を-1に置き換えている。データを一括読み込みすると、ビデオメモリーが足りなくなるので、要求されるごとに読み込むようにしてある。

    結果画像


     以下の左の列に、テスト画像に対する予測結果をテスト画像にオーバレイした画像を示す。下に記した数値はaccuracyである。ここでaccuracyとは、1枚の画像の中で何%の画素が正解したかを計算したものである(境界線上の画素は除く)。右側の列は、Ground Truthである。

    このあと、fcn16s/fcn8sを行えば、物体の輪郭はもう少し正確になるだろう。

    ダウンロード


     ここの、タグが2016-07-09のものです。今となっては不要なファイルがたくさん残っています。