C言語で多層パーセプトロンを書いてみる

#はじめに
多層パーセプトロン（以下、MLPと略します）は、所謂ニューラルネットワークの1番基本となる形をしたものです。基礎知識に関する説明は下記の記事に丸投げします。
https://qiita.com/nishiy-k/items/1e795f92a99422d4ba7b
ほとんどの記事がPythonと機械学習ライブラリを用いた実装方法であり、c言語などといった低レベルな言語で書いたものは少なかったので本ページを投稿してみました。今回は隠れ層が１層のMLPを用いてXOR問題を解いてみます。

#シグモイド関数とその微分
MLPで用いる活性化関数は基本シグモイド関数です。シグモイド関数を用いる理由をシンプルに説明すると微分をするためです。微分は、ニューラルネットワークの基本となる誤差伝搬法というアルゴリズムにおいて重要なものになります。ここらへんの基礎知識は下記を参考にしてください。
https://qiita.com/Ugo-Nama/items/04814a13c9ea84978a4c
シグモイド関数と、その微分関数の実装は下記の通りです。

```c:
double sigmoid(double x) {
	return 1/(1+exp(-x));
}

//derivative of sigmoid function
double d_sigmoid(double x) {
	double a = 0.1;
	return a*x*(1-x);
}
```

#入力ベクトルとその重み付け
今回、学習データと教師データは以下のように定義しました

```c:
double train_x[4][NUM_INPUT+1] = {{0, 0, -1},{0, 1, -1},{1, 0, -1},{1, 1, -1}};
double d[4] = {0, 1, 1, 0};
```

NUM_INPUTに入力ノード数を定義し、それにバイアス-1を加えています。また、wという重みベクトルを定義して入力ベクトルとの内積を下記のように計算します。

```c:
for(i=0; i<NUM_INPUT+1; i++) {
    dot += train_x[i] * w[i];
}
```

#重みの更新
ニューラルネットワークにおいての学習とは、重みベクトルを更新させるために行うものですが、その重みを更新させるには出力データが教師データとどのくらい誤差があるのかを各層に伝える必要があります。そのために行うのが誤差伝搬法という学習アルゴリズムです。誤差とは単純に、出力データから教師データを引いた数です。
細かい理屈や理論は抜きにして、出力層から隠れ層までの重みの更新は、

```c:
for(i=0; i<NUM_HIDDEN+1; i++) {
    v[i] = v[i] - eta * y[i] * d_sigmoid(z) * (z - d);
}
```
という形で行います。NUM_HIDDENとは隠れ層のノード数で、zとdは出力層データと教師データで、etaは学習率です。また、隠れ層から入力層への重みの更新は、

```c:
for(j=0; j<NUM_INPUT+1; j++) {
    for(i=0; i<NUM_HIDDEN+1; i++) {
        w[i][j] = w[i][j] - eta * train_x[j] * d_sigmoid(y[i]) * d_sigmoid(z) * (z - d) * v[i];
    }
}
```

という形で行います。先ほどは、出力層のノード数が1個だったのでfor分が1回だけでしたが、他2層は2個以上のノード数を持っているので、2回分のfor文を回しています。
ただ、どこの層の式にも必ずetaがかかっていることがわかるかと思います。つまりこの学習率の大きさで重みの更新の変動の大きさが変わってくるわけです。

#誤差伝搬法によるXOR問題の学習の実装
実際に上記と誤差伝搬法を加えて実装してみたものが下記になります。XORとは日本語では、排他的論理和と呼ばれてます。詳しい説明は省きますが、要は入力に{01, 10}を入れたときは出力が1に、{00, 11}を入れたときは出力が0になってほしいわけです。実はこの問題、普通のパーセプトロンでは解けません。理由は非線形だからです。これが人工知能の第1の冬の時代を迎える原因になったそうですが（違ったらごめんなさい）、そこらへんは詳しい内容が書いたページでも見てください。
コードの説明にうつります。初期値として学習率は0.1、学習回数は1000000にして、重みには最初に適当な乱数を入れています。所謂ディープラーニングではないので、学習にそこまで時間はかかりません。

```c:mlp.c
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
//num of units
#define NUM_INPUT 2
#define NUM_HIDDEN 20

double sigmoid(double x) {
	return 1/(1+exp(-x));
}

//derivative of sigmoid function
double d_sigmoid(double x) {
	double a = 0.1;
	return a*x*(1-x);
}

int main(void) {
	srand((unsigned)time(NULL));
//train data
	double train_x[4][NUM_INPUT+1] = {{0, 0, -1},{0, 1, -1},{1, 0, -1},{1, 1, -1}};
	double d[4] = {0, 1, 1, 0};
//net
	double w[NUM_HIDDEN+1][NUM_INPUT+1];
	double v[NUM_HIDDEN+1];
	double y[4][NUM_HIDDEN+1];
	double z[4];
	double eta = 0.1;
	int epoch = 1000000;
//other
	int i, j, k, l;
	double tmp = 0;

//update weights using rand()
	for(l=0; l<NUM_HIDDEN+1; l++) {
		for(i=0; i<NUM_INPUT+1; i++) {
			w[l][i] = ((double)rand() / ((double)RAND_MAX + 1));
		}
	}
	for(i=0; i<NUM_HIDDEN+1; i++) {
		v[i] = ((double)rand() / ((double)RAND_MAX + 1));
	}

//tain
	for(k=0; k<epoch; k++) {
		//feedforward
		for(j=0; j<4; j++) {
			//hidden
			for(l=0; l<NUM_HIDDEN; l++) {
				for(i=0; i<NUM_INPUT+1; i++) {
					tmp += train_x[j][i] * w[l][i];
				}
				y[j][l] = sigmoid(tmp);
				tmp = 0;
			}
			y[j][NUM_HIDDEN] = -1;
			//output
			for(i=0; i<NUM_HIDDEN+1; i++) {
				tmp += y[j][i] * v[i];
			}
			z[j] = sigmoid(tmp);
			tmp = 0;
		
		//backward
			//output
			for(i=0; i<NUM_HIDDEN+1; i++) {
				v[i] = v[i] - eta * y[j][i] * d_sigmoid(z[j]) * (z[j] - d[j]);
			}
			
			//hidden
			for(l=0; l<NUM_INPUT+1; l++) {
				for(i=0; i<NUM_HIDDEN+1; i++) {
					w[i][l] = w[i][l] - eta * train_x[j][l] * d_sigmoid(y[j][i]) * d_sigmoid(z[j]) * (z[j] - d[j]) * v[i];
				}
			}
		}
		//print detail
		printf("z=");
		for(i=0; i<4; i++) {
			printf("%f ", z[i]);
		}
		printf("epoch:%d\n",k);
	}

//predict
	for(j=0; j<4; j++) {
		//hidden
		for(l=0; l<NUM_HIDDEN; l++) {
			for(i=0; i<NUM_INPUT+1; i++) {
				tmp += train_x[j][i] * w[l][i];
			}
			y[j][l] = sigmoid(tmp);
			tmp = 0;
		}
		y[j][NUM_HIDDEN] = -1;
		//output
		for(i=0; i<NUM_HIDDEN+1; i++) {
			tmp += y[j][i] * v[i];
		}
		z[j] = sigmoid(tmp);
		tmp = 0;
	}

//print result
	printf("z=");
	for(i=0; i<4; i++) {
		printf("%f ", z[i]);
	}
	printf("epoch:%d\n",k);
	return 0;
}
```

#結果
最終的に得られた結果は、

```
z=0.014707 0.984199 0.984157 0.017386 epoch:1000000
```

でした。まあ、教師データの{0, 1, 1, 0}に近い値になったのではないでしょうか。