C++で作るDeepLearning - パート2

# はじめに
この記事は[C++で作るDeepLearning - パート1](https://qiita.com/Luke02561/items/e68948da075893cff43d)のパート2です!

1. [C++で作るDeepLearning - パート1](https://qiita.com/Luke02561/items/e68948da075893cff43d)
2. C++で作るDeepLearning - パート2　←今回
3. [C++で作るDeepLearning - パート3](https://qiita.com/Luke02561/items/ee98d6e742febef04731)

パート1~3を全て実装したコードは[GitHub](https://github.com/Luke256/CppNN)にあります。(多少コードが違ってるかも)

このパートでは、ニューラルネットワークの学習を解説・実装していきます

# 学習の流れ
まずは、ニューラルネットワークの学習を解説していきましょう
学習には、**損失関数**というものを使います。これは、「予測した値と正解の値がどれだけ違うか」というものを求めるためのものです。つまり、ネットワークの学習というのは、この損失関数の値が小さくなるように重みやバイアスを調整していくことになります。

では、どのように重みやバイアスを調整すれば損失関数が(効率よく)小さくなるのかというと、「*そのパラメータをどのくらい変化させたら損失関数がどのくらい変化するのか*」というものを求めればよいです。これからは、それを「**勾配**」と呼びます。

つまり微分ですね

ただ、数学の微分のような、

```math
\frac{d}{dx}(x^2) = 2x
```

こんな微分はできません。なぜなら、数式がなかったり、とても複雑だからです。そこで、「そのパラメータ$P$の勾配」は、次のように求めます。

1. $P$の値を元の値から$h$だけ大きくする
2. そのときの損失関数$l_1$を求める
3. $P$の値を元の値から$h$だけ小さくする
4. そのときの損失関数$l_2$を求める
5. 最後に、$\frac{l_1-l_2}{2h}$を計算すれば$P$の勾配が求まる

勾配は損失関数を小さくするために最適な方向を示すため、これにしたがって微小量だけパラメータを更新すれば良いです。具体的には、

```math
w_{ij}=w_{ij} - grad_{ij} lr
```

こんな感じで更新します。$w_{ij}$は前のレイヤーのニューロンjから現在のニューロンjの重み、$grad_{ij}$はその勾配、$lr$は先ほどでいう「微小量」で、**学習率**と言います。
この学習率は、大きすぎると変化が大胆すぎて学習にならなかったり、小さすぎても学習がほとんど進まないので、適切な値である必要があります。

ちなみに、このようにして損失関数の値を減少させていく手法を「**確立的勾配降下法**(**SGD**)」といいます。他にも「Adam」や「Momentum」といった手法があります。

では実装に移っていきましょう

# 実装
## 損失関数
損失関数はいくつか種類がありますが、今回はそのうちの二つ「二乗和誤差」と「交差エントロピー誤差」を実装します。

### ミニバッチ学習
とその前に、ミニバッチ学習について少し解説しておきます。
学習時、学習用のデータが少ない時は、それら全てのデータを使って学習ができますが、通常はデータ数は10,000個以上だったりします。そうなると、学習1回ごとに全てのデータを使うと重くてしょうがありません。
そこで行われるのが**ミニバッチ学習**です。これは、1回の学習にデータ全てを使うのではなく、全データからいくつか選んで、それで学習するというものです。こうすることで大きいデータ数でも効率的に学習ができます。

### 二乗和誤差
二乗和誤差は、その名の通り、「教師データとの差を2乗して足し合わせた値」です。こうすれば、誤差が大きいほど損失関数の値も大きくなるし、2乗してるので絶対に値は正になります。
式はこんな感じ

```math
E = \frac{1}{2}\sum_k (y_k-t_k)^2
```

$y$が推論結果、$t$が教師データです。
これを、ミニバッチ学習に対応させてみるとこうなります

```math
E = \frac{1}{2}\frac{1}{N}\sum^N_{n=0}\sum_k (y_k-t_k)^2
```

真ん中にちょっと増えましたね。ここで重要なのは、ミニバッチ学習時の損失関数が渡すべき値は、*平均して1個当たりの損失である*ということです。(案外忘れやすい)

では実装していきましょう。前のパートと同じく、最初からバッチに対応させてしまいます。

```C++
long double mean_squared_error(std::vector<std::vector<long double>> &y, std::vector<std::vector<long double>> &t){
        long double res = 0;
        for(int i=0; i < y.size(); ++i){
            long double sum=0;
            for(int j=0; j<y[0].size(); ++j){
                sum += (y[i][j] - t[i][j]) * (y[i][j] - t[i][j]);
            }

            sum /= 2*y.size();
            res+=sum;
        }

        return res;
    }
```

### 交差エントロピー誤差
交差エントロピー誤差は以下のような式で表されます

```math
E = - \sum_kt_k\,log\,y_k
```
この誤差は(自然)対数の性質上、y_kが1の時この値は0となり、逆に0に近づけば近づくほど大きくなります。
さて、この誤差がどんな場面で使われるかというと、「**いくつか選択肢があって、その中から正解の一つを選ぶ**」問題でよく使われます。なぜなら、このような問題の教師データは正解ラベルが1、それ以外は0と表される(ont-hot表現)からです。このことがどう影響してくるかというと、この交差エントロピー誤差は、実質的に正解ラベルと対応する出力にのみ働くというわけです。
ちなみに、この損失関数は最後がsoftmax関数の時に相性がいいです

そして、これをバッチに対応させると式は下のようになります。


```math
E = -\frac{1}{N} \sum_n \sum_kt_{nk}\,log\,y_{nk}
```

二乗和誤差と同じく、バッチ内の誤差の平均になりました。では、実装していきましょう

```C++
    long double mean_cross_entropy_error(std::vector<std::vector<long double>> &y, std::vector<std::vector<long double>> &t){
        long double res = 0;
        for(int i=0; i < y.size(); ++i){
            long double sum=0;
            long double delta = 1e-7;
            for(int j=0; j<y[0].size(); ++j){
                // y[i][j]が0だとバグが起きるのでdeltaを加算して計算
                sum += t[i][j] * std::log(y[i][j] + delta);
            }
            sum *= -1;
            res+=sum;
        }
        return res / y.size();
    }
```

### まとめる
最後に、この二つの損失関数を一つの関数にまとめてしまいましょう。これは[前のパート](https://qiita.com/Luke02561/items/e68948da075893cff43d)で実装したModelクラスの新しいメンバ関数として実装します。

ここで、Modelクラスに新しくメンバ変数`m_loss`を追加しましょう。これは、学習時に使用する損失関数の種類を表します。

受け取った訓練データと教師データを元に、モデルの損失関数の値を求めます。

```C++
long double caluculate_loss(std::vector<std::vector<long double>>&batch_x, std::vector<std::vector<long double>>&batch_y){
    // cen:考査エントロピー誤差
    // mse:二乗和誤差
    std::vector<std::vector<long double>>t = predict(batch_x);
    if(m_loss == "cen") return mean_cross_entropy_error(t, batch_y);
    return mean_squared_error(t, batch_y);
}
```

## 数値微分
ここまで、学習の指標となる損失関数が実装できました。次は、それを使って、「そのパラメータをどのくらい変化させたら損失関数がどのくらい変化するのか」という勾配を求めるプログラムを実装していきましょう。
ここで、もう一度勾配を求める手順を確認しましょう

1. $P$の値を元の値から$h$だけ大きくする
2. そのときの損失関数$l_1$を求める
3. $P$の値を元の値から$h$だけ小さくする
4. そのときの損失関数$l_2$を求める
5. 最後に、$\frac{l_1-l_2}{2h}$を計算すれば$P$の勾配が求まる

(ちなみに、このような、小さな差分で微分を求めることを「**数値微分**」と言います)
この作業を全てのパラメータについて行えば良さそうです

### 重みの数値微分
まずは、重みについて微分する関数を組んでいきます。

この関数もModelクラスの新しいメンバ関数にします

~~この記事書いてる時に気づいたけど、数値微分のコード、もっと綺麗に効率よくかけたなぁ~~

```C++
std::vector<std::vector<long double>>numerical_gradient_layer(std::vector<std::vector<long double>>&batch_x, std::vector<std::vector<long double>>&batch_y, int index){
// index:微分対象のレイヤー
    
    long double h = 1e-4;
    std::vector<std::vector<long double>> grad(model[index].neuron.size(), std::vector<long double>(model[index].neuron[0].size()));
    for(int i = 0; i < grad.size(); ++i){
        for(int j = 0; j < grad[i].size(); ++j){
            long double tmp = model[index].neuron[i][j];

            model[index].neuron[i][j] = tmp + h; // 1
            long double fxh1 = caluculate_loss(batch_x, batch_y); // 2

            model[index].neuron[i][j] = tmp - h; // 3
            long double fxh2 = caluculate_loss(batch_x, batch_y); // 4

            grad[i][j] = (fxh1 - fxh2) / (2 * h); // 5
            model[index].neuron[i][j] = tmp;
        }
    }

    return grad;
}
```

### バイアスの微分
バイアスの方は重みの微分のループが一つ減った程度でほとんど変わりません。

```C++
std::vector<long double>numerical_gradient_bias(std::vector<std::vector<long double>>&batch_x, std::vector<std::vector<long double>>&batch_y, int index){
    long double h = 1e-4;
    std::vector<long double> grad(model[index].bias.size());
    for(int i = 0; i < grad.size(); ++i){
        long double tmp = model[index].bias[i];

        model[index].bias[i] = tmp + h; // 1
        long double fxh1 = caluculate_loss(batch_x, batch_y); // 2

        model[index].bias[i] = tmp - h; // 3
        long double fxh2 = caluculate_loss(batch_x, batch_y); // 4

        grad[i] = (fxh1 - fxh2) / (2 * h); // 5
        model[index].bias[i] = tmp;
    }

    return grad;
}
```

## モデルの学習
では、学習の最後のステップ、パラメータの更新を実装して、モデルの学習ができるようにしましょう
やることは

1. ミニバッチ作成
2. ミニバッチを元に勾配を求める
3. 求めた勾配を元にパラメータ更新

これを指定されたステップ数繰り返すだけです

```C++
    std::vector<long double> fit(int step, long double learning_rate, std::vector<std::vector<long double>> &x, std::vector<std::vector<long double>> &y, int batch_size, std::string loss)
    {

        std::vector<long double> history;
        m_loss = loss;

        for (int current_step = 0; current_step < step; ++current_step)
        {
            std::vector<std::vector<long double>> batch_x, batch_y;
            // diff_error.clear();
            // ミニバッチの作成
            for (int i = 0; i < batch_size; ++i)
            {
                batch_x.push_back(x[(batch_size * current_step + i) % x.size()]);
                batch_y.push_back(y[(batch_size * current_step + i) % y.size()]);
            }

            // // 逆伝播の準備
            // auto test_y = predict(batch_x);
            // for (int i = 0; i < test_y.size(); ++i)
            // {
            //     std::vector<long double> t;
            //     for (int j = 0; j < test_y[i].size(); ++j)
            //     {
            //         t.push_back((test_y[i][j] - batch_y[i][j]) / batch_size);
            //     }
            //     diff_error.push_back(t);
            // }

            // // 逆伝播
            // backward();

            // long double loss_step;
            // if (m_loss == "mse")
            //     loss_step = Loss::mean_squared_error(test_y, y);
            // else
            //     loss_step = Loss::mean_cross_entropy_error(test_y, y);

            for (int index = 0; auto &layer : model)
            {
                // 勾配を計算
                // 数値微分(遅い)
                std::vector<std::vector<long double>>layer_grad = numerical_gradient_layer(batch_x, batch_y, index);
                std::vector<long double>bias_grad = numerical_gradient_bias(batch_x, batch_y, index);

                // // 誤差逆伝播法(早い)
                // std::vector<std::vector<long double>> layer_grad = layer.grad_layer;
                // std::vector<long double> bias_grad = layer.grad_bias;

                // 計算した勾配に従って重みを調整
                // 1.layerの調整
                for (int i = 0; i < layer_grad.size(); ++i)
                {
                    for (int j = 0; j < layer_grad[i].size(); ++j)
                    {
                        layer.neuron[i][j] -= learning_rate * layer_grad[i][j];
                    }
                }

                // 2.biasの調整
                for (int i = 0; i < bias_grad.size(); ++i)
                {
                    layer.bias[i] -= learning_rate * bias_grad[i];
                }

                ++index;
            }

            // 学習経過の記録
            // history.push_back(loss_step);
            history.push_back(caluculate_loss(batch_x, batch_y));
        }

        return history;
    }
```

プログラム内で、コメントアウトされたソースコードがたくさんあります。これは、パート3の誤差逆伝播法の時に使うもので、今回は必要ありません

# まとめ
今回は、パラメータの勾配を求めることが分かってしまえばあとは比較的簡単な内容だったと思います
次回はニューラルネットワーク学習の要、誤差逆伝播法について解説していきます。
(あまりうまく動作するとは言い難いけど、~~そこは試行回数でごり押せばいい。速いし。~~改善しないとなぁ)

# ソースコード(おまけ)
最後に、今回実装したコードとおまけ程度にmain関数も置いておきます。
今回実装したネットワークにXORを計算させてます。

また、パート一と違うところには`// NEW!`というコメントをつけてます。
(パート1のソースコードと多少違う点があるので、コピペよりはNEWの部分を追記/削除することをお勧めします)

```C++
#include <bits/stdc++.h>

// NEW!
long double mean_squared_error(std::vector<std::vector<long double>> &y, std::vector<std::vector<long double>> &t)
{
    long double res = 0;
    for (int i = 0; i < y.size(); ++i)
    {
        long double sum = 0;
        for (int j = 0; j < y[0].size(); ++j)
        {
            sum += (y[i][j] - t[i][j]) * (y[i][j] - t[i][j]);
        }

        sum /= 2 * y.size();
        res += sum;
    }

    return res;
}

// NEW!
long double mean_cross_entropy_error(std::vector<std::vector<long double>> &y, std::vector<std::vector<long double>> &t)
{
    long double res = 0;
    for (int i = 0; i < y.size(); ++i)
    {
        long double sum = 0;
        long double delta = 1e-7;
        for (int j = 0; j < y[0].size(); ++j)
        {
            // y[i][j]が0だとバグが起きるのでdeltaを加算して計算
            sum += t[i][j] * std::log(y[i][j] + delta);
        }
        sum *= -1;
        res += sum;
    }
    return res / y.size();
}

class Activation
{
    /* sigmoid関数の定義 */
    std::vector<std::vector<long double>> sigmoid(std::vector<std::vector<long double>> &x)
    {
        std::vector<std::vector<long double>> res;

        for (auto batch : x)
        {
            // バッチごとに処理
            std::vector<long double> t;
            for (auto i : batch)
            {
                t.push_back(1.0 / (1.0 + exp(-i)));
            }
            res.push_back(t);
        }

        return res;
    }

    /* relu関数の定義 */
    std::vector<std::vector<long double>> relu(std::vector<std::vector<long double>> &x)
    {
        std::vector<std::vector<long double>> res;

        for (auto batch : x)
        {
            // バッチごとに処理
            std::vector<long double> t;
            for (long double i : batch)
            {
                t.push_back((i >= 0) * i);
            }
            res.push_back(t);
        }
        return res;
    }

    /* softmax関数の定義 */
    std::vector<std::vector<long double>> softmax(std::vector<std::vector<long double>> &x)
    {
        std::vector<std::vector<long double>> res;

        for (auto batch : x)
        {
            // バッチごとに処理
            std::vector<long double> t;
            long double c = *max_element(batch.begin(), batch.end());
            long double sum = 0;
            for (long double i : batch)
            {
                sum += exp(i - c);
            }
            for (long double i : batch)
            {
                t.push_back(exp(i - c) / sum);
            }
            res.push_back(t);
        }
        return res;
    }

    std::string m_name;

public:
    Activation();
    Activation(std::string name) : m_name(name) {}
    /* forward関数の定義 */
    std::vector<std::vector<long double>> forward(std::vector<std::vector<long double>> &x)
    {
        if (m_name == "sigmoid")
            return sigmoid(x);
        else if (m_name == "relu")
            return relu(x);
        else if (m_name == "softmax")
            return softmax(x);
        else
            return x;
    }
};

class Dense
{
public: // public:の位置変更 // NEW!
    Activation m_activation;
    std::vector<long double> bias;
    std::vector<std::vector<long double>> neuron;

    // コンストラクタ
    Dense(int input_unit, int unit, std::string activation) : bias(unit),
                                                              neuron(unit, std::vector<long double>(input_unit)),
                                                              m_activation(activation)
    {

        double sigma = 0.05;
        if (activation == "relu")
            sigma = std::sqrt(2.0 / (double)input_unit);
        else if (activation == "sigmoid" || activation == "leaner")
            sigma = std::sqrt(1.0 / (double)input_unit);
        else
            sigma = 0.05;

        // 重みとバイアスを初期化
        std::random_device seed;
        std::mt19937 engine(seed());
        std::normal_distribution<> generator(0.0, sigma);

        for (int i = 0; i < unit; ++i)
        {
            bias[i] = generator(engine);
            for (int j = 0; j < input_unit; ++j)
            {
                neuron[i][j] = generator(engine);
            }
        }
    }

    // forward関数
    std::vector<std::vector<long double>> forward(std::vector<std::vector<long double>> &data)
    {

        std::vector<std::vector<long double>> ans;

        // データごとに処理
        for (int index = 0; auto &i : data)
        {

            std::vector<long double> res;
            for (int j = 0; j < neuron.size(); ++j)
            {
                long double t = 0;
                // 入力 * 重み
                for (int k = 0; k < neuron[j].size(); ++k)
                {
                    t += i[k] * neuron[j][k];
                }

                // バイアスの適用
                t -= bias[j];

                res.push_back(t);
            }

            ans.push_back(res);

            ++index;
        }

        ans = m_activation.forward(ans); // 活性化関数を適用

        return ans;
    }
};

class Model
{
private:
    int m_input_size, m_output_size;
    std::vector<Dense> model;

    std::string m_loss; // NEW!

    // NEW!
    long double caluculate_loss(std::vector<std::vector<long double>> &batch_x, std::vector<std::vector<long double>> &batch_y)
    {
        // cen:考査エントロピー誤差
        // mse:二乗和誤差
        std::vector<std::vector<long double>> t = predict(batch_x);
        if (m_loss == "cen")
            return mean_cross_entropy_error(t, batch_y);
        return mean_squared_error(t, batch_y);
    }

    // NEW!
    std::vector<std::vector<long double>> numerical_gradient_layer(std::vector<std::vector<long double>> &batch_x, std::vector<std::vector<long double>> &batch_y, int index)
    {
        // index:微分対象のレイヤー

        long double h = 1e-4;
        std::vector<std::vector<long double>> grad(model[index].neuron.size(), std::vector<long double>(model[index].neuron[0].size()));
        for (int i = 0; i < grad.size(); ++i)
        {
            for (int j = 0; j < grad[i].size(); ++j)
            {
                long double tmp = model[index].neuron[i][j];

                model[index].neuron[i][j] = tmp + h;                  // 1
                long double fxh1 = caluculate_loss(batch_x, batch_y); // 2

                model[index].neuron[i][j] = tmp - h;                  // 3
                long double fxh2 = caluculate_loss(batch_x, batch_y); // 4

                grad[i][j] = (fxh1 - fxh2) / (2 * h); // 5
                model[index].neuron[i][j] = tmp;
            }
        }

        return grad;
    }

    // NEW!
    std::vector<long double> numerical_gradient_bias(std::vector<std::vector<long double>> &batch_x, std::vector<std::vector<long double>> &batch_y, int index)
    {
        long double h = 1e-4;
        std::vector<long double> grad(model[index].bias.size());
        for (int i = 0; i < grad.size(); ++i)
        {
            long double tmp = model[index].bias[i];

            model[index].bias[i] = tmp + h;                       // 1
            long double fxh1 = caluculate_loss(batch_x, batch_y); // 2

            model[index].bias[i] = tmp - h;                       // 3
            long double fxh2 = caluculate_loss(batch_x, batch_y); // 4

            grad[i] = (fxh1 - fxh2) / (2 * h); // 5
            model[index].bias[i] = tmp;
        }

        return grad;
    }

public:
    Model(int input_size) : m_input_size(input_size),
                            m_output_size(input_size)
    {
    }

    void AddDenseLayer(int unit, std::string activation)
    {
        Dense dense(m_output_size, unit, activation);
        model.push_back(dense);
        m_output_size = unit;
    }
    std::vector<std::vector<long double>> predict(const std::vector<std::vector<long double>> &data)
    {
        std::vector<std::vector<long double>> res = data;
        for (auto &layer : model)
        {
            res = layer.forward(res);
        }

        return res;
    }

    // NEW!
    std::vector<long double> fit(int step, long double learning_rate, std::vector<std::vector<long double>> &x, std::vector<std::vector<long double>> &y, int batch_size, std::string loss)
    {

        std::vector<long double> history;
        m_loss = loss;

        for (int current_step = 0; current_step < step; ++current_step)
        {
            std::vector<std::vector<long double>> batch_x, batch_y;
            // diff_error.clear();
            for (int i = 0; i < batch_size; ++i)
            {
                batch_x.push_back(x[(batch_size * current_step + i) % x.size()]);
                batch_y.push_back(y[(batch_size * current_step + i) % y.size()]);
            }

            // // 逆伝播の準備
            // auto test_y = predict(batch_x);
            // for (int i = 0; i < test_y.size(); ++i)
            // {
            //     std::vector<long double> t;
            //     for (int j = 0; j < test_y[i].size(); ++j)
            //     {
            //         t.push_back((test_y[i][j] - batch_y[i][j]) / batch_size);
            //     }
            //     diff_error.push_back(t);
            // }

            // // 逆伝播
            // backward();

            // long double loss_step;
            // if (m_loss == "mse")
            //     loss_step = Loss::mean_squared_error(test_y, y);
            // else
            //     loss_step = Loss::mean_cross_entropy_error(test_y, y);

            for (int index = 0; auto &layer : model)
            {
                // 勾配を計算
                // 数値微分(遅い)
                std::vector<std::vector<long double>>layer_grad = numerical_gradient_layer(batch_x, batch_y, index);
                std::vector<long double>bias_grad = numerical_gradient_bias(batch_x, batch_y, index);

                // // 誤差逆伝播法(早い)
                // std::vector<std::vector<long double>> layer_grad = layer.grad_layer;
                // std::vector<long double> bias_grad = layer.grad_bias;

                // 計算した勾配に従って重みを調整
                // 1.layerの調整
                for (int i = 0; i < layer_grad.size(); ++i)
                {
                    for (int j = 0; j < layer_grad[i].size(); ++j)
                    {
                        layer.neuron[i][j] -= learning_rate * layer_grad[i][j];
                    }
                }

                // 2.biasの調整
                for (int i = 0; i < bias_grad.size(); ++i)
                {
                    layer.bias[i] -= learning_rate * bias_grad[i];
                }

                ++index;
            }

            // 学習経過の記録
            // history.push_back(loss_step);
            history.push_back(caluculate_loss(batch_x, batch_y));
        }

        return history;
    }
};

int main()
{
    Model model(2);

    model.AddDenseLayer(10, "relu");
    model.AddDenseLayer(2, "softmax");

    std::vector<std::vector<long double>> x = {
        {0, 1},
        {0, 0},
        {1, 1},
        {1, 0}
    };
    
    std::vector<std::vector<long double>> y = {
        {0, 1},
        {1, 0},
        {1, 0},
        {0, 1}
    };

    auto t = model.predict(x);

    for(auto i : t){
        for(auto j : i) std::cout << j << " ";
        std::cout << std::endl;
    }

    auto history = model.fit(1000, 0.1, x, y, 4, "cen");

    t = model.predict(x);

    std::cout << std::endl;

    for(auto i : t){
        for(auto j : i) std::cout << j << " ";
        std::cout << std::endl;
    }
}
```