ニューラルネットワークを理解するためにPHPで実装してみる

こんにちは、みなさん

機械学習が盛んになってきている今ですが、なかなか自分の中でしっくりこないというか、なじまない感覚がしていました。
そもそも機械学習ってどうやって動いてるんだ？とか、ニューラルネットってなんだ？とか。
理論の解説をしている方も多くおり、機械学習のフレームワークなども多く出回っているというのに、なぜなのか？と思っていたのですが。。。

結局のところ、根っこの実装を自分でやっていないから、というのが私の思い至った結論です。
また、PHPerなんで、他の言語で実装されていても「PHPでok」状態になるので、PHPで実装して、それを通してニューラルネットの動作を理解しようと考えました。

(FANNというのがあるよと聞いたのですが、あれ、C言語。。。)

# 結果

先に結果を示しておきます。
https://github.com/niisan-tokyo/phpnn

packagistにも登録してあるので、composerで取得することもできます。

```
composer require niisan-tokyo/phpnn:dev-master
```

`src/layer/Base.php`でほとんどの処理をしています。
結果としては、ニューラルネットと行っても、単純なループを使って実装できました。
以下で書いてあるのは、自分がやったことを忘れないようにするために作ったメモです。

# ニューラルネット

## 関数展開

先に機械学習でよく使われているニューラルネットについて基礎的な部分を先に知っておく必要がありました。
ニューラルネットの概念ですが、単純に言えば関数です。
ある入力$x$に対する応答$y$とすると、

```math
y = f(x)
```
こんな感じに書けるわけです。
この$f(x)$がはじめからわかっているのであれば、何も悩む必要はありません。
その関数$f(x)$を実装してしまえば、あらゆる入力$x$に対して、正しい応答$y$を返せるからです。

まあ、そんな関数がすぐ見つかればいいのですが、多くの場合見つかりません。
そこで関数展開という考え方ができました。
関数展開というのは、よくわからない関数を、別の既知の関数で表してみようというもので、有名なものではテイラー展開とかフーリエ級数展開とかがあります。
例えばテイラー展開は以下のように表されます。

```math
f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \frac{f^{(3)}(a)}{3!}(x-a)^3 + ... = \sum_{k=0}^{\infty}\frac{f^{(k)}(a)}{k!}(x-a)^k
```
ここで$f^{(n)}(x)$は関数$f(x)$のn回微分です。

ニューラルネットというのは、実際にはこれらの関数展開の仲間と捉えることができます。

## 概念

以下のリンクも参照してください
http://d.hatena.ne.jp/ura_ra/20111026/1319642014
http://www.sist.ac.jp/~kanakubo/research/neuro.html

ニューラルネット、とくに階層型のニューラルネットを考える際は、よく層状のネットワーク図が使われていますが、それはどこにでも書いてありますし、私には絵のセンスが無いので、あえて書きません。

今、ある入力$\bf{x}$に対し出力$\bf{y}$を得る時、関数を

```math
\bf{y} = \rm{F}(\bf{x})
```
ここで、$\bf{y}$と$\bf{x}$は幾つもの情報を持っている可能性があるため、多次元のベクトルであると考えてよいです。つまり、$\bf{x} = \rm{(x_1, x_2, ..., x_i)}$と表すことができます。

さて、この関数ですが、当然内容がわからないわけです。そこで、別の関数で表し直すわけですが、ニューラルネットは次のように関数を近似しています。

```math
F(x) \sim f_n(f_{n-1}(f_{n-2}(... (f_1(\bf{x})) ... )))
```
ここで各関数$f_k$もまた、複数の情報を持つ可能性があるため、多次元のベクトルだと考えられます。
更に難しくなっていないか？とも思いますが、実は数値計算的にはそんなに難しくなかったりします。

## 順伝搬

まず$f_1(x)$を考えて見ます。$f_1$を次のようにおいてみましょう。

```math
f_1(x) = (\phi _1(u_1), \phi _2(u_2),..., \phi _m(u_m))
```
```math
u_i = \sum_jw_{ij}x_j
```
先に述べたとおり、$f_1$は複数の要素からなるベクトルでしたので、ひとつ目の式ができました。
ここで$\phi_i(u_i)$はただのスカラ関数で$u_i$もただの数値です。
そして、$u_i$はふたつ目の式にあるように、入力値$x$の各要素を適当な重み$w$を使って足しあわせて作られた新しい数です。
こうして$f_1$が求められたので、同様に$f_2, f_3...$と求めていきます。
こうすると、$f_k$は次のように表すことができます。

```math
f_k(x) = (\phi ^k_1(u_1), \phi ^k_2(u_2),..., \phi ^k_m(u_m))
```

```math
u_i = \sum_j w^{k}_{ij}x^{k-1}_j
```

```math
x^k = f_k(x^{k-1})
```

各層にある関数$\phi ^k_i (u_i)$はこちらで好きな関数を使ってやればよいです。その代わり重み$w$をちゃんと選ぶことで、近似が可能になるということになります。
ちなみに、各層の重み$w^k$を適切に選んでやることで、あらゆる関数を近似できるということらしいです。
http://qiita.com/Ugo-Nama/items/04814a13c9ea84978a4c

## 実装

このようなニューラルネットを考えると、実装は割とシンプルに行けそうです。

```php
public function something($x)
{
    $ret = $x;
    foreach ($this->layers as $func) {
        $ret = $func->prop($ret);
    }

    return $ret;
}
```
かなり大雑把ですが、こんな感じに作れば良いと考えました。
あるレイヤーから出力された結果を、次のレイヤーに入力し、そこから出てきた出力をまた。。。という形で次々に出力されていきます。
これが順伝播です。

ロジックを見る限り、各層ごとにオブジェクトを作ってやればよいでしょう。

```php
abstract class Base
{
    private $input_dim = 0;
    private $output_dim = 0;
    private $matrix = [];

    public function prop($states)
    {
        $ret = [];
        for ($i = 0; $i < $this->output_dim; $i++) {
            $temp = 0;
            for ($j = 0; $j < $this->input_dim; $j++) {
                $temp += $this->matrix[$i][$j] * $states[$j];
            }
            $ret[$i] = $this->activate($temp);
        }

        return $ret;
    }

    protected abstract function activate($val);
}
```
こんなもんでいいでしょう。
ここで関数$\phi _i$に当たるものを`activate`というプログラム上の関数に置き換えています。
別に各要素$i$で違う関数を使わなければならないという制約もないので、同じ関数を使用することにしています。
また、各層が持つ重みに関しては`$this->matrix`という２次元配列のプロパティで表現しています。

例えば、中間層の関数によく使われるreluという関数がありますが、これは入力値が0以下であれば0, そうでなければ入力値をそのまま返すという関数です。
実装はこんな感じになります。

```php
class Relu extends Base
{
    protected function activate($val)
    {
        return ($val > 0) ? $val : 0;
    }
}
```
数式で表すとなんだかごちゃごちゃしていて難しそうですが、分解して実装に落としこむと、随分と簡単に見えます。
流石、PHPです。

# 学習

## 重み

こうしてニューラルネットをPHPで実装できたわけですが、当然このままでいいわけではありません。
層ごとにオブジェクトを作ることも、各層が持つ関数も定めることができましたが、一つ決まっていないものがあります。
各層に対して、入力値を各関数の変数に変換するための重み$w$です。
これが **適切であれば** ニューラルネットが関数の近似として良い物になるということですが、逆にそうでなければあまり役には立たないということです。

とはいえ、重み$w$が無いと、そもそも順伝播の計算ができません。
そこで、まず適当に重みを設定して、後で正しい値と照らし合わせ、重みを適宜修正することで、正しい関数形に近づけることにします。

## 重みの修正 = 学習

入力と出力がわかっている場合、こちらで適当に作ったニューラルネットの出力との間に、どれだけの誤差があるかを計算してやることで、その誤差を修正するように重みを修正してやることで、関数形を正しい形に近づけることができます。

こちらで用意した仮の関数$F(x)$と実際の値$z$との間の誤差を$L(F, z)$とします。
この誤差関数$L$が小さくなるように$w$を修正するのですが、そこで数学の力を借ります。

```math
\frac{\partial L}{\partial w^n_{ij}} =  \frac{\partial L}{\partial F}
 \frac{\partial F}{\partial w^n_{ij}} =
 \frac{\partial L}{\partial F} \frac{\partial f_{n}(x^n)}{\partial w^n_{ij}}
```
ここで、$x^n$とは$x$のn乗というわけではなく、n層目への入力という意味です。そして、$\frac{\partial L}{\partial w^n_{ij}}$と言うのは$L$を$w^n_{ij}$で偏微分したもので、$w^n_{ij}$をちょっとだけ変化させた時、$L$がどれだけ動くかを表しています。
ところで、$f_n(x^n) = (\phi ^n_1 (u_1), \phi ^n_2 (u_2), ... , \phi ^n_m (u_m))$の各要素を見ると、

```math
\phi^n_{i}(u_i) = \phi ^n_{i}\left( \sum_k w_{ik}x^{n-1}_k \right)
```
となっているので、連鎖率( 合成関数の微分 )を使うと

```math
\frac{\partial \phi _i(u_i)}{\partial w^n_{ij}} = \phi _i'(u_i) \frac{\partial }{\partial w^n_{ij}} \left( \sum_k w^n_{ik}x^{n-1}_k \right)
= \phi _i'(u_i) x^{n-1}_j
```
と書けます。

ここまで来ると、数値計算可能になります。
誤差関数Lはこっちが勝手に決めることができるので、その微分も作れますし、$\phi_i$もやはり微分可能な関数を選んでやれば、微分も作れます。$ u_i, x^n $も順伝播の時に取得済みです。

では第$l$層も以下のように書けます

```math
\frac{\partial L}{\partial w^l_{ij}} = \frac{\partial L}{\partial F}
 \frac{\partial f_{n}}{\partial f_{n-1}}
 \frac{\partial f_{n-1}}{\partial f_{n-2}} ... 
\frac{\partial f_{l+1}}{\partial f_{l}}
\frac{\partial f_{l}}{\partial w^l_{ij}}
```
という感じになります。ところで、

```math
f_l(u^l) =f_l (w^l f_{l-1}(x^{l-1}))
```
だったので、

```math
\frac{\partial f_{l}}{\partial f_{l-1}}
= f'_l(u^l)w^l
```
ふむ。。。まだこれでは計算方法がよくわかりませんね。ところが、以下のように置くと計算方針が見えてきます。

```math
\delta _{n-l} = \delta _{n-(l+1)}\frac{\partial f_{l+1}}{\partial f_{l}} = \delta _{n-(l+1)}f'_{l+1}(u^{l+1})w^{l+1}, \ \ \  \delta _0 = \frac{\partial L}{\partial F}
```
とすると、

```math
\frac{\partial L}{\partial w^l_{ij}} = \delta _{n-l} \ \ \phi _i'(u_i) x^{l-1}_j
```
となります。
さて, $\delta_{n-l}$は$\delta_{n-(l+1)}$つまり、 一つ後の層から現在の層の $\delta$が決まることが判明しました。
ところで、先に述べたとおり、$\frac{\partial L}{\partial w^l_{ij}}$は$w_{ij}$を少しだけ動かした時に誤差がどれだけ変化するかを表しています。
この値がマイナスであれば、$w^l_{ij}$をちょっとだけ増やすことで、より誤差を減らすことができますし、その逆にプラスであれば, $w^l_{ij}$をちょっとだけ減らせば、誤差を減らすことができると予想できます。

というわけで、ようやくですが、これで重みの修正方法がわかりました

```math
w^l_{ij} \gets w^l_{ij} - \alpha \frac{\partial L}{\partial w^l_{ij}}
```
ここで$\alpha $は学習率というもので、適当な値です(適当ですが、ちゃんと選ばないとうまく収束しません。。。)。

## 実装

長くなりましたが、いよいよ実装に移りましょう。
出力層、つまり第n層で発生した誤差からn-1層、n-2層、... で使用する$\delta$を作り、その都度重みを修正していけばいいわけです。
ということで、実装は次のような感じになります。

```php
public function somethingBack($x)
{
    $ret = $x;
    foreach (array_reverse($this->layers) as $func) {
        $ret = $func->backProp($ret);
    }

    return $ret;
}
```
結局ループで使用する関数の順番を変えて、使用する関数を変えただけで、順伝播と同じ形になりました。
で、当の関数`backProp`の中身ですが、

```php
    public function backProp($states)
    {
        $ret = array_fill(0, $this->input_dim, 0);
        $back_state = $this->state_history;
        $delta = [];
        for ($i = 0; $i < $this->output_dim; $i++) {
            $delta[$i] = $this->defferential($back_input[$i]) * $states[$i];

            for ($j = 0; $j < $this->input_dim; $j++) {
                $ret[$j] += $this->matrix[$i][$j] * $delta[$i];
                $this->matrix[$i][$j] -= $this->effect * $delta[$i] * $back_state[$j];
            }
        }

        return $ret;
    }
```
こんな感じで、自身の持っている重みを修正しつつ、前の層に送るための$\delta $を作っています。
微分`defferential`は関数の微分値です。例えば次のうような関数になっています。

```php
    public function defferential($val)
    {
        if ($val > 0) {
            return 1;
        }

        return 0;
    }
```
これはReluの微分値になります。

# まとめ

ニューラルネットの資料とかは結構あるのですが、なんというか、PHPで実装したものが見当たらなかったので、作ってみた次第です。
PHPerの母国語はPHPなのです。
PythonとかC++で書かれていても、馴染みがないのです。
今まで漠然とPythonのフレームワークとかを写経していましたが、やはりこの方法が一番しっくり来ます。

PHPになって自由度が上がったので、これでいろんな問題を解いてみようと思います。

# 参考資料

http://d.hatena.ne.jp/ura_ra/20111026/1319642014
http://www.sist.ac.jp/~kanakubo/research/neuro.html
http://qiita.com/Ugo-Nama/items/04814a13c9ea84978a4c