こんにちは。@Rn86222と申します。この記事は、ISer Advent Calendar 2022 のために書きました（3投稿目です）。他の記事も面白いのでぜひ見てみてください。前回のC言語版『ゼロから作るDeep Learning』①の続きです。前回は自動微分を書きやすくして第1ステージが終わりました。今回は第2ステージから始めます。

なお掲示するPythonのコード（*.py）はすべて次のGitHubレポジトリで公開されているものの引用である ことをここに明示します。

第2ステージ自然なコードで表現する

ステップ11 可変長の引数（順伝播編）

これまで扱ってきたFunctionクラス（構造体）は、入力も出力も1つであることを想定していました。実際SquareもExpも入出力はそれぞれ1つです。しかし当然のことながら、一般の関数には入力か出力あるいはその両方が複数であるものがあります。そこでまずは、Functionクラスへの入力と出力をVariableのリスト（Cでは配列）に変更します。
Pythonでは次のようにします。

step11.py

class Function:
    def __call__(self, inputs):
        xs = [x.data for x in inputs]  # Get data from Variable
        ys = self.forward(xs)
        outputs = [Variable(as_array(y)) for y in ys]  # Wrap data

        for output in outputs:
            output.set_creator(self)
        self.inputs = inputs
        self.outputs = outputs
        return outputs

Pythonならではの内包表記を使うことで綺麗に書けています。また内包表記に限らず、Pythonはリストの要素をすべて順に取り出すということが標準でできるのでだいぶ書きやすいです。しかしCの配列で同じことをするためには配列の要素数を知っておく必要があります。そこでまずはVariable構造体のメンバ変数にinput_numおよびoutput_numを追加します。

step11/function.h

typedef struct function {
  const struct functionmethods *p_methods;
  Variable ***p_io;
  int input_num;
  int output_num;
} Function;

そしてFunction_initを呼び出す際にこれらを引数として受け取るようにします。

step11/function.c

void Function_init(Function* p_self, const int input_num, const int output_num) {
  p_self->p_io = (Variable***)malloc(2 * sizeof(Variable**));
  p_self->p_io[0] = (Variable**)malloc(input_num * sizeof(Variable*));
  p_self->p_io[1] = (Variable**)malloc(output_num * sizeof(Variable*));
  p_self->input_num = input_num;
  p_self->output_num = output_num;
  // 省略
}

これに伴い、Square_initやExp_initに変更を加える必要があります。

step11/square.c

void Square_init(Square* p_self) {
  ((Function*)p_self)->p_methods = &SQUARE_METHODS;
  Function_init((Function*)p_self, 1, 1); // input_num = 1, output_num = 1
}

こんな感じで、Function構造体を継承した構造体のinit内において入出力数を指定することになります。
では準備ができたのでCのFunction_callを変更しましょう。

step11/function.c

Variable** Function_call(Function* p_self, Variable** p_inputs) {
  float* xs;
  xs = (float*)malloc(p_self->input_num * sizeof(float))
  for (int i = 0; i < p_self->input_num) {
    p_self->p_io[0][i] = p_inputs[i];
    xs[i] = p_inputs[i].data;
  }
  float* ys;
  ys = (float*)malloc(p_self->output_num * sizeof(float));
  ys = p_self->p_methods->forward(p_self, xs);
  for (int i = 0; i < p_self->output_num; i++) {
    Variable_init(p_self->p_io[1][i], ys[i]);
    Variable_set_creator(p_self->p_io[1][i], p_self);
  }
  free(xs);
  free(ys);
  return p_self->p_io[1];
}

やはりCで書くなら多少の処理の煩雑さは受け入れなければなりませんね...。ところで、コードを見ると分かるように、Function_callの返り値の型をVariable**（意図としてはVariable*の配列）に変更するとともに、forwardメソッドの第2引数を配列に変更しています。したがってまずFunctionMethodsを以下のように修正します（backwardの第2引数も配列にしておきます）。

step11/function.h

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

そしてSquareなども再び変更する必要があるのですが、それはもう少し可変長引数の実装を改良してからにすることにして、一旦新しくAdd構造体を作ります。これはその名の通り足し算をするもので、入力数が2になります。

step11/add.h

#include "function.h"

#ifndef _ADD_H_
#define _ADD_H_

typedef struct add {
  Function function;
} Add;

void Add_init(Add* p_self);

#endif // _ADD_H_

step11/add.c

#include "function.h"
#include "variable.h"
#include "add.h"
#include <stdlib.h>

static float* add_forward(Function* const p_self, const float* xs) {
  float* ys;
  ys = (float*)malloc(sizeof(float));
  ys[0] = xs[0] + xs[1];
  return ys;
}

static float* add_backward(Function* const p_self, const float* gy) {
 // あとで実装
}

static const FunctionMethods ADD_METHODS = {
  add_forward,
  add_backward
};

void Add_init(Add* p_self) {
  ((Function*)p_self)->p_methods = &ADD_METHODS;
  Function_init((Function*)p_self, 2, 1);
}

一旦forwardのみ作りました。Add_initで入力数として2を指定していることに注意してください。ではこのAdd構造体を使ってみましょう。

step11/step11.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "variable.h"
#include "function.h"
#include "add.h"

Variable** add(Variable** xs){
  Add* ad;
  ad = (Add*)malloc(sizeof(Add));
  Add_init(ad);
  return Function_call((Function*)ad, xs[0], xs[1]);
}

int main() {
  Variable x0;
  Variable x1;
  Variable_init(&x0, 2.0);
  Variable_init(&x1, 3.0);
  Variable xs[2] = {x0, x1};

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(&xs);

  printf("%.2f\n", y[0]->data); // 5.00

  return 0;
}

ただ足し算をするだけでこんなに書かなきゃいけないとは...。とりあえずステップ11はここまでです。

ステップ12 可変長の引数（改善編）

では可変長変数をもう少し使いやすくしましょう。先ほどFunction_callでは入力を配列にしていましたが、入力として与えたい変数をわざわざ配列として与えるのは面倒です。つまりSquareであればsquare(&x)、Addであればadd(&x0, &x1)のようにして呼び出すのが自然でしょうし便利です。ではもっと可変長らしい書き方にしましょう。Pythonでは次のような記法があります。

step12.py

class Function:
    def __call__(self, *inputs):
        xs = [x.data for x in inputs]

このように、引数にアスタリスク*をつけると任意の数の引数を受け取れます。これがPythonにおける本当の可変長引数です。
ではCではどうするのか。「可変長引数なんてなくないか...？」と最初思ったのですが、そういえば一番よく使う関数が可変長引数をばりばり使っていました。そう、printfです。printfは次のようにして任意の数の引数を受け取れます。

printf("%f %f\n", x, y);

scanfも同様ですね。つまり、Cにも引数を可変長にする機能はあるわけです。ではそんな関数をどのように定義するのかというと、stdarg.hに便利なマクロが用意されています。以下のような感じです。

#include <stdio.h>
#include <stdarg.h>

int sample_sum(int input_num, ...) {  // `...`で可変長引数にする
  va_list va_ptr;  // va_list型の変数を宣言
  va_start(va_ptr, input_num);  // `...`の直前の引数の名前を与えて初期化
  int sum = 0;
  for (int i = 0; i < input_num; i++)
    sum += va_arg(va_ptr, int);  // 型を指定して可変長引数を1つずつ受け取る
  va_end(va_ptr);  // 使い終わったらva_end
  return sum;
}

int main() {
  printf("%d\n", sample_sum(3, 1, 2, 3)); // 6
  return 0;
}

とりあえず最低限の使い方はコメントを見ていただければわかるかと思います。ただ、Pythonの可変長引数とは異なる点がいくつかあって、

すべての引数を可変長にはできない（例. sample_sum(...)）
可変長部分でいくつの変数を受け取ったのかを直接知ることは（恐らく）できない

などの制限があります。特に2つ目には注意しなければならず、実際に受け取った可変長部分の引数の数よりも多くva_argを呼び出すと未定義動作になるようです。それを防ぐために、上の例では引数にinput_numを与えています（これによち1つ目の方の問題も解消されます）。
ではこれをFunction_callに使ってみましょう。次のようになります。

step12/function.c

Variable** Function_call(Function* p_self, ...) {
  va_list va_ptr;
  va_start(va_ptr, p_self);
  float* xs;
  xs = (float*)malloc(p_self->input_num * sizeof(float));
  for (int i = 0; i < p_self->input_num; i++) {
    p_self->p_io[0][i] = va_arg(va_ptr, Variable*);
    xs[i] = p_self->p_io[0][i]->data;
  }
  va_end(va_ptr);
  // 省略

今回の場合p_selfが可変長部分の引数の数の情報をメンバ変数にもっているので、va_argを使う際にはそれを利用しています。Pythonよりもやたら書くことが多くなりましたが、これで立派な可変長引数の実装ができました。
ところでPythonでは以下のような改良も加えています。

step12.py

class Function:
    def __call__(self, *inputs):
        # 省略
        return outputs if len(outputs) > 1 else outputs[0]

このようにすることで、出力が1個の場合は長さ1のリストではなくその最初の要素を返すようにでき、使いやすくなります。しかしCでは、これをやろうとすると出力の個数に応じて返り値の型が変わることになってしまうので、ここの工夫は断念します（というかできたとしてもCではそのようなことはやらない方がいいでしょう）。
さらに、SquareやAddのforwardメソッドの実装を楽にするために以下のような改善もしています。

step12.py

class Function:
    def __call__(self, *inputs):
        xs = [x.data for x in inputs]
        ys = self.forward(*xs)
        if not isinstance(ys, tuple):
            ys = (ys,)
        # 省略

xsはリストですが*xsとするとその要素が展開されます。またforwardの出力がタプルでない場合はタプルに変換しています。このようにしておくことで、例えばAddであればforwardメソッドを次のように定義できます。

step12.py

class Add(Function):
    def forward(self, x0, x1):
        y = x0 + x1
        return y

Pythonの機能を・書き方を有効に使っていますね。ちょっとCで再現するのは難しそうです...。まあやはりどうせCで書くなら書きやすさよりも安全性を優先しましょう。¹とりあえず1つ目の改善のみを施したら、今度は次のようにAddを使えます。

step12/step12.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "variable.h"
#include "function.h"
#include "add.h"

Variable** add(Variable* x, Variable* y){
  Add* ad;
  ad = (Add*)malloc(sizeof(Add));
  Add_init(ad);
  return Function_call((Function*)ad, x, y);
}

int main() {
  Variable x0;
  Variable x1;
  Variable_init(&x0, 2.0);
  Variable_init(&x1, 3.0);

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(&x0, &x1);

  printf("%.2f\n", y[0]->data); // 5.00

  return 0;
}

前ステップよりはましになりました。ステップ12は以上です。

ステップ13 可変長の引数（逆伝播編）

これまでで実装してきたのは主に可変長引数の順伝播でしたが、続いて可変長引数に対応した逆伝播を実装していきます。まずAddですが、$f(x,y)=x+y$ に対し$\frac{\partial{f}}{\partial{x}}=\frac{\partial{f}}{\partial{y}}=1$であるので、逆伝播は次のように実装できます。

step13/add.c

static float* add_backward(Function* const p_self, const float* gys) {
  float* gxs;
  gxs = (float*)malloc(p_self->input_num * sizeof(float));
  gxs[0] = gys[0];  gxs[1] = gys[0];
  return gxs;
}

またVariable_backwardを修正する必要があります。Functionの入出力の型がVariable**になったためです。次のようにします。

step13/variable.c

void Variable_backward(Variable* p_self) {
  if (!(p_self->creator_exists)) {
    return;
  }
  if (!(p_self->grad_exists)) {
    p_self->grad = 1.0;
    p_self->grad_exists = TRUE;
  }
  FunctionStack fs;
  FunctionStack_init(&fs, FUNCTION_STACK_SIZE);
  FunctionStack_push(&fs, *(p_self->p_creator));
  Function f;
  while (fs.last != -1) {
    f = FunctionStack_pop(&fs);
    float* gys;
    gys = (float*)malloc(f.output_num * sizeof(float));
    for (int i = 0; i < f.output_num; i++) {
      gys[i] = f.p_io[1][i]->grad;
    }
    float* gxs;
    gxs = (float*)malloc(f.input_num * sizeof(float));
    gxs = f.p_methods->backward(&f, gys);

    for (int i = 0; i < f.input_num; i++) {
      f.p_io[0][i]->grad = gxs[i];
      if (f.p_io[0][i]->creator_exists) {
        FunctionStack_push(&fs, *(f.p_io[0][i]->p_creator));
      }
    }
  }
}

だいぶ複雑になったようにも見えますが、単に入出力が複数になったことに対応するためにfor文を回しているだけです。ちなみに今ステップから定数系はconstant.h(constant.c)にまとめて書くことにしています。²

step13/constant.h

#ifndef _CONSTANT_H_
#define _CONSTANT_H_

enum bool {
  FALSE,
  TRUE
};

const int FUNCTION_STACK_SIZE;

#endif // _CONSTANT_H_

step13/constant.c

#include "constant.h"

const int FUNCTION_STACK_SIZE = 100;

では放置していたSquareとExpの修正も行いましょう。

step13/square.c

static float* square_forward(Function* const p_self, const float* xs) {
  float* ys;
  ys = (float*)malloc(sizeof(float));
  ys[0] = pow(xs[0], 2);
  return ys;
}

static float* square_backward(Function* const p_self, const float* gys) {
  float* gxs;
  gxs = (float*)malloc(sizeof(float));
  gxs[0] = 2 * p_self->p_io[0][0]->data * gys[0];
  return gxs;
}

step13/exp.c

static float* exp_forward(Function* const p_self, const float* xs) {
  float* ys;
  ys = (float*)malloc(sizeof(float));
  ys[0] = exp(xs[0]);
  return ys;
}

static float* exp_backward(Function* const p_self, const float* gys) {
  float* gxs;
  gxs = (float*)malloc(sizeof(float));
  gxs[0] = p_self->p_io[1][0] * gys[0];
  return gxs;
}

ここもやはり引数が配列になっているだけなので特に問題はないでしょう。³ ということで可変長引数の逆伝播も実装できました。実際に使ってみましょう。

step13/step13.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "square.h"
#include "variable.h"
#include "function.h"
#include "add.h"

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

Variable** add(Variable* x, Variable *y){
  Add* ad;
  ad = (Add*)malloc(sizeof(Add));
  Add_init(ad);
  return Function_call((Function*)ad, x, y);
}

int main() {
  Variable x0;
  Variable x1;
  Variable_init(&x0, 2.0);
  Variable_init(&x1, 3.0);

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(square(&x0)[0], square(&x1)[0]);

  Variable_backward(y[0]);

  printf("%.2f\n", y[0]->data); // 13.00
  printf("%.2f\n", x0.grad); // 4.00
  printf("%.2f\n", x1.grad); // 6.00

  return 0;
}

上では$(x_0, x_1)=(2.0, 3.0)$における$f(x_0, x_1) = x_0^2 + x_1^2$の値およびにおける偏微分を計算していて、$\frac{\partial{f}}{\partial{x_0}}=2x_0$、$\frac{\partial{f}}{\partial{x_1}}=2x_1$ であるので、正しい結果が得られています。複数の入出力に対応した自動微分がしっかりできていそうですね。⁴ステップ13はこれで終わりです。

ステップ14 同じ変数を繰り返し使う

さて、前回はめでたく自動微分ができるようになり、また前ステップでは複数の入出力に対応したわけですが、実はいくつか問題が残っています。次のコードを見てください。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "square.h"
#include "variable.h"
#include "function.h"
#include "add.h"

Variable** add(Variable* x, Variable *y){
  Add* ad;
  ad = (Add*)malloc(sizeof(Add));
  Add_init(ad);
  return Function_call((Function*)ad, x, y);
}

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

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(&x, &x);

  Variable_backward(y[0]);

  printf("x.grad %.2f\n", x.grad); // 1.00!

  return 0;
}

$y=x+x=2x$の$x=3.0$における微分の値を求めるプログラムになっていそうですが、正しい値である2.00にならず、1.00と出力されています。なぜかというと、次の部分に問題があるからです。

step13/variable.c

void Variable_backward(Variable* p_self) {
  // 省略
  while (fs.last != -1) {
    // 省略
    gxs = f.p_methods->backward(&f, gys);

    for (int i = 0; i < f.input_num; i++) {
      f.p_io[0][i]->grad = gxs[i]; // 間違い！
      if (f.p_io[0][i]->creator_exists) {
        FunctionStack_push(&fs, *(f.p_io[0][i]->p_creator));
      }
    }
  }
}

実際このようなコードでは、ある変数に対する逆伝播の値は、同じ変数に対する逆伝播によって上書きされてしまいます。正しくは和を求めるべきなのです。したがって、まずVariableにgrad_existsというメンバ変数を追加したうえで、Variable_backwardを次のように変更します。

step14/variable.h

typedef struct variable {
  float data;
  float grad;
  bool creator_exists;
  bool grad_exists; // 追加
  Function* p_creator;
} Variable;

step14/variable.c

void Variable_backward(Variable* p_self) {
  if (p_self->creator_exists == FALSE) {
    return;
  }
  if (p_self->grad_exists == FALSE) {
    p_self->grad = 1.0;
    p_self->grad_exists = TRUE;
  }
  // 省略
  while (fs.last != -1) {
    // 省略
    gxs = f.p_methods->backward(&f, gys);

    for (int i = 0; i < f.input_num; i++) {
      if (f.p_io[0][i]->grad_exists) {
        f.p_io[0][i]->grad += gxs[i];
      } else {
        f.p_io[0][i]->grad = gxs[i];
        f.p_io[0][i]->grad_exists = TRUE;
      }
      if (f.p_io[0][i]->creator_exists) {
        FunctionStack_push(&fs, *(f.p_io[0][i]->p_creator));
      }
    }
  }
}

これで、既に勾配が存在する変数に対してさらに逆伝播で勾配が伝わってきた際には、それらが加算されるようになりました。⁵ これで上のような逆伝播を行うと微分の値は2.00と正しく求まります。
さて、上記のような変更を行った場合にはある注意点が出てきます。それは次のような場面です。

// 省略
int main() {
  Variable x;
  Variable_init(&x, 3.0);

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(&x, &x);

  Variable_backward(y[0]);

  printf("x.grad %.2f\n", x.grad); // 2.00（OK）

  Variable** z;
  z = (Variable**)malloc(sizeof(Variable*));
  z = add(&x, add(&x, &x)[0]);

  Variable_backward(z[0]); // 5.00（？）

  printf("x.grad %.2f\n", x.grad);

  return 0;
}

ここでは$x=3.0$における$y=2x$という計算の微分値を求めたあとに、同じ$x$を使って$z=3x$の微分値を求めようとしています。後者は正しい値は5.00です。しかし2回目の逆伝播のとき既に$x$には勾配があるので、先ほどの実装のせいで（？）$x$の勾配は加算され5.00という結果になってしまいます。「同じ変数を使わなきゃいいのでは？」と思うかもしれませんが、後々（はるか先？）ニューラルネットワークなどの計算を行う場合にはメモリの節約というのがかなり大事になってくるので、同じ変数を使い回すことができてほしいわけです。そこで勾配の値をリセットするためにVariable_cleargradを実装します。

step14/variable.c

void Variable_cleargrad(Variable* p_self) {
  p_self->grad = 0.0;
  p_self->grad_exists = FALSE;
}

やっていることは単純ですね。ではこれを使ってみましょう。

step14/step14.c

// 省略
int main() {
  Variable x;
  Variable_init(&x, 3.0);

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(&x, &x);

  Variable_backward(y[0]);

  printf("x.grad %.2f\n", x.grad); // 2.00（OK）

  Variable_cleargrad(&x);

  Variable** z;
  z = (Variable**)malloc(sizeof(Variable*));
  z = add(&x, add(&x, &x)[0]);

  Variable_backward(z[0]); // 3.00（OK）

  printf("x.grad %.2f\n", x.grad);

  return 0;
}

2回の逆伝播の間でVariable_cleargrad(&x)とすることで、正しく微分値を求めることができました。
というわけで自動微分の問題点を1つ解決することができました。しかしまだもっと深刻な問題点が残っているので、それをこれから解決していきます。

ステップ15 複雑な計算グラフ（理論編）

これまで扱ってきた計算グラフは、トポロジー的にやや単純な構造（一直線など）をしていました。しかし、次のような少し複雑な計算グラフを扱う必要が当然ながらこの先出てきます。

現在のFunctionおよびVariableを用いれば、上記のような計算グラフにおける順伝播はちゃんと計算できますが、逆伝播は正しく行えません。なぜかというとVariable_backwardにおいて逆伝播の順番はFunctionStack、つまりLIFOのデータ構造に依存しているからです。正しい順番は、上のグラフでいえばD、B、C、Aの順で逆伝播すべきです（あるいはD、C、B、A）。しかしFunctionStackを使っている現在の実装ではD、B、A、C、Aのような順序になってしまいます。これでは正しく逆伝播ができず、ひいては微分の値も間違えることになってしまいます。
ではどのような実装をすればよいのかというと、それぞれのFunction（およびVariable）に逆伝播における適切な優先順位をつければよさそうです。それは順伝播において行われた計算の順序を逆にすればよいので、順伝播で何番目に計算されたかという情報を「世代」として持つことが出来れば、そこそこ簡単に実装できます。
なおゼロつくでは図などを用いてもっとわかりやすく解説されているので詳しく知りたい方はぜひそちらを見てみてください。

ステップ16 複雑な計算グラフ（実装編）

では前ステップでの議論を実装してみましょう。まずVariable構造体のメンバ変数にgenerationを追加します。

step16/step16.h

typedef struct variable {
  float data;
  float grad;
  int num;
  bool creator_exists;
  bool grad_exists;
  int generation; // 追加
  Function* p_creator;
} Variable;

そしてgenerationを扱うために以下のように変更します。

step16/variable.c

void Variable_init(Variable* p_self, const float data) {
  // 省略
  p_self->creator_exists = FALSE;
  p_self->grad_exists = FALSE;
  p_self->generation = 0; // 追記
}

void Variable_set_creator(Variable* p_self, Function* func) {
  p_self->p_creator = func;
  p_self->creator_exists = TRUE;
  p_self->generation = func->generation + 1; // 追記
}

このようにすることで、Variable_initしたVariableの世代は0に設定され、またクリエイターをセットしたときにはそのクリエイター（Function）の世代に1を足したものをVariableの世代とすることにします。
続いてFunctionにもgenerationを追加し、さらにその値としては入力変数の世代の最大値（入力が複数の場合はこのようにするのが適切です）を入れるようにします。

step16/function.h

typedef struct function {
  const struct functionmethods *p_methods;
  Variable ***p_io;
  int input_num;
  int output_num;
  int generation; // 追加
} Function;

step16/function.c

 va_list va_ptr;
  va_start(va_ptr, p_self);
  float* xs;
  xs = (float*)malloc(p_self->input_num * sizeof(float));
  int max_gen = 0;

  for (int i = 0; i < p_self->input_num; i++) {
    p_self->p_io[0][i] = va_arg(va_ptr, Variable*);
    xs[i] = p_self->p_io[0][i]->data;
    if (max_gen < p_self->p_io[0][i]->generation) {
      max_gen = p_self->p_io[0][i]->generation;
    }
  }
  p_self->generation = max_gen;
  va_end(va_ptr);
  // 省略

では実際に逆伝播を行ってみましょう。一旦参考のためPythonの場合の実装例を見てみましょう。

step16.py

class Variable:
    # 省略
    def backward(self):
        if self.grad is None:
            self.grad = np.ones_like(self.data)

        funcs = []
        seen_set = set()

        def add_func(f):
            if f not in seen_set:
                funcs.append(f)
                seen_set.add(f)
                funcs.sort(key=lambda x: x.generation)

        add_func(self.creator)

        while funcs:
            f = funcs.pop()
            gys = [output.grad for output in f.outputs]
            gxs = f.backward(*gys)
            if not isinstance(gxs, tuple):
                gxs = (gxs,)

            for x, gx in zip(f.inputs, gxs):
                if x.grad is None:
                    x.grad = gx
                else:
                    x.grad = x.grad + gx

                if x.creator is not None:
                    add_func(x.creator)

このように、クリエイターをfuncsに追加するたびにfuncsを世代によってソートし、Functionを取り出す際には最後尾から取り出すことで、常に世代が大きい方から取り出すことができます。またseen_setというsetを用いることで同じ関数に対して複数逆伝播が行われることを防いでいます。
ところで書籍にもある通り、generationが最大のものを取り出しさえすればよいわけなので、ソートをする必要はありません。優先度つきキュー（ヒープ）を使えば十分です。そこでFunctionHeapを実装しました。

step16/fuctionheap.h

#include "function.h"

#ifndef _FUNCTIONHEAP_H_
#define _FUNCTIONHEAP_H_

typedef struct functionheap {
  Function* heap;
  int last;
  int size;
} FunctionHeap;

void FunctionHeap_init(FunctionHeap *p_self, const int size);
void FunctionHeap_insert(FunctionHeap *p_self, Function func);
Function FunctionHeap_deletemax(FunctionHeap *p_self);
void FunctionHeap_destroy(FunctionHeap* p_self);

#endif // _FUNCTIONHAP_H_

具体的な実装については省略します（ただのヒープなので）。ではこれを用いてVariable_backwardを修正します。

step16/variable.c

void Variable_backward(Variable* p_self) {
  if (p_self->creator_exists == FALSE) {
    return;
  }
  if (p_self->grad_exists == FALSE) {
    for (int i = 0; i < p_self->p_grad->size; i++)
      p_self->p_grad->array[i] = 1.0;
    p_self->grad_exists = TRUE;
  }
  FunctionHeap fh;
  FunctionHeap_init(&fh, FUNCTION_HEAP_SIZE);
  FunctionHeap_insert(&fh, *(p_self->p_creator));
  Function f;
  Function** seen_fs;
  seen_fs = (Function**)malloc(FUNCTION_HEAP_SIZE * sizeof(Function*));
  seen_fs[0] = p_self->p_creator;
  int cnt = 1;
  while (fh.last != -1) {
    f = FunctionHeap_deletemax(&fh);
    Ndarray* gys;
    gys = (float*)malloc(f.output_num * sizeof(float));
    for (int i = 0; i < f.output_num; i++) {
      gys[i] = f.p_io[1][i]->p_grad;
    }
    float* gxs;
    gxs = (float*)malloc(f.input_num * sizeof(float));
    gxs = f.p_methods->backward(&f, gys);

    for (int i = 0; i < f.input_num; i++) {
      if (f.p_io[0][i]->grad_exists) {
        f.p_io[0][i]->grad += gxs[i];
      } else {
        f.p_io[0][i]->grad = gxs[i];
        f.p_io[0][i]->grad_exists = TRUE;
      }
      if (f.p_io[0][i]->creator_exists) {
        bool found = FALSE;
        for (int j = 0; j < cnt; j++) {
          if (seen_fs[j] == f.p_io[0][i]->p_creator) {
            found = TRUE;
            break;
          }
        }
        if (!found) {
          FunctionHeap_insert(&fh, *(f.p_io[0][i]->p_creator));
          seen_fs[cnt++] = f.p_io[0][i]->p_creator;
        }
      }
    }
    free(gxs);
    free(gys);
  }
  free(seen_fs);
  FunctionHeap_destroy(&fh);
}

こんな感じで、（ソートの代わりにヒープを用いて）先ほどのPythonの実装と同様のことができました。なおCにはsetがないので、代わりにFunctionへのポインタの配列を使いました。
これで上手くいくはずです！　では次のような計算グラフの逆伝播を行ってみましょう。

step16/step16.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "square.h"
#include "variable.h"
#include "function.h"
#include "add.h"

Variable** add(Variable* x, Variable *y){
  Add* ad;
  ad = (Add*)malloc(sizeof(Add));
  Add_init(ad);
  return Function_call((Function*)ad, x, y);
}

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

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

  Variable** a;
  a = (Variable**)malloc(sizeof(Variable*));
  a = square(&x);

  Variable** y;
  y = (Variable**)malloc(sizeof(Variable*));
  y = add(square(a[0])[0], square(a[0])[0]);

  Variable_backward(y[0]);

  printf("y.data %.2f\n", y[0]->data); // 32.00
  printf("x.grad %.2f\n", x.grad); // 64.00

  return 0;
}

上の計算グラフが表しているのは$y=2x^4$です。このとき$y'=8x^3$となるので、$x=2.0$のときは$y=32.00$、$y'=64.00$となるはずです。実際上のコードを実行するとそのような結果が得られます。ちゃんと逆伝播できています！これで上のような少し複雑な計算グラフのみならず、どんなに複雑な場合においても正しく微分が自動で求まります。⁶ ここまできてようやく自動微分の完成といえるでしょう。これでステップ16は終わりです。

ステップ11～16総括

ちゃんとCでまともな自動微分が実装できました、嬉しいです。もう少し進める予定でしたが割と長くなってきたので今回はここまでにします（~~アドカレとしても遅刻してますし~~）⁷。結局ndarrayとメモリ管理(freeなど)がまだできませんでした...。次回こそはndarrayもどきを作ることから始めるつもりです。ここまで読んでいただきありがとうございました。
みなさまよきクリスマスをお過ごしください。メリークリスマス。

うすうす気づいた方もいるかもしれませんが、可変長引数を使った場合コンパイル時の型チェックが（通常は）できなくなります。そういう意味では既に型の安全性をある程度捨てているのですが、このあたりは使いやすさとのトレードオフで考える必要がありそうです。 ↩
Cのenum使うことって個人的にあまりないんですがbool使うときとかはこう書くと便利ですよね。ちなみにstdbool.hに_Bool型があります。また、C23では標準でbool型が搭載されるみたいです。 ↩
細かいことではありますが、最近指数関数の計算は比較的重い処理であると聞いたのでexp_backwardの実装においては出力の値を用いることにしました。 ↩
まあ複数の出力をもつ関数は当分出てきませんが。 ↩
PythonではgradがNoneかどうかで判定できます。 ↩
無論ですがそれぞれの細かい計算(Squareなど)のbackwardが正しく実装されている必要があります。 ↩
ところで記事を書くよりもコード書く方が楽しいせいで実装の方が既にステップ45まで進んでいます。 ↩

第2ステージ 自然なコードで表現する

ステップ11 可変長の引数（順伝播編）

ステップ12 可変長の引数（改善編）

ステップ13 可変長の引数（逆伝播編）

ステップ14 同じ変数を繰り返し使う

ステップ15 複雑な計算グラフ（理論編）

ステップ16 複雑な計算グラフ（実装編）

ステップ11～16総括

第2ステージ自然なコードで表現する