C言語でニューラルネットワークの実装（５）〜モデルの構造と順伝播の実装〜

# 前回まで
[前回](https://qiita.com/tky1117/items/7cdebb370e0da7ba6f71)はデータセットを準備しました。今回は，モデルの構造と順伝播を実装していきます。

# モデルの構造

モデルの構造についての情報を表すための`MODEL_PARAMETER`型を定義します。

```C:_struct.h
typedef struct {
    int L; //number of hidden layers
    int K; //number of classes
    int dim_in; //input dimension
    int *N_unit; //number of units
    float *(**activation)(float *a, int D); // activation function
    float (*loss)(float *y, float *t, int K); //loss function
} MODEL_PARAMETER;
```

`MODEL_PARAMETER` 型では，
1. 隠れ層の数（ `L` ）
2. 分類するクラスの数（ `K` ）
3. 入力の次元（ `dim_in` ）
4. 各層のユニットの数（ `N_unit` ）
5. 各層の活性化関数（ `activation` ）
6. 損失関数（ `loss` ）
を持っています。特に**活性化関数と損失関数は，条件分岐を少なくするために関数ポインタで持ちます**。次の `build_model` 関数では，これらのパラメータを決めていきます。隠れ層の数（ `L` ），各層のユニットの数（ `N_unit[l]` ），出力層以外の活性化関数（ `activation` ）はいくらか変更の余地があり，この値を変えることでモデルの構造が変化します。

```C:forward_sample.c
#define DATA_LINE 5

MODEL_PARAMETER build_model(void){
    int l;
    MODEL_PARAMETER model_parameter;
    
    model_parameter.L = 2;
    model_parameter.K = 3;
    model_parameter.dim_in = DATA_LINE * DATA_LINE;
    
    model_parameter = setup_model_parameter(model_parameter);
    
    for(l = 1; l <= model_parameter.L; l++) {
        model_parameter.N_unit[l] = 16;
    }
    
    for(l = 1; l <= model_parameter.L; l++) {
        model_parameter.activation[l] = ReLU;
    }
    
    model_parameter.activation[model_parameter.L + 1] = softmax;
    model_parameter.loss = cross_entropy;
    
    setup_parameters(model_parameter);
    
    return model_parameter;
}
```

なお， `setup_parameters(model_parameter);` は，モデルの重み $w$ や，逆伝播で重要な $\partial E / \partial a$ の配列の領域を確保する関数です。

# 活性化関数の再実装
順伝播の実装に入る前に，活性化関数についてもう一度考えたいと思います。[第1回](https://qiita.com/tky1117/items/ffc367e52c4a04cffa1d)で実装した活性化関数はどれも $\boldsymbol{a}$ を引数にとり， $\boldsymbol{y}$ を返す関数でした。しかし，逆伝播ではどれも $\boldsymbol{a}$ を引数にとるものの，ソフトマックス関数は行列を，それ以外の関数はベクトルを返します。これでは扱いづらいので**全て行列を返す**関数に変更します。他の活性化関数もソフトマックス関数のように

```math
\frac{\partial z_{d}}{\partial a_{d'}} = \begin{cases}
    h'(a_{d})  & (d = d')\\
    0 & (d \neq d')
\end{cases}
```

とすることで行列のように扱うことができます。
また，活性化関数の引数に順伝播のときの通常の計算か，それとも逆伝播のときの微分の計算かを指定する引数 `derivative` を与えます。`derivative == 1` であれば逆伝播を意味します。
逆伝播の結果 $\partial \boldsymbol{z} / \partial \boldsymbol{a}$ は， 二次元配列の`derivative_matrix` に書き込みます。この配列 `derivative_matrix` の領域の確保は関数の外でしておきます。

```math
\frac{\partial \boldsymbol{z}}{\partial \boldsymbol{a}} = \left(\begin{matrix}
    \dfrac{\partial z_{1}}{\partial a_{1}} & \dfrac{\partial z_{1}}{\partial a_{2}} & \cdots & \dfrac{\partial z_{1}}{\partial a_{d'}} & \cdots & \dfrac{\partial z_{1}}{\partial a_{D}}\\
    \dfrac{\partial z_{2}}{\partial a_{1}} & \dfrac{\partial z_{2}}{\partial a_{2}} & \cdots & \dfrac{\partial z_{2}}{\partial a_{d'}} & \cdots & \dfrac{\partial z_{2}}{\partial a_{D}} \\
    \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\
    \dfrac{\partial z_{d}}{\partial a_{1}} & \dfrac{\partial z_{d}}{\partial a_{2}} & \cdots & \dfrac{\partial z_{d}}{\partial a_{d'}} & \cdots  & \dfrac{\partial z_{d}}{\partial a_{D}} \\
    \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\
    \dfrac{\partial z_{D}}{\partial a_{1}} & \dfrac{\partial z_{D}}{\partial a_{2}} & \cdots & \dfrac{\partial z_{D}}{\partial a_{d'}} & \cdots & \dfrac{\partial z_{D}}{\partial a_{D}}
\end{matrix}\right) \in \mathbb{R}^{D \times D}
```

つまり，各活性化関数の名前で名付けられた関数は，

1. 順伝播の場合（ `derivative == 0` ）は， $\boldsymbol{z} = h(\boldsymbol{a})$ を返し， `derivative_matrix` には何も書き込まず，
2. 逆伝播の場合（ `derivative == 1` ）は， $\partial\boldsymbol{z} / \partial \boldsymbol{a}$ を， `derivative_matrix` に書き込み， `NULL` を返す

というような実装に変更します。

## 1. ステップ関数
ステップ関数の微分は次のように表されます。

```math
\frac{\partial z_{d}}{\partial a_{d'}} = 0
```

[第1回](https://qiita.com/tky1117/items/ffc367e52c4a04cffa1d)のステップ関数の実装に，偏微分の要素を追加していきます。
※ `z` の領域確保は省略しているので注意。

```C:activation_function.h
float *step(float *a, int D, int derivative, float **derivative_matrix) {
    int d, d_dash;
    float *z = NULL;
    
    if(derivative) {
        for(d = 1; d <= D; d++) {
            for(d_dash = 1; d_dash <= D; d_dash++) {
                derivative_matrix[d][d_dash] = 0.0f;
            }
        }
        return NULL;
    }
    
    z[0] = 1.0f;
    
    for(d = 1; d <= D; d++) {
        if(a[d] < 0.0f) {
            z[d] = 0.0f;
        } else {
            z[d] = 1.0f;
        }
    }
    
    return z;
}
```

## 2. シグモイド関数
シグモイド関数の微分は次のように表されます。

```math
\frac{\partial z_{d}}{\partial a_{d'}} = \begin{cases}
    h(a_{d})(1 - h(a_{d})) & (d = d')\\
    0 & (d \neq d')
\end{cases}
```

```C:activation_function.h
float *sigmoid(float *a, int D, int derivative, float **derivative_matrix) {
    int d, d_dash;
    float *z = NULL;
    z[0] = 1.0f;
    
    for(d = 1; d <= D; d++) {
        z[d] = 1.0f / (1.0f + expf(- a[d]));
    }
    
    if(derivative) {
        for(d = 1; d <= D; d++) {
            for(d_dash = 1; d_dash <= D; d_dash++) {
                derivative_matrix[d][d_dash] = 0.0f;
            }
        }
        
        for(d = 1; d <= D; d++) {
            derivative_matrix[d][d] = z[d] * (1.0f - z[d]);
        }
        
        return NULL;
    }
    
    return z;
}
```

## 3. ReLU関数
ReLU関数の微分は次のように表されます。

```math
\frac{\partial z_{d}}{\partial a_{d'}} = \begin{cases}
    0 & (d \neq d') \\
    0 & (d = d' \mathrm{かつ} a_{d'} < 0) \\
    1 & (d = d' \mathrm{かつ} a_{d'} \geqq 0)
\end{cases}
```

```C:activation_function.h
float *ReLU(float *a, int D, int derivative, float **derivative_matrix) {
    int d, d_dash;
    float *z = NULL;

    if(derivative) {
        for(d = 1; d <= D; d++) {
            for(d_dash = 1; d_dash <= D; d_dash++) {
                derivative_matrix[d][d_dash] = 0.0f;
            }
        }
        
        for(d = 1; d <= D; d++) {
            if(a[d] > 0.0f) {
                derivative_matrix[d][d] = 1.0f;
            }
        }
        
        return NULL;
    }
    
    z[0] = 1.0f;
    
    for(d = 1; d <= D; d++) {
        if(a[d] < 0.0f) {
            z[d] = 0.0f;
        } else {
            z[d] = a[d];
        }
    }
    
    return z;
}
```

## 4. ソフトマックス関数

ソフトマックス関数は指数関数 $e^{x}$ を扱うため，オーバーフローが起きる可能性があります。これを防ぐために， $a_{k}$ から定数 $C$ を引くという方法があります[^1]。 $C > 0$ ならば， $e^{a_{k}} > e^{a_{k} - C}$
であるため，オーバーフローを防ぐことができます。

```math
y_{k} = \frac{e^{a_{k} - C}}{{\displaystyle \sum_{k'=1}^{K}}e^{a_{k'} - C}}
```

この定数 $C$ は，プログラムでは

```math
C = \max_{k}a_{k}
```

としています。

```C:activation_function.h
float *softmax(float *a, int K, int derivative, float **derivative_matrix) {
    int k, k_dash;
    float *y = NULL;
    float sum_exp = 0.0f;
    float max_a = a[1];

    for(k = 2; k <= K; k++) {
        if(a[k] > max_a) {
            max_a = a[k];
        }
    }
    
    y[0] = 1.0f;
    
    for(k = 1; k <= K; k++) {
        y[k] = expf(a[k] - max_a);
        sum_exp += y[k];
    }
    
    for(k = 1; k <= K; k++) {
        y[k] /= sum_exp;
    }
    
    if(derivative) {
        for(k = 1; k <= K; k++) {
            for(k_dash = 1; k_dash < k; k_dash++) {
                derivative_matrix[k][k_dash] = -y[k] * y[k_dash];
                derivative_matrix[k_dash][k] = -y[k] * y[k_dash];
            }
            derivative_matrix[k][k] = y[k] * (1.0f - y[k]);
        }
        
        return NULL;
    }
    
    return y;
}
```

# 損失関数の再実装
損失関数も，活性化関数と同様に逆伝播の際の微分について考えます。こちらも，引数の `derivative` で順伝播に利用されるのか，それとも逆伝播の微分として利用されるのかを指定していきます。
逆伝播の微分の場合，引数で受け取る `dE_dy` に計算結果を書き込み， `0` を返します。

## 1. 平均二乗誤差
平均二乗誤差は，

```math
    E = \frac{1}{D}\sum_{d=1}^{D}(t_{d} - y_{d})^{2}
```

で，その微分は，

```math
    \frac{\partial E}{\partial y_{d}} = - \frac{2(t_{d} - y_{d})}{D}
```

でした。

```C:loss_function.h
float mean_squared_error(float *y, float *t, int D, int derivative, float *dE_dy) {
    int d;
    float delta;
    float E = 0.0f;

    if(derivative) {
        for(d = 1; d <= D; d++) {
            dE_dy[d] = - 2.0f * (t[d] - y[d]) / (float) D;
        }
        return 0.0f;
    }

    for(d = 1; d <= D; d++) {
        delta = t[d] - y[d];
        E += delta * delta;
    }
    
    return E / (float) D;
}
```

## 2. 交差エントロピー
交差エントロピーは，

```math
    E = -\sum_{k=1}^{K}t_{k}\log y_{k}
```

で，その微分は，

```math
    \frac{\partial E}{\partial y_{k}} = - \frac{t_{k}}{y_{k}}
```

でした。

```C:loss_function.h
float EPS = 1e-12;

float cross_entropy(float *y, float *t, int K, int derivative, float *dE_dy) {
    int k;
    float E = 0.0f;
    
    if(derivative) {
        for(k = 1; k <= K; k++) {
            dE_dy[k] = - t[k] / (y[k] + EPS);
        }
        return 0.0f;
    }
    
    for(k = 1; k <= K; k++) {
        E -= t[k] * logf(y[k] + EPS);
    }
    
    return E;
}
```

# モデルの構造の再定義
活性化関数と損失関数の引数が変わったため， `MODEL_PARAMETER` 型を再定義します。

```C:struct.h
typedef struct {
    int L; //number of hidden layers
    int K; //number of classes
    int dim_in; //input dimension
    int *N_unit; //number of units
    float *(**activation)(float *a, int D, int derivative, float **derivative_matrix); // activation function
    float (*loss)(float *y, float *t, int K, int derivative, float *dE_dy); //loss function
} MODEL_PARAMETER;
```

活性化関数と損失関数を関数ポインタで持っているため， 

`model_parameter.activation[3](a[3], D[3], 0, NULL)` とすれば，

```math
\boldsymbol{z}^{(3)} = h^{(3)}\left(\boldsymbol{a}^{(3)}\right)
```

というように計算できます（逆伝播については[次回]()）。

# 順伝播の実装
モデルの構造が決定したため，順伝播を実装していきます。各層で全結合（full connection）と活性化関数を通す操作を行っていきます。
コードは[こちら](https://github.com/tky1117/Neural_Network_in_C/tree/master/multi_layer_perceptron/code)

## 0. グローバル変数
次の変数をグローバル変数で用意します。

|変数|意味|
|:---:|:---:|
| `x` |入力データの各成分の値|
| `D` |各層のユニット数|
| `w` |各全結合における重み|
| `a` |各層の入力|
| `z` |各層の出力|
| `y` |出力層の出力|
| `t` |正解のone-hot表現|

```C:parameters.h
float *x = NULL;
int *D = NULL;
float ***w = NULL;
float **z = NULL;
float **a = NULL;
float *y = NULL;
float *t = NULL;
```

## 1. 入力層・隠れ層

|変数|意味|
|:---:|:---:|
|$i$|入力層のユニットのインデックス（ $i = 1, 2, \cdots , D_{0}$ ）|
|$j$|隠れ層（第 $1$ 層）のユニットのインデックス（ $j = 1, 2, \cdots , D_{1}$ ）|
|$D_{0}$|入力データの次元|
|$x_{i}$|入力データの各成分の値|
|$D_{1}$|隠れ層（第 $1$ 層）のユニットの個数|
|$a_{j}^{(1)}$|隠れ層（第 $1$ 層）のユニット $j$ が受け取る値|
|$h^{(1)}(\cdot)$|隠れ層（第 $1$ 層）の活性化関数|
|$z_{j}^{(1)}$|隠れ層（第 $1$ 層）のユニット $j$ に活性化関数を通した出力値|

![input_layer.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/297735/70ee0596-19d8-90e4-c103-045e593c781c.png)

|変数|意味|
|:---:|:---:|
|$l$|隠れ層のインデックス（ $l = 1, 2, \cdots , L$ ）|
|$i$|隠れ層（第 $l$ 層）のユニットのインデックス（ $i = 1, 2, \cdots , D_{l}$ ）|
|$a_{i}^{(l)}$|隠れ層（第 $l$ 層）の入力|
|$h^{(l)}(\cdot)$|隠れ層（第 $l$ 層）の活性化関数|
|$z_{i}^{(l)}$|隠れ層（第 $l$ 層）で活性化関数を通した値|
|$j$|隠れ層（第 $l + 1$ 層）のユニットのインデックス（ $j = 1, 2, \cdots , D_{l+1}$ ）|
|$a_{j}^{(l+1)}$|隠れ層（第 $l + 1$ 層）の入力|

![hidden_layer.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/297735/4e3e3fbe-d589-7f86-97a2-d234a26a0ec5.png)

```C:forwardsample.c
void forward(MODEL_PARAMETER model_parameter) {
    int l, i, j, k;
    int L = model_parameter.L;
    int K = model_parameter.K;
    float tmp_a;
    
    // layer: input layer
    z[0] = x;
    
    // layer: hidden layer
    for(l = 1; l <= L; l++) {
        int l_prev = l - 1;
        int d_prev_layer = D[l_prev];
        int d_max = D[l];
        
        for(j = 1; j <= d_max; j++) {
            tmp_a = w[l_prev][0][j]; //bias
            
            for(i = 1; i <= d_prev_layer; i++) {
                tmp_a += w[l_prev][i][j] * z[l_prev][i];
            }
            
            a[l][j] = tmp_a;
        }
        
        z[l] = model_parameter.activation[l](a[l], D[l], 0, NULL);
        // z[l][0] = 1.0f;
    }
    ...
```

## 2. 出力層

|変数|意味|
|:---:|:---:|
|$i$|隠れ層（第 $L$ 層）のユニットのインデックス（ $i = 1, 2, \cdots , D_{L}$ ）|
|$z_{i}^{(L)}$|隠れ層（第 $L$ 層）の出力|
|$k$|出力層のユニットのインデックス（ $k = 1, 2, \cdots , K$ ）|
|$a_{k}^{(L+1)}$|出力層の入力|
|$y_{k}$|出力層の出力|

![output_layer.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/297735/f57d8241-5d0d-11cf-1b55-d7a421f68ad6.png)

```C:forwardsample.c
void forward(MODEL_PARAMETER model_parameter) {
    ...
    // layer: output layer
    int l_prev = L;
    int d_prev_layer = D[l_prev];
    
    for(k = 1; k <= K; k++) {
        tmp_a = w[L][0][k];
        
        for(i = 0; i <= d_prev_layer; i++) {
            tmp_a += w[L][i][k] * z[L][i];
        }
        
        a[L + 1][k] = tmp_a;
    }

    y = model_parameter.activation[L + 1](a[L + 1], K, 0, NULL);
}
```

`forward` 関数を実行することで，推測した確率 $\boldsymbol{y}$ が得られます。 `forward` 関数の中で，損失 $E$ を計算しても良いのですが，あえて `forward` 関数の外で損失の計算を行います[^2]。

```C:forwardsample.c
int main(int argc, char **argv) {
    int N;
    char *file_path = NULL;
    DATA *train_data = NULL;

    ...

    file_path = argv[1];
    train_data = loda_data(file_path, &N);

    MODEL_PARAMETER model_parameter = build_model();

    forward(model_parameter);
    E = model_parameter.loss(y, t, model_parameter.K, 0, NULL);
    ...
}
```

訓練データの1個目のデータを読み込んで損失を計算すると次のようになりました。

```shell
$ gcc forward_sample.c
$ ./a.out ../data/train_data.txt
> E: 1.102192
```

# 続き
[C言語でニューラルネットワークの実装（１）〜多層パーセプトロンの構造と活性化関数〜](https://qiita.com/tky1117/items/ffc367e52c4a04cffa1d)
[C言語でニューラルネットワークの実装（２）〜順伝播と損失関数〜](https://qiita.com/tky1117/items/d2847ac2c84858c7307b)
[C言語でニューラルネットワークの実装（３）〜誤差逆伝播法〜](https://qiita.com/tky1117/items/c2c69a5f69d9cf8ca751)
[C言語でニューラルネットワークの実装（４）〜データの準備〜](https://qiita.com/tky1117/items/7cdebb370e0da7ba6f71)
C言語でニューラルネットワークの実装（５）〜モデルの構造と順伝播の実装〜　←現在の記事
[C言語でニューラルネットワークの実装（６）〜逆伝播の実装〜](https://qiita.com/tky1117/items/83403794280e74357660)
[C言語でニューラルネットワークの実装（７）〜オンライン学習と重みの初期値〜](https://qiita.com/tky1117/items/bbf250bcbd4ddae2d32e)
[C言語でニューラルネットワークの実装（８）〜ミニバッチ学習〜](https://qiita.com/tky1117/items/867f311b91fa927cacd2)

### 更新履歴
- 2020/1/28 ゼロ除算を防ぐようにコードを修正しました。

[^1]: $ \dfrac{e^{a_{k} - C}}{{\displaystyle \sum_{k'=1}^{K}}e^{a_{k'} - C}} = \dfrac{e^{-C}e^{a_{k}}}{e^{-C}{\displaystyle \sum_{k'=1}^{K}}e^{a_{k'}}} = \dfrac{e^{a_{k}}}{{\displaystyle \sum_{k'=1}^{K}}e^{a_{k'}}}$
[^2]: $E$ の値は標準出力などに出力したいが，関数内で計算すると暗黙的で，計算されているのかどうかがわかりにくくなってしまうため。