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

# はじめに
ディープラーニングなどの機械学習には、多くの場合強力なライブラリ「numpy」を持つPythonが使用されます。ただPythonには一つ欠点が...

**実行速度が遅い**

そこで僕はこう考えました。「じゃあ実行速度の速いC++で作れないものか」と。
というわけで、オライリー・ジャパンさんの『[ゼロから作るDeep Learning](https://www.oreilly.co.jp/books/9784873117584/)』を参考に実際にC++で作ってみたので記事にしてみました

一応、この記事内でもディープラーニングの仕組みの解説はしていきますが、ある程度ディープラーニングの仕組みがわかっていた方が読みやすいと思います。

追記(2023/2/25)
~~正直Pythonでやった方が速い~~

# この記事について
この記事は、3つのパートに分けて書こうと思います。

- パート1：ネットワークの構築・推論
- [パート2](https://qiita.com/Luke02561/items/74995d9d5f0e20cbf2a0)：数値微分を使った学習
- [パート3](https://qiita.com/Luke02561/items/ee98d6e742febef04731)：誤差逆伝播法を使った学習

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

追記(2023/10/25)
新しく実装しなおしたものが[「NNCpp」(GitHub)](https://github.com/Luke256/NNCpp)にあります。こちらの方は比較的使いやすいようになっていると思うので、もし、使う場合はこちらを使ってください。

# ニューラルネットワークの推論の方法
## ニューロン間の重み
実装をしていく前に、軽くニューラルネットワークがどうやって推論・予測をするのかについて解説します。

ニューラルネットワークは、ざっと下のような構造になっています
![アートボード 1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/596925/8b2e6e6f-ba68-67a5-c82d-0a6aac0a6d00.png)

入力層と中間層に着目すると、まず入力層に入ったデータは、水色の線を通る時にその「**重み**」が乗算されて中間層に足されます(ちなみに、図の中の丸は「**ニューロン**」と呼ばれます)。つまり、i番目の入力データを$x_i$、中間層のi番目のニューロンに来るデータを$y_i$、j番目の入力データから中間層のi番目のニューロンの間の重みを$w_{ij}$とすると、

```math
y_i = \sum^{n}_{j=0}y_j w_{ij}
```

となります。言葉にすると、「入力値とそれに対応した重みを掛け合わせたものの和」と言ったところでしょうか。

## バイアス
図にはもう一つ、黄緑色の線があります。これは「**バイアス**」というもので、一般にそのニューロンの発火のしやすさを表し、赤丸の部分からは常に1が出力されます。つまりそのバイアスの重みがそのまま、つながっているニューロンに届くわけです。  
そしてこのバイアスは、「入力値とそれに対応した重みを掛け合わせたものの和」と足し合わされる(引かれる)ことになります。中間層のi番目のニューロンのバイアスを$b_i$と定義すると、先程の式は最終的にこのようになります。

```math
y_i = \sum^{n}_{j=0}y_j w_{ij} - b_i
```

## 活性化関数
さて、これだけでは終わりません。ニューロンは、前の層のニューロンから重みが乗算された数値に自身のバイアスが足された数をそのまま次のニューロンに渡すのではなく、**活性化関数**というものを通してから次のニューロンに渡します。  
活性化関数、なんか名前からしていかにも重要そうですよね、そうですよね!  
実際、活性化関数はニューラルネットワーク、特にニューラルネットワークをたくさん重ねたディープラーニングでとても大きな役割を担います。
もう少し丁寧にいうと、活性化関数はレイヤーそれぞれが持っていて、受け取った数値をどう活性化させるか、どのように発火するかというのを決定させる役割があります。

## 小まとめ
ここまでをまとめると、ニューラルネットワークの推論は、

1. 出力した数値を対応した重みと掛け合わせて・足し合わせて次の層に渡す
2. 前の層から数値を受け取った層は、それに自身のバイアスを足し合わせる
3. 最後に、その値を自身の活性化関数に渡す
4. 1.に戻る

これの繰り返しになります。

# 実装

## 活性化関数
まずは活性化関数から実装していきます。方針としては、クラスのインスタンス化時に活性化関数の種類を文字列として受け取り、forward関数が呼び出されるたびにその受け取った文字列から適切なメンバ関数を呼び出す形にします。

あとのことを考えて、最初からバッチ対応にしましょう

### メンバ変数

| 型 | 変数名 | 意味・役割 |
| --- | --- | --- |
| ActivationType | m_name | 活性化関数の種類 |
※ 2021/8/31 追記
m_nameの型を`std::string`から`ActivationType`に変更しました
ActivationTypeは以下のように定義されています

```C++
enum ActivationType{
    Sigmoid,
    Linear,
    SoftMax,
    Relu
};
```

### コンストラクタ
このクラスでのコンストラクタは、使用する活性化関数の種類を表す`name`を引数として受け取って、それをメンバ変数の`m_name`に代入するだけです。(ここでの実装は省きます)

### sigmoid関数
sigmoid関数は、以下のような式で表されます

```math
h(x)=\frac{1}{1+exp(-x)}
```

expというのは、ネイピア数の$e$を底とした指数関数です。これを使って、sigmoid関数を実装しましょう

```C++
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関数
ReLU関数はsigmoid関数よりは単純で、以下のように表されます

```math
h(x)=\left\{
\begin{array}{ll}
x & (x≥0) \\
0 & (x<0)
\end{array}
\right.
```
では実装していきましょう

```C++
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関数
次はsoftmax関数です。この関数は主に出力層に使われます。式はこんな感じ

```math
y_k=\frac{exp(a_k)}{\sum^n_{i=1} exp(a_i)}
```
ちょっとわかりにくいかもしれないですが、これは「出力された数値のうち、どのくらいの大きさを占めているか」のようなものを表します。したがって、この関数は確率的な推論(1である確率は23%、2である確率は46%...のような推論)をする場合によく用いられます。  
　ただ、これだとオーバーフローするかもしれないため、通常はこれを変形して以下のような式を使うらしいです

```math
y_k=\frac{exp(a_k-C')}{\sum^n_{i=1} exp(a_i-C')}
```
ここでいう$C'$は入力された信号の最大値が一般的みたいです(マイナスのオーバーフローってしないのかな)  
では実装していきましょう

```C++
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;
}
```

### forward関数
そして、これら3つを実行するためのforward関数を実装します。

```C++
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;
}
```

### Activationクラス(まとめ)
最後に、今までの4関数をまとめたActivationクラスを作成します(全て書くと長くなるので多少省略しています)

```C++
class Activation
{
    /* sigmoid関数の定義 */

    /* relu関数の定義 */

    /* softmax関数の定義 */

    ActivationType m_name;

public:
    Activation();
    Activation(ActivationType name) : m_name(name) {}
    /* forward関数の定義 */
};
```

## レイヤー
さて、次にネットワークのレイヤー1層分となるレイヤークラスを設計していきます

### メンバ変数

| 型 | 変数名 | 意味・役割 |
| --- | --- | --- |
| Activation | m_activation | 活性化関数|
| std::vector\<long double> | bias | バイアスの重み |
| std::vector\<std::vector\<long double>> | neuron | ニュローンの重み |

`activation`はさっき実装したActivationクラスです  
`bias`と`neuron`の添字は以下のように定めます

- bias[i] := ニューロン$i$のバイアス
- neuron[i][j] := 前の層のニューロン$j$から現在の層のニューロン$i$の重み

`neuron`の方は添字が少しややこしいかもしれませんが、こうした方が都合がいいのでこうしてます。

### コンストラクタ
#### 引数
| 型 | 変数名 | 意味・役割 |
| --- | --- | --- |
| int | input_unit | 入力数 |
| int | unit | 出力数 |
| ActivationType | activation | 使用する活性化関数の種類 |

#### 処理
ここで行うのは、レイヤーの重みとバイアスを表す配列の初期化です。これらは、平均0、標準偏差$\sigma$の正規分布の乱数で初期化します

さて、この$\sigma$ですが、使用される活性化関数によって適した数値が変わってきます。

| 活性化関数 | 適した$\sigma$ |
|---|---|
|sigmoid | $\frac{1}{\sqrt{n}}$
|relu| $\sqrt{\frac{2}{n}}$

このほかの活性化関数について適切な$\sigma$は僕は知らないので、とりあえず0.05を使うことにします。

#### 実装

```C++
Dense(int input_unit, int unit, Activation activation):
    bias(unit),
    neuron(unit, std::vector<long double>(input_unit)),
    m_activation(activation)
    {

        double sigma = 0.05;
        if(activation == ActivationType::Relu) sigma = std::sqrt(2.0 / (double)input_unit);
        else if(activation == ActivationType::Sigmoid || activation == ActivationType::Softmax) 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関数
次に、forward関数を実装していきます。これは、推論時に前のレイヤーから受け取った数値を処理して、次のレイヤーに渡す(**順伝播**させる)関数です

#### 仕組み
では、もう一度順伝搬の仕組みについて見ていきます。ここでは、下の図のOutput-1について見てみます。　　
Output-1が受け取る数値は、「Input-1の値と、Input-1とOutput-1の重みをかけあわせたもの」、「Input-2の値と、Input-2とOutput-1の重みをかけあわせたもの」、そして自身のバイアスです。式にするとこんな感じです

```math
Out_i = \Sigma^{input}_{j=0}Input_j w_{ij} - Bias_i
```

計算したら、最後にこれを自身の活性化関数に通して処理は完了です

![アートボード 2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/596925/74bec3a1-04b0-c752-9771-6d4faf452de3.png)

#### 実装
```C++
    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;
    }
```

### Denseクラス(まとめ)
では、ここまでをまとめて、Denseクラスを実装していきましょう。(関数などはもちろん省略)
ちなみに、Denseというのは「全結合層」という意味(意訳)です。レイヤーには他の種類もあるので、それらとの差別化のためにこういった名前にしています

```C++
class Dense{
    Activation m_activation;
    std::vector<long double>bias;
    std::vector<std::vector<long double>> neuron;

public:
    // コンストラクタ

    // forward関数
};
```

## モデル
さて、ようやく材料が揃ったのでネットワークの本体となるモデルクラスを実装していきましょう

### メンバ変数
|型|変数名|意味・役割|
|---|---|---|
|int | m_input_size | 入力サイズ|
|int | m_output_size | 出力サイズ |
|std::vector<Dense>|model|ネットワーク本体|

### コンストラクタ
インスタンスかの時点では何も特別なことはしません。  
入力サイズを受け取って、`m_input_size`と`m_output_size`に代入するだけです。(最初はレイヤーはないので入力=出力になります)

```C++
Model(int input_size):
m_input_size(input_size),
m_output_size(input_size)
{}
```

### 層の追加
層の追加をしていきます。
ここは特にこれといった説明はないです

```C++
void AddDenseLayer(int unit, ActivationType activation){
    Dense dense(m_output_size, unit, activation);
    model.push_back(dense);
    m_output_size = unit;
}
```

### 推論
ようやくメインの推論に入ります。といっても、推論処理のほとんどはレイヤーの時に実装できているので、ここでやるのはレイヤー同士の橋渡し程度です

```C++
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;
}
```

これだけ。とっても簡単。(難しいのはパート2から...)

### Modelクラス(まとめ)

ではModelクラスも実装していきましょう(例のごとく関数などは省略)

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

public:
    /*コンストラクタ*/

    /*層の追加*/

    /*推論*/
};
```

## まとめ
ここまで、ニューラルネットワークの構築と推論の実装をしてきました。次回はこのネットワークの学習を実装していきましょう

## ソースコード(おまけ)
最後に、今回実装したコードとおまけ程度にmain関数も置いておきます。現時点では学習とかもないのでエラーが出てなければ多分正常に動いてます。

bitsを使ってますが、嫌いな方は各自で置き換えてください


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

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
{
    Activation m_activation;
    std::vector<long double> bias;
    std::vector<std::vector<long double>> neuron;

public:
    // コンストラクタ
    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;

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;
    }
};

int main(){
    Model model(2);

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

    std::vector<std::vector<long double>>x={
        {1, 2},
        {3, 4},
        {5, 6}
    };

    auto y = model.predict(x);

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