numpyだけでCNN実装

# はじめに

pythonでCNNを実装しました．
深層学習ライブラリは使用せず，numpyだけで実装しました．
教科書として[『深層学習』](http://www.amazon.co.jp/深層学習-機械学習プロフェッショナルシリーズ-岡谷-貴之/dp/4061529021)を使いました．

本記事の構成

* はじめに
* CNN
    * 畳込み層
    * プーリング層
* 学習
    * 重みの更新
    * 誤差逆伝播
* pythonでの実装
    * 畳込み層の実装
    * プーリング層の実装
* MNISTデータセットでの実験
    * 学習
    * 結果
* おわりに


# CNN

CNNとは，畳込み演算を用いた順伝播型ネットワークであり，主に画像認識に応用されています．
一般的なニューラルネットワークは，隣接層のユニットが全結合されたものですが，
CNNは，隣接層間の特定のユニットのみが結合した特別な層を持ちます．
これらの特殊な層では，**畳込み** および **プーリング** という演算を行います．
以下では，畳込みとプーリングについて説明します．

## 畳込み層

畳込みは，画像にフィルタの対応する画素同士の積をとり，その総和をとる演算です．
フィルタの濃淡パターンと類似した濃淡パターンを検出する働きがあります．
画像サイズを $W \times W$，インデックスを $(i, j)$，画素値を $x_{ij}$ で表します．
フィルタサイズを $H \times H$，インデックスを $(p, q)$，フィルタの画素値を $h_{pq}$ で表します．
畳込みは，以下の式で表されます．

```math
u_{ij} = \sum^{H-1}_{p=0} \sum^{H-1}_{q=0} x_{i+p, j+q} \, h_{pq}
```

画像内に収まる範囲内でフィルタを動かしたとき，畳込みの結果の画像サイズは以下のようになります．
ただし $\lfloor \cdot \rfloor$ は小数点以下を切り捨て，整数化する演算子です．

```math
(W - 2 \lfloor H / 2 \rfloor) \times (W - 2 \lfloor H / 2 \rfloor)
```

畳込み層では，下図のように畳込み演算を行います．
入力画像のサイズを $W \times W \times K$，畳込みフィルタのサイズを $H \times H \times K \times M$ とします．
$K$ は画像のチャネル数を，$M$ はフィルタの種類の数を表します．
$l - 1$ 層の $k$ チャネルの画像を $m$ 番目のフィルタで畳込んだ結果は以下になります．
ただし $b_{ijm}$ はバイアス，$f$ は活性化関数を表します．

```math
\begin{align}
u_{ijm} &= \sum^{K-1}_{k=0} \sum^{H-1}_{p=0} \sum^{H-1}_{q=0} z^{(l-1)}_{i+p, j+q, k} \, h_{pqkm} + b_{ijm} \\
\\
z_{ijm} &= f(u_{ijm})
\end{align}
```

$z_{ijm}$ を順番に並べたものを $M$ チャネルの画像とみなし，次の層への入力 $z^{(l)}_{ijm}$ とします．
また，フィルタをずらしながら適用するため，同一の重みが繰り返し使われます．
これを **重み共有** と呼びます．

<img width="1000" alt="conv_layer.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/beefb47a-1ed2-b110-9dc5-57e40fc71fd4.png">

上の図では，入力画像のチャネル数 $K = 3$，フィルタの種類 $M = 2$ とし，$2$ チャネルの画像を出力しています．

## プーリング層

プーリングは，画像の局所領域を1つの値にまとめる演算です．
畳込み層で抽出された特徴の位置感度を低下させ，多少の位置ずれに対してプーリング層の出力が不変になるようにします．
サイズ $W \times W \times K$ の画像上で画素 $(i, j)$ を中心とする $H \times H$ 正方領域をとり，領域内の画素の集合を $P_{ij}$ で表します．
プーリングによって得られる値 $u_{ijk}$ は以下のように表せます．

```math
u_{ijk} = \biggl(\frac{1}{H^2} \sum_{(p, q) \in P_{ij}} z^{P}_{pqk} \biggr)^{\frac{1}{P}}
```

$P = 1$ のとき，領域内の画素の平均をとるため，**平均プーリング** と呼ばれます．
$P = \infty$ のとき，領域内の画素の最大値をとるため，**最大プーリング** と呼ばれます．
下図に最大プーリングの図を示します．
領域サイズ $2 \times 2$，ストライド $2$ で $4 \times 4$ の入力画像をプーリングをしています．

<img width="400" alt="pool.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/98b4812a-decc-d2c8-c29f-db18d7df97c3.png">


# 学習

CNNの学習について説明します．
以前書いた[『pythonでニューラルネットワーク実装』](http://qiita.com/ta-ka/items/bcdfd2d9903146c51dcb)の誤差逆伝播の章が参考になると思います．

## 重みの更新

学習データから計算される出力値を教師ラベルに近づけるために，誤差関数 $E$ の最小化を考えます．
誤差関数 $E$ を重み $w$ で偏微分し，$0$ に近づくように重みを更新します．

```math
w_{new} = w_{old} - \varepsilon \frac{\partial E}{\partial w_{old}}
```

## 誤差逆伝播

誤差逆伝播の計算方法は，一般的なニューラルネットワークと変わりません．
$l+1$ 層の誤差 $\delta^{(l+1)}$ と重み $w^{(l+1)}$ および $l$ 層目からの入力の微分値の積で，$l$ 層での誤差を求めることができます．
ただし，以下の2点がCNNと全結合のニューラルネットワークと異なります．

* 畳込み層およびプーリング層で特定のユニットのみが結合している
* 畳込み層で重みが共有されている


# pythonでの実装

今回[ Chainerによる畳み込みニューラルネットワークの実装 ](http://aidiary.hatenablog.com/entry/20151007/1444223445)で利用されているネットワークを実装しました．
実装したコードは[ ここ ](https://github.com/shota-takayama/cnn)にあげてあります．
畳込み層，プーリング層に分けて実装上のポイントを紹介します．

## 畳込み層の実装

畳込み層では，まず画像の局所領域を切り出し，順番に並べたものを入力として順伝播しています．
例えば，$20 \times 12 \times 12$ の入力画像を $50 \times 20 \times 5 \times 5$ のフィルタで畳込む場合，下図のように入力画像を成形します．
成形された入力のサイズは，$64 \times 20 \times 5 \times 5$ となります．
$64$ という数値は，以下のようにして計算されます．
$(12 - \lfloor 5 / 2 \rfloor \times 2) \times (12 - \lfloor 5 / 2 \rfloor \times 2) = 64$

<img width="700" alt="im2patch.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/3e902509-25bf-2def-bae9-5b2c4acf708f.png">

成形した入力画像とフィルタの畳込みを計算します．
繰り返しになりますが，入力およびフィルタのサイズは，それぞれ $64 \times 20 \times 5 \times 5$，$50 \times 20 \times 5 \times 5$ です．
```np.tensordot(X, weight, ((1, 2, 3), (1, 2, 3)))``` により $64 \times 50$ の出力が得られます．

<img width="1000" alt="tensordot.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/1f8aaf87-bc78-a0f2-e798-888dec58990b.png">

以下のコードが，畳込み層における順伝播のキーとなる部分です．
ミニバッチで学習できるよう入力の次元数が一つ増えているので，```tensordot``` の引数は一つずれて ```axes = ((2, 3, 4), (1, 2, 3))``` となっています．

```py
    def __forward(self, X):
        s_batch, k, xh, xw = X.shape
        m = self.weight.shape[0]
        oh, ow = xh - self.kh / 2 * 2, xw - self.kw / 2 * 2
        self.__patch = self.__im2patch(X, s_batch, k, oh, ow)
        return np.tensordot(self.__patch, self.weight, ((2, 3, 4), (1, 2, 3))).swapaxes(1, 2).reshape(s_batch, m, oh, ow)


    def __im2patch(self, X, s_batch, k, oh, ow):
        patch = np.zeros((s_batch, oh * ow, k, self.kh, self.kw))
        for j in range(oh):
            for i in range(ow):
                patch[:, j * ow + i, :, :, :] = X[:, :, j:j+self.kh, i:i+self.kw]
        return patch
```

畳込み層における逆伝播では，一つ先の層の $\delta$ およびフィルタの係数の積が逆伝播させる誤差となります．
以下のように実装しました．
局所領域ごとに得られた誤差を入力画像の形に再成形しています．
これは，順伝播のとき入力画像を局所領域に切り出した処理と逆の処理になります．

```py
    def backward(self, delta, shape):
        s_batch, k, h, w = delta.shape
        delta_patch = np.tensordot(delta.reshape(s_batch, k, h * w), self.weight, (1, 0))
        return self.__patch2im(delta_patch, h, w, shape)

    def __patch2im(self, patch, h, w, shape):
        im = np.zeros(shape)
        for j in range(h):
            for i in range(w):
                im[:, :, j:j+self.kh, i:i+self.kw] += patch[:, j * w + i]
        return im

```

## プーリング層の実装

プーリング層での順伝播も入力画像を成形し，局所領域を順番に並べます．
畳込み層と異なるのは，最大値を与えた画素のインデックスを保存している点です．
これは，どの画素から値が伝播されてきたかという情報を使って，誤差を逆伝播するためです．
以下のコードが，プーリング層における順伝播のキーとなる部分です．

```py
    def forward(self, X):
        s_batch, k, h, w = X.shape
        oh, ow = (h - self.kh) / self.s + 1, (w - self.kw) / self.s + 1
        val, self.__ind = self.__max(X, s_batch, k, oh, ow)
        return val

    def __max(self, X, s_batch, k, oh, ow):
        patch = self.__im2patch(X, s_batch, k, oh, ow)
        return map(lambda _f: _f(patch, axis = 3).reshape(s_batch, k, oh, ow), [np.max, np.argmax])


    def __im2patch(self, X, s_batch, k, oh, ow):
        patch = np.zeros((s_batch, oh * ow, k, self.kh, self.kw))
        for j in range(oh):
            for i in range(ow):
                _j, _i = j * self.s, i * self.s
                patch[:, j * ow + i, :, :, :] = X[:, :, _j:_j+self.kh, _i:_i+self.kw]
        return patch.swapaxes(1, 2).reshape(s_batch, k, oh * ow, -1)

```

上述のように，プーリング層での逆伝播は，最大値を与えた画素に対してそのまま誤差を伝播しています．
以下のように実装しました．

```py
    def backward(self, X, delta, act):
        s_batch, k, h, w = X.shape
        oh, ow = delta.shape[2:]
        rh, rw = h / oh, w / ow
        ind = np.arange(s_batch * k * oh * ow) * rh * rw + self.__ind.flatten()
        return self.__backward(delta, ind, s_batch, k, h, w, oh, ow) * act.derivate(X)

    def __backward(self, delta, ind, s_batch, k, h, w, oh, ow):
        _delta = np.zeros(s_batch * k * h * w)
        _delta[ind] = delta.flatten()
        return _delta.reshape(s_batch, k, oh, ow, self.kh, self.kw).swapaxes(3, 4).reshape(s_batch, k, h, w)
```


# MNISTデータセットでの実験

MNISTデータセットを使って，手書き数字を学習しました．
層の構造は，[ Chainerによる畳み込みニューラルネットワークの実装 ](http://aidiary.hatenablog.com/entry/20151007/1444223445)と同じです．

## 学習 
各種パラメータは以下になります．
入力データ：$28 \times 28$ のグレースケール画像 $10000$ 枚
学習率：$\varepsilon = 0.005$
正則化項の係数：$\lambda = 0.0001$
学習率の減衰係数：$\gamma = 0.9$
バッチサイズ：$5$
エポック：$50$
テストデータ：入力データと同じサイズのグレースケール画像 $100$ 枚

## 結果

各エポックにおけるロスを描画したグラフが以下になります．
最終的にロスは，$0.104299490259$ まで下がりました．
また，テスト画像 $100$ 枚での識別精度は，$\boldsymbol{0.96}$ でした．

<img width="600" alt="loss.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/392123ef-78bb-b50e-173c-2cc239ed0832.png">


# おわりに

CNNを実装することができました．
今回パディングやbatch normalizationを実装できていませんが，疲れたので終わりにします．
ライブラリ作っている人たち凄いと感じました．