C言語版『ゼロから作るDeep Learning』①

こんにちは。@Rn86222と申します。この記事は、[ISer Advent Calendar 2022](https://adventar.org/calendars/8285) のために書きました。他の記事も面白いのでぜひ見てみてください。

前回の[C言語版『ゼロから作るDeep Learning』⓪](https://qiita.com/Rn86222/items/58017e2e3bf8d4d20794)の続きです。前回はとりあえず`Function`構造体と`Variable`構造体を用いて自動微分を実装し始めたところでした。各計算における逆伝播は手動で行っていましたが、今回はそこも自動化できるようにしていきます。ステージ1ステップ7から始めます。

https://qiita.com/Rn86222/items/58017e2e3bf8d4d20794

なお**掲示するPythonのコード（*.py）はすべて次のGitHubレポジトリで公開されているものの引用である** ことをここに明示します（前回は基本的にPythonのコードとCのコードを並列する形でしたが、[本](https://www.oreilly.co.jp/books/9784873119069/)の内容説明が主題ではないので今回からはコードの引用は控えめにしておきます）。

https://github.com/oreilly-japan/deep-learning-from-scratch-3
# ステージ1 微分を自動で求める
## ステップ7 バックプロパゲーションの自動化
ステップ6では逆伝播を次のようにしていました。
```c:step06/step06.c
  y.grad = 1.0;
  b.grad = Function_backward((Function*)&C, y.grad);
  a.grad = Function_backward((Function*)&B, b.grad);
  x.grad = Function_backward((Function*)&A, a.grad);
```
このくらいの単純な関数(`Square(Exp(Square(x)))`)であれば上のようにすればよいのですが、もっと複雑な関数を考えるとその計算をいちいち逆からたどるのは大変ですし、そもそも`forward`で行った順伝播の時点でどのように逆伝播すべきかは決まっているのですから、この作業を自動化したくなります。つまり、最後の出力変数`y`に対して`y.backward()`のようにしさえすれば自動で`x`の勾配が求まるようにしたいのです。
ここで以下の計算の流れを可視化したものを見てみましょう。
![graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2726341/d13eee3b-4ea2-2c31-d3d2-9b4690664386.png)
これは計算グラフといって、これ以降も頻繁に出てきます。グラフ中の矢印は順伝播を表していますが、逆伝播はこの矢印を逆にたどっていくことであるので、それを自動化するとなると以下の2つを順伝播の際に把握させることが必要になってきます。
* それぞれの関数の入力はどの変数か
* ぞれぞれの変数はどの関数から作られたか

1つ目の方は既に前回実装しているので、2つ目を実装しましょう。`Variable`構造体のメンバ変数に`Function* p_creator`を追加し、どの関数から作られたのかをセットする関数を定義します。また、後のためにその変数を作った関数が存在するかどうか（存在しない場合とは、ユーザーが手動で`Variable`構造体を定義した場合などであり、上の例で言えば`x`が該当します）を`creator_exists`で判断できるようにしてあります。
```c:step07/variable.h
typedef struct function Function;

typedef struct variable {
  float data;
  float grad;
  int creator_exists;
  Function* p_creator;
} Variable;
// (省略)
void Variable_set_creator(Variable* p_self, Function* func);
```
```c:step07/variable.c
void Variable_set_creator(Variable* p_self, Function* func) {
  p_self->p_creator = func;
  p_self->creator_exists = 1;
}
```
ところでなぜ`step07/variable.h`で`typedef struct function Function;`としているのか疑問に思う方がいるかもしれません。「普通に`#include "function.h"`とすればよいのではないか」と。これには循環参照の問題が絡んでいます。今、`Function`構造体は`Variable input`というメンバ変数を持っていますが、`variable.h`と`function.h`でお互いに`#include`してしまうと、お互いを参照し続けてしまうため上手く行きません。そこで、`typedef struct function Function;`と前方宣言しておけば`Variable`構造体のメンバ変数に`Function*`を用いることが出来ます。もちろんこれだけではこのメンバ変数`Function*`がどんな変数へのポインタなのかわからないので、その実体は後から（コンパイル時に）判明する形です。またなぜポインタにしているのかというと、そうしないと`Variable`構造体のメモリ上でのサイズがわからない（したがってコンパイルできない）ため、というのが1つの理由です[^1]。
したがって、`Function`構造体のメンバ変数の`Variable input`についても`Variable`の前方宣言の上ポインタに直す必要があります。また、ステップ8のときのことを考えて出力の方も覚えさせておくようにします。
```c:step07/function.h
typedef struct variable Variable;
// (省略)
typedef struct function {
  const struct functionmethods *p_methods;
  Variable **p_io;
} Function;
```
```c:step07/function.c
Variable Function_call(Function* p_self, Variable* p_input) {
  p_self->p_io[0] = p_input;
  float y = p_self->p_methods->forward(p_self, p_input->data);
  Variable_init(p_self->p_io[1], y);
  Variable_set_creator(p_self->p_io[1], p_self);
  return *(p_self->p_io[1]);
}
```
`Variable** p_io`というのは、`p_io[0]`で入力変数へのポインタ、`p_io[1]`で出力変数へのポインタにアクセスできるものです。[^2]この場合メモリを動的に確保しておかねばなりませんので、面倒ですが例えば`Square_init`は以下のように変更します。
```c:step07/square.c
void Square_init(Square* p_self) {
  ((Function*)p_self)->p_methods = &SQUARE_METHODS;
  ((Function*)p_self)->p_io = (Variable**)malloc(2 * sizeof(Variable*));
  ((Function*)p_self)->p_io[0] = (Variable*)malloc(sizeof(Variable));
  ((Function*)p_self)->p_io[1] = (Variable*)malloc(sizeof(Variable));
}
```
では以上の準備の下、逆伝播の流れを自動化しましょう！ 出力変数に対する`backward`メソッド（のようなもの）を実装したいので、`Variable_backward`として作ります。
```c:step07/variable.h
void Variable_backward(Variable* p_self);
```
```c:step07/variable.c
void Variable_backward(Variable* p_self) {
  if (p_self->creator_exists) {
    Function* f = p_self->p_creator;
    Variable* x = p_self->p_creator->p_io[0];
    x->grad = f->p_methods->backward(f, p_self->grad);
    Variable_backward(x);
  }
}
```
これは、流れとしては
1. 変数のクリエイターへのポインタを取得
2. 変数のクリエイターへの入力へのポインタを取得
3. 変数のクリエイターへの入力の勾配を`Function`構造体の`backward`から求める
4. 変数のクリエイターへの入力へのポインタに対し`Variable_backward`を適用

という再帰呼び出しになっています。こうすることで、冒頭の`if`文の条件式が偽になる、すなわち関数から作られたのではない変数（入力変数）に至るまで、逆順に変数をたどっていく逆伝播が行われます。
では実際に使ってみましょう。
```c:step07/step07.c
#include <stdio.h>
#include <math.h>
#include "variable.h"
#include "square.h"
#include "function.h"
#include "exp.h"

int main() {
  Variable x;
  Variable_init(&x, 0.5);

  Square A;
  Square_init(&A);
  Exp B;
  Exp_init(&B);
  Square C;
  Square_init(&C);

  Variable a = Function_call((Function*)&A, &x);
  Variable b = Function_call((Function*)&B, &a);
  Variable y = Function_call((Function*)&C, &b);

  y.grad = 1.0;

  Variable_backward(&y);

  printf("%.10f\n", x.grad); // 3.29...

  return 0;
}
```
かなりすっきり書けました！ 自動微分らしくなってきましたね。Pythonの場合と見比べてみましょう。
```py:step07/step07.py
A = Square()
B = Exp()
C = Square()

x = Variable(np.array(0.5))
a = A(x)
b = B(a)
y = C(b)

# backward
y.grad = np.array(1.0)
y.backward()
print(x.grad)
```
やはりPythonの方が大分簡潔になっていますが、Cでも頑張ればいい感じに書けるんだなあと個人的には思います。ステップ7は以上です。
## ステップ8 再帰からループへ
先ほどは逆伝播を再帰関数によって自動化していました。ここではそれを`while`ループを使って書き直します。これは、ステップ15以降でもう少し複雑な計算グラフを扱う場合にも逆伝播が自動で行えるようにするための布石です。Pythonでは以下のように書けます。
```py:step08/step08.py
class Variable:
    #(省略)
    def backward(self):
        funcs = [self.creator]
        while funcs:
            f = funcs.pop()  # 1. Get a function
            x, y = f.input, f.output  # 2. Get the function's input/output
            x.grad = f.backward(y.grad)  # 3. Call the function's backward

            if x.creator is not None:
                funcs.append(x.creator)
```
当然のように`pop`、`append`メソッドが使われていますが、Cにはそんなものはデフォルトでは存在しないので自分で作ります。`append`は`push`のことなのでスタックを実装すればよさそうです。`FunctionStack`構造体として以下のように実装しました。
```step08/functionstack.h
#include "function.h"

#ifndef _FUNCTIONSTACK_H_
#define _FUNCTIONSTACK_H_

typedef struct functionstack {
  Function* stack;
  int last;
  int size;
} FunctionStack;

void FunctionStack_init(FunctionStack* p_self, const int size);
void FunctionStack_push(FunctionStack* p_self, Function func);
Function FunctionStack_pop(FunctionStack* p_self);

#endif // _FUNCTIONSTACK_H_
```
```step08/functionstack.c
#include "functionstack.h"
#include "function.h"
#include <stdio.h>
#include <stdlib.h>

void FunctionStack_init(FunctionStack* p_self, const int size) {
  p_self->last = -1;
  p_self->size = size;
  p_self->stack = (Function* )malloc(size * sizeof(Function));
}

void FunctionStack_push(FunctionStack* p_self, Function func) {
  if (p_self->last >= p_self->size - 1) {
    printf("FunctionStack_push error: stack overflow.\n");
    return;
  }
  p_self->stack[++p_self->last] = func;
}

Function FunctionStack_pop(FunctionStack* p_self) {
  if (p_self->last == -1) {
    printf("FunctionStack_pop error: no function in stack.\n");
    Function tmp;
    return tmp;
  }
  return p_self->stack[p_self->last--];
}
```
まあただスタックを作っているだけなのでここはいいでしょう。では`FunctionStack`を用いてループで逆伝播の自動化を書いてみます。
```c:step08/variable.c
void Variable_backward(Variable* p_self) {
  if (p_self->creator_exists == 0) {
    return;
  }
  const int stack_size = 100;
  FunctionStack fs;
  FunctionStack_init(&fs, stack_size);
  FunctionStack_push(&fs, *(p_self->p_creator));

  Function f;
  while (fs.last != -1) {
    f = FunctionStack_pop(&fs);
    Variable* x = f.p_io[0];
    Variable* y = f.p_io[1];
    x->grad = f.p_methods->backward(&f, y->grad);
    if (x->creator_exists) {
      FunctionStack_push(&fs, *(x->p_creator));
    }
  }
}
```
前ステップで関数に出力を覚えさせておいたのは、このようにスタックから取り出した関数が入力側に伝えるべき勾配を知るためです。この`Variable_backward`は内部実装が変わっただけなので前ステップと同様に使うことができて、同じ結果が得られることが確かめられます。ステップ8はこれだけです。
## ステップ9 関数をより使いやすく
ステップ9は、まず次のようにすることで関数呼び出しを楽にしています。
```py:step09/step09.py
def square(x):
    return Square()(x)

def exp(x):
    return Exp()(x)

x = Variable(np.array(0.5))
y = square(exp(square(x)))
y.backward()
print(x.grad)
```
Cで同じようなことをするのは少し難しそうですが、やればできます。以下のようにします。
```c:step09/step09.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "variable.h"
#include "square.h"
#include "function.h"
#include "exp.h"

Variable* square(Variable* x){
  Square* sq;
  sq = (Square*)malloc(sizeof(Square));
  Square_init(sq);
  return Function_call((Function*)sq, x);
}

Variable* expo(Variable* x){
  Exp* ex;
  ex = (Exp*)malloc(sizeof(Exp));
  Exp_init(ex);
  return Function_call((Function*)ex, x);
}

int main() {
  Variable x;
  Variable_init(&x, 0.5);

  Variable* y = square(expo(square(&x)));

  Variable_backward(y);

  printf("%.10f\n", x.grad); // 3.29...
  
  return 0;
}
```
すごくいい感じに書けました！[^3]　`exp`ではなく`expo`としているのは関数名の衝突を回避するためです。また、前ステップからの変更として`Variable`と`Variable_backward`に少し手を加えています。
```c:step09/variable.h
typedef struct variable {
  float data;
  float grad;
  int creator_exists;
  int grad_exists;
  Function* p_creator;
} Variable;
```
```c:step09/variable.c
void Variable_backward(Variable* p_self) {
  if (p_self->creator_exists == 0) {
    return;
  }
  if (p_self->grad_exists == 0) {
    p_self->grad = 1.0;
    p_self->grad_exists = 1;
  }
  // (省略)
}
```
このようにしておけばいちいち最後の出力変数に対し`y.grad = 1.0`などとしていたのを省略できます（上の`step09.c`は実際に省略しています）。
この他DeZeroでは`Variable.data`としてndarrayのみを扱うように色々と工夫をしています。...そういえば未だにndarrayを実装していないのでした。ですがそれはステップ16あたりまで待ってください。とりあえずしばらくは代わりに`float`を扱うということで問題ないはずです。というわけでステップ9は以上です。
## ステップ10 テストを行う
ステップ10にはDeZeroのテストを行うやり方がいくつか示されています。例えば逆伝播が上手くいっているか判断するために数値微分（ステップ4参照）を行ってその結果と比較する、などです。ただこれらはPythonのunittestモジュールを使って実装されているのでCでは実装しません（期待していた方々がいたら申し訳ありません）。逆伝播が上手くいっているかどうかは（DeZeroを信頼して）DeZeroの結果と一致するかどうかで判断することにします。

# ステップ7～10総括
というわけで（やや消化不良感は否めませんが）ステージ1が終わりました！ キリがいいので今回の記事はここまでにします。最初始めた頃は本当にできるのか不安でしたが、だんだん軌道に乗ってきたように感じています。ステージ2では`Function`構造体を可変長引数に対応させたり、次のような少し複雑な計算グラフに対しても正しく逆伝播が行えるようにしていったりします。
![graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2726341/a9152d7d-f709-6f16-1f3b-3ac38b975d2d.png)
実装自体は既に済んでいるので、~~[アドカレ](https://adventar.org/calendars/8285)を埋めるためにも~~早く続きを書きたいと思います。
ここまで読んでくださりありがとうございました。

[^1]: ただ、今回の例で言えば仮にサイズがわかったとしてもポインタにすべきです。なぜなら、変数と関数の間につながりを作ることが大事だからです。ポインタでなく単に`Function creator`に`Function`構造体を代入するとなると、それは単にコピーをしているだけであり、本当のクリエイターとのつながりを作ることができません。それで実用上しばらく問題ないとしても、ステップ26で計算グラフを自動で可視化する際に問題になってきます。
[^2]: 普通に`Variable* p_input`と`Variable* p_output`をメンバ変数にしておいた方がいいのでしょう。私もそう思って最初はそのようにしたのですが、どうやら実体のない構造体(`Variable`)へのポインタを2つメンバ変数に持たせるのはメモリ確保の際上手くいかないようです（私が間違っているだけかもしれないので詳しい方いたら教えてください）。
[^3]: まあこれだとmallocで確保したメモリ領域を解放するのが難しくなってきますね。そのあたりもステップ16あたりで対応していければと思っています。