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

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

# 『ゼロから作るDeep Learning』とは
[『ゼロから作るDeep Learning』](https://www.oreilly.co.jp/books/9784873117584/)といえば、[斎藤 康毅さん](https://twitter.com/saitohkoki)が書かれた言わずと知れたディープラーニングの入門書シリーズで、よく『ゼロつく』と呼ばれています。昨今流行りに流行っている深層学習を、PyTorchやTensorflowなどの外部ライブラリを用いることなくなるべく自分の手でやってみようというもので、現在シリーズで4巻まで刊行されています。私はこのシリーズがとても好きで、機械学習・深層学習の勉強をする際大変お世話になっています。

# 「ゼロ」から？
ところで、『ゼロつく』のサンプルコードはいずれもPython 3系（以下単にPython）で書かれています。実際深層学習をしようとなればPythonを使うのがかなりメジャーな感じなので、当然といえば当然です。また『ゼロつく』において用いられるライブラリはNumPy（数値計算用ライブラリ）やMatplotlib（グラフ等の可視化用ライブラリ）あたりに限られているので、かなりゼロに近いところからの深層学習の実装が体験できます。しかし、以下のような感想を抱く人がいるかもしれません。[^1]
* NumPyずるい（数値計算などが大分楽になっている）
* 動的型付けずるい（関数定義等の際に苦労しない）
* Pythonずるい（もともと機能が豊富）
* Pythonが嫌い（私情）

そのような人のために（？）、より「ゼロ」から深層学習を実装してみようと思いました。
# C言語で『ゼロつく』
というわけで、最近私が大学の勉強で触れることの多いC言語（以下単にC）でやってみようということになりました。Cを少しでもやったことのある人はすでにお察しのように、Cで深層学習をやるというのは縛りプレイも甚だしいです。というか常人には無理なんじゃないかという感じもします。まあ無理だったとしてもある程度のところまでいければいいやという精神で、始めていこうと思います。
## 何をするのか
まず『ゼロつく』をやるといっても何をやるのかという話ですが、[『ゼロから作るDeep Learning ❸ ―フレームワーク編』](https://www.oreilly.co.jp/books/9784873119069/)をやります（以下単に『ゼロつく』といったらこれのことを指す）。これはシリーズ3巻にあたり、DeZeroという深層学習のフレームワークをゼロから作ろうというものです。ステップが1～60までに細かく分かれており、初心者にもわかりやすく構成されています。それでもCでやるには充分に大変なのですが...。

## ルール
一応実装におけるルールを定めておきます。が、いかんせんやってみないことには先が見えないので、後から変更する可能性もあります。
* C言語の規格はC11で、対応する標準Cライブラリはすべて使用可能
* コンパイラはgcc (MinGW.org GCC Build-2) 9.2.0
* gnuplotとGraphvizは使用可能
* 上記とエディタ等を除くすべての外部ライブラリおよびソフトウェアは使用禁止

まあ上2つは変更されることはないでしょう。問題はグラフの可視化です。『ゼロつく』ではMatplotlibライブラリ、Graphvizが使われています。まあこれらが使われるのはステップ25以降なので、一旦上のようにしておきます。

前置きが長くなりましたが、いよいよ実装を始めていきます。Cはまだまだ初心者なので何か間違い等あれば何でも教えてください。また **掲示するPythonのコード（*.py）はすべて次のGitHubレポジトリで公開されているものの引用である** ことをここに明示します。

https://github.com/oreilly-japan/deep-learning-from-scratch-3

# 第1ステージ 微分を自動で求める（ステップ1～10）
## ステップ1 箱としての変数
最初はDeZeroの変数を`Variable`クラスとして実装します。Pythonではこうです。
```python:steps/step01.py
class Variable:
    def __init__(self, data):
        self.data = data
```
...いきなりCには存在しないクラスという概念が出てきました。まあCにも構造体があるので何とか耐えてほしいところです。とりあえず`variable.h`と`variable.c`を次のようにしました。メソッドはとりあえず構造体の外部で構造体の名前を最初につけて定義することにします。
```C:step01/variable.h
#ifndef _VARIABLE_H_
#define _VARIABLE_H_

typedef struct variable {
  float data;
} Variable;

void Variable_init(Variable* p_self, const float data);

#endif // _VARIABLE_H_
```
```step01/variable.c
#include "variable.h"

void Variable_init(Variable* p_self, const float data) {
  p_self->data = data;
}
```
ところでDeZeroでは`Variable`クラスのインスタンス変数`data`に入れるものとしてはNumPyの多次元配列(ndarray)のみを扱うことにしています。以下のような感じです。
```python:steps/step01.py
data = np.array(1.0)
x = Variable(data)
print(x.data)
```
「え、じゃあいきなりNumPyみたいなものを実装しなきゃいけないのでは...」となりましたが、それだとちょっといつまでも先に進めないので、一旦`Variable`構造体のメンバ変数`data`の型は`float`としておきました。使い方は以下の通りです。
```step01/step01.c
#include <stdio.h>
#include "variable.h"

int main() {
  Variable x;
  Variable_init(&x, 1.0);
  printf("%.2f\n", x.data); // 1.00
  return 0;
}
```
ndarrayがないという問題はさておいて、とりあえずステップ1は完了です。

## ステップ2 変数を生み出す関数
ステップ2では、ステップ1で作った`Variable`クラスのインスタンスを入力として受け取り、`Variable`インスタンスを出力する関数のクラスである`Function`クラスを実装します。
```python:steps/step02.py
class Function:
    def __call__(self, input):
        x = input.data
        y = self.forward(x)
        output = Variable(y)
        return output

    def forward(self, in_data):
        raise NotImplementedError()
```
`Function`クラスは上のようにしておいて、実際の計算は`Function`クラスを継承した具象クラスで`forward`メソッドをオーバーライド（上書き）する形です。例えば`Square`クラスは以下のようにします。
```python:steps/step02.py
class Square(Function):
    def forward(self, x):
        return x ** 2
```
Cの構造体にはメソッドすらないのに、継承とかオーバーライドとか言われるともう投げ出したくなりますが、まだ諦めるには早いです。Cでオブジェクト指向を表現するやり方を先人が開拓してくださっているので、それを参考に進めていくことにします。参考にしたものは以下の2つの記事です。

https://qiita.com/qiita_kuru/items/8f3441bce9b2f53d1d62#c%E8%A8%80%E8%AA%9E%E3%81%AE%E5%A0%B4%E5%90%88-1

https://qiita.com/qiita_kuru/items/19d501b364c7e32b3576

ではまず`Function`構造体については以下のようにします（新しく構造体を作るたびにそれ用の`.h`ファイルと`.c`ファイルを作ることにしています）。
```c:step02/function.h
#include "variable.h"

#ifndef _FUNCTION_H_
#define _FUNCTION_H_

struct functionMethods;

typedef struct function {
  const struct functionmethods *p_methods;
} Function;

typedef struct functionmethods {
  float (*forward)(Function* const, const float);
} FunctionMethods;

Variable Function_call(Function* p_self, const Variable input);

#endif // _FUNCTION_H_
```
```step02/function.h
#include "function.h"
#include "variable.h"

Variable Function_call(Function* p_self, const Variable input) {
  float y = p_self->p_methods->forward(p_self, input.data);
  Variable output;
  Variable_init(&output, y);
  return output;
}
```
まず`Function`構造体に`FunctionMethods`という関数ポインタ[^2]の配列`p_methods`をメンバ変数として持たせておくことで、`p_self->p_methods->forward`のようにして、構造体がまるでメソッドを持つかのように見せかけることが出来ます。さらに、こうしておくことでクラスの継承のみならずメソッドのオーバーライドまでもが再現（？）できます。例えば`Square`構造体を以下のようにすればよいです。
```c:step02/square.h
#include "function.h"

typedef struct square {
  Function function;
} Square;

void Square_init(Square* p_self);
```
```c:step02/square.c
#include "square.h"
#include "function.h"
#include "math.h"

static float square_forward(Function* const p_self, const float x) {
  return pow(x, 2);
}

static const FunctionMethods SQUARE_METHODS = {
  square_forward
};

void Square_init(Square* p_self) {
  ((Function*)p_self)->p_methods = &SQUARE_METHODS;
}
```
重要なポイントは2つあります。まず`Square`構造体のメンバ変数として最初に`Function`構造体を持たせておくことです。このようにすれば、`Square`構造体へのポインタに対し`(Function*)`とキャストすることができ、それができれば`Function`構造体用に作った関数（メソッド）を安全に用いることができます（継承）。2つ目のポイントは、`Square_init`で関数ポインタ`p_methods`を書き換えていることです。これでオーバーライドのようなことができます。
DeZeroでは`Function`クラスを以下のように使います。
```python:steps/step02.py
x = Variable(np.array(10))
f = Square()
y = f(x)
print(type(y))
print(y.data)
```
Cの場合はこうです。
```c:step02/step02.c
#include <stdio.h>
#include "variable.h"
#include "square.h"
#include "function.h"

int main() {
  Variable x;
  Variable_init(&x, 10);
  Square f;
  Square_init(&f);
  Variable y = Function_call((Function*)&f, x);
  printf("%.2f\n", y.data); // 100.00
  return 0;
}
```
**色々煩わしいですね...。** 構造体の変数を作る際に「宣言→init関数を呼び出す」という二度手間を取らなきゃいけないだけでなく、いちいち`(Function*)`とキャストしなければならないなど、可読性はかなり落ちています。まあ縛りプレイなので仕方ないです。とりあえずステップ2は完了です。
## ステップ3 関数の連結
ステップ3では`Square`クラスに続いて`Exp`クラスを作ります。
```python:steps/step03.py
class Exp(Function):
    def forward(self, x):
        return np.exp(x)
```
Cでは次のようにします。
```C:step03/exp.h
#include "function.h"

typedef struct exp {
  Function function;
} Exp;

void Exp_init(Exp* p_self);
```
```C:step03/exp.c
#include "exp.h"
#include "function.h"
#include <math.h>

static float exp_forward(Function* const p_self, const float x) {
  return exp(x);
}

static const FunctionMethods EXP_METHODS = {
  exp_forward
};

void Exp_init(Exp* p_self) {
  ((Function*)p_self)->p_methods = &EXP_METHODS;
}
```
`Square`構造体とほとんど同じで`forward`の部分だけ変更していることが分かると思います。さて、今`Function`クラス(構造体)の入力と出力はいずれも`Variable`であるので、`Function`を「連結」することが可能です。DeZeroならば次のようになります。
```steps/step03.py
A = Square()
B = Exp()
C = Square()

x = Variable(np.array(0.5))
a = A(x)
b = B(a)
y = C(b)
print(y.data)
```
Cの場合は...かなり面倒になってきます。
```c:step03/step03.c
#include <stdio.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);

  printf("%.10f\n", y.data); // 1.648...
  return 0;
}
```
たかが $y=(e^{x^2})^2$ を計算するだけでこんなに書かなきゃいけないとは...。先が思いやられます。
ステップ3はこれだけです。
## ステップ4 数値微分
さて、いよいよ微分のパートに入っていきます。まずは簡単な数値微分です。数値微分とは、微分の定義式 $f'(x)=\lim _{h\rightarrow0}\frac{f(x+h)-f(x)}{h}$ を近似して、微小な値$\varepsilon$を用いることで $f'(x)\approx\frac{f(x+\varepsilon)-f(x)}{\varepsilon}$ とするやり方です。この式をそのまま用いるよりも中心差分近似$f'(x)\approx\frac{f(x+\varepsilon)-f(x-\varepsilon)}{2\varepsilon}$ を用いた方が精度がよいことが知られているので、こちらを用います。
では数値微分を実装してみます。Pythonでは次のようになります。
```python:steps/step04.py
def numerical_diff(f, x, eps=1e-4):
    x0 = Variable(x.data - eps)
    x1 = Variable(x.data + eps)
    y0 = f(x0)
    y1 = f(x1)
    return (y1.data - y0.data) / (2 * eps)
```
引数fはFunctionクラスのインスタンスです。次のように用いることが出来ます。
```python:steps/step04.py
def f(x):
    A = Square()
    B = Exp()
    C = Square()
    return C(B(A(x)))


x = Variable(np.array(0.5))
dy = numerical_diff(f, x)
print(dy)
```
では同じことをCでやってみましょう。まず`numerical_diff`は次のように書けます。
```c:step04/step04.c
float numerical_diff(Function* f, Variable x) {
  float eps = pow(10, -4);
  Variable x0;
  Variable_init(&x0, x.data - eps);
  Variable x1;
  Variable_init(&x1, x.data + eps);
  Variable y0 = Function_call(f, x0);
  Variable y1 = Function_call(f, x1);
  return (y1.data - y0.data) / (2 * eps);
}
```
本質的な話ではないのですが、ここでもCとPythonの差があって、Pythonでは引数の`eps`にデフォルト値を与えておくことで必ずしも引数として与える必要はなくなっているものの、Cではそのようなことはできません。可変長引数を用い、引数をいくつ与えたかを判断できるようにすれば何とか再現できるのですが、関数の引数および内容が冗長になることが多いです。今回はとりあえず引数の`eps`はなくして関数内で値を決めておくことにしました。使い方は次のようになります。
```c:step04/step04.c
#include <stdio.h>
#include <math.h>
#include "variable.h"
#include "square.h"
#include "function.h"
#include "exp.h"

float numerical_diff(Function* f, Variable x) {
  // (省略)
}

typedef struct abc {
  Function function;
} ABC;

static float abc_forward(Function* const p_self, const float x) {
  Variable input;
  Variable_init(&input, x);
  Square A;
  Square_init(&A);
  Exp B;
  Exp_init(&B);
  Square C;
  Square_init(&C);

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

  return y.data;
}

static const FunctionMethods ABC_METHODS = {
  abc_forward
};

void ABC_init(ABC* p_self) {
  ((Function*)p_self)->p_methods = &ABC_METHODS;
}

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

  ABC f;
  ABC_init(&f);

  float dy = numerical_diff((Function*)&f, x);

  printf("%.10f\n", dy); // 3.29...
  return 0;
}
```
`numerical_diff`を使うためだけに新しく`ABC`構造体を`Function`構造体の「継承」として作らなければならないという面倒臭さ...。まあともかくPythonの場合と実行結果が大体一致したので実装に問題はなさそうです。とりあえず数値微分を用いていくつかの関数の組み合わせに対しその微分を自動で求められるようになりました。ステップ4は以上なのですが、数値微分には誤差が大きかったり計算コストが大きかったりするため、実際に機械学習で用いられることは滅多にないと思います。そこでいよいよ合成関数の微分、特に逆伝播法を用いた自動微分に移っていきます。
## ステップ5　バックプロパゲーションの理論
ここは逆伝播の理論についてなので省略します。『ゼロつく』ではわかりやすく解説されているのでぜひ本を買って見てみてください。
## ステップ6　手作業によるバックプロパゲーション
では逆伝播の実装です。DeZeroではまず`Variable`クラスのインスタンス変数に`grad`（勾配）を追加します。
```python:steps/step06.py
class Variable:
    def __init__(self, data):
        self.data = data
        self.grad = None
```
そして`Function`クラスについては、まず逆伝播において必須の「入力値を覚えておく機能」としてインスタンス変数に`input`を追加し、さらに逆伝播を行うための`backward`メソッドを追加します。
```python:steps/step06.py
class Function:
    def __call__(self, input):
        x = input.data
        y = self.forward(x)
        output = Variable(y)
        self.input = input
        return output

    def forward(self, x):
        raise NotImplementedError()

    def backward(self, gy):
        raise NotImplementedError()
```
そして`backward`メソッドの中身を`Square`クラスや`Exp`クラスについて書きます。
```python:steps/step06.py
class Square(Function):
    def forward(self, x):
        y = x ** 2
        return y

    def backward(self, gy):
        x = self.input.data
        gx = 2 * x * gy
        return gx


class Exp(Function):
    def forward(self, x):
        y = np.exp(x)
        return y

    def backward(self, gy):
        x = self.input.data
        gx = np.exp(x) * gy
        return gx
```
各々のクラスの微分については手計算で求めておく必要があります（例えば$f(x)=x^2$ であれば$f'(x)=2x$ ）が、このようにすることでこれらを組み合わせてできる複雑な関数についても自動で微分ができるわけです。次のようにします。
```python:steps/step06.py
A = Square()
B = Exp()
C = Square()

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

y.grad = np.array(1.0)
b.grad = C.backward(y.grad)
a.grad = B.backward(b.grad)
x.grad = A.backward(a.grad)
print(x.grad)
```
出力側から入力側の順番で勾配を求め、伝播させていきます（トップダウン方式といいます）。
    ではCで再現していきましょう。まず`Variable`構造体のメンバ変数として`grad`を追加しますが、DeZeroでは初期値は`None`でした。Cには`None`にあたるものは恐らくないので、とりあえず初期値は`0.0`としておきます。
```c:step06/variable.h
typedef struct variable {
  float data;
  float grad;
} Variable;
```
```c:step06/variable.c
void Variable_init(Variable* p_self, const float data) {
  p_self->data = data;
  p_self->grad = 0.0;
}
```
続いて`Function`構造体ですが、こちらは色々と変更点があります。
```c:step06/function.h
#include "variable.h"

#ifndef _FUNCTION_H_
#define _FUNCTION_H_

struct functionMethods;

typedef struct function {
  const struct functionmethods *p_methods;
  Variable input;
} Function;

typedef struct functionmethods {
  float (*forward)(Function* const, const float);
  float (*backward)(Function* const, const float);
} FunctionMethods;

Variable Function_call(Function* p_self, const Variable input);
float Function_backward(Function* p_self, const float y);

#endif // _FUNCTION_H_
```
```c:step06/function.c
#include "function.h"
#include "variable.h"

Variable Function_call(Function* p_self, const Variable input) {
  p_self->input = input;
  float y = p_self->p_methods->forward(p_self, input.data);
  Variable output;
  Variable_init(&output, y);
  return output;
}

float Function_backward(Function* p_self, const float y) {
  float dy = p_self->p_methods->backward(p_self, y);
  return dy;
}
```
まずメンバ変数に`input`を加えます。さらに関数ポインタの配列である`p_methods`に`backward`を追加します。そして`Function_init`の中で`input`を覚えさせています。最後に、`Function_backward`関数を作ります。これは`Function_call`と同様に`p_self->p_methods->backward`のようにして「`backward`メソッド」を呼び出すものです。
そして`Square`構造体および`Exp`構造体で具体的な逆伝播を実装し、「メソッド」に追加します。
```c:step06/square.c
static float square_backward(Function* const p_self, const float gy) {
  float gx = 2 * p_self->input.data * gy;
  return gx;
}

static const FunctionMethods SQUARE_METHODS = {
  square_forward,
  square_backward
};
```
```c:step06/exp.c
static float exp_backward(Function* const p_self, const float gy) {
  float gx = exp(p_self->input.data) * gy;
  return gx;
}

static const FunctionMethods EXP_METHODS = {
  exp_forward,
  exp_backward
};
```
では実際に逆伝播を行ってみます。
```c:step06/step06.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;
  b.grad = Function_backward((Function*)&C, y.grad);
  a.grad = Function_backward((Function*)&B, b.grad);
  x.grad = Function_backward((Function*)&A, a.grad);

  printf("%.10f\n", x.grad); // 3.29...
  return 0;
}
```
そこそこ簡単に書けましたね（錯覚）。再び結果が一致したので上手くいってそうです。Cだけで逆伝播まで実装できました！ステップ6は以上です。

## ステップ1～6総括
長くなってきたので一旦ここまでとします。縛りプレイはやはり楽しいですね。一応手元ではステップ20あたりまで実装できているので、また時間あるときに続きを書いていきます。よければ「いいね」してくださるととても喜びます。
読んでくださりありがとうございました。

続きはこちら。

https://qiita.com/Rn86222/items/72ec8669937a1ac171fd

[^1]:念のため言いますが私は『ゼロつく』はかなりの名著だと思っています。
[^2]:関数ポインタについては[こちらの記事](https://monozukuri-c.com/langc-funcpointer/)がわかりやすいです。