libtorchでXOR (3層バックプロパゲーション)

# 概要
libtorchでXORを学習させるニューラルネットワークを構築する。構築するニューラルネットワークは3層バックプロパゲーション。XORの学習はニューラルネットワークの基礎編でよく取り上げられるのでその作り方を記載する。

今回作成するのは以下の形のニューラルネットワーク。
![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161139/a0034b12-6845-5c96-a63e-1eb13dfc12c5.png)

# テンソル
これ読んでいる人は知っていると思うが、テンソルというのはスカラー・ベクトル・行列を扱うデータ型になる。プログラムで表現するときは多次元配列で扱う。数学的にはいろいろ定義があるらしいが、libtorch（というよりニューラルネットワークライブラリ全般）では多次元配列って扱いでOK。
スカラー・ベクトル・行列って書いたがそれ以上もある。ちなみに
0階テンソル = スカラー
1階テンソル = ベクトル
2階テンソル = 行列
である。3階テンソル以上の呼び名はない（はず）。

# ほぼ自力編
とりあえず、ほぼ自力でソースコードを書いてみる。自力じゃない部分は下記の2点。
- 行列計算にテンソルを使用している
- 勾配計算に自動微分を使用している

## ソースコード

``` cpp:example-app.cpp
#include <torch/torch.h>
#include <iostream>

int main() {

	auto op = torch::TensorOptions().requires_grad( true ).dtype( torch::ScalarType::Float );
	std::pair<torch::Tensor, torch::Tensor> data[] = {
		{ torch::tensor( { 0.0, 0.0 }, op ), torch::tensor( { 0.0 }, op ) },
		{ torch::tensor( { 0.0, 1.0 }, op ), torch::tensor( { 1.0 }, op ) },
		{ torch::tensor( { 1.0, 0.0 }, op ), torch::tensor( { 1.0 }, op ) },
		{ torch::tensor( { 1.0, 1.0 }, op ), torch::tensor( { 0.0 }, op ) },
	};

	auto weight1 = torch::rand( { 3, 2 }, op );
	auto bias1 = torch::rand( { 3 }, op );
	auto weight2 = torch::rand( { 1, 3 }, op );
	auto bias2 = torch::rand( { 1 }, op );
	for( size_t epoch = 0; epoch < 100001; epoch++ ) {
		float loss_mean = 0;

		for( auto& batch : data ) {
			auto hidden1 = batch.first.matmul( weight1.t() ).add( bias1 );
			auto hidden2 = torch::sigmoid( hidden1 );
			auto prediction = hidden2.matmul( weight2.t() ).add( bias2 );

			auto loss = torch::pow( ( prediction - batch.second ), 2 ).mean();
			loss.backward();

			{
				torch::NoGradGuard noguard;
				weight1 -= 0.01 * weight1.grad();
				bias1 -= 0.01 * bias1.grad();
				weight2 -= 0.01 * weight2.grad();
				bias2 -= 0.01 * bias2.grad();
				weight1.grad().zero_();
				bias1.grad().zero_();
				weight2.grad().zero_();
				bias2.grad().zero_();
			}

			loss_mean += *loss.data<float>();
			if( epoch % 100 == 0 )
				std::cout << batch.first[0].item<float>() << " " << batch.first[1].item<float>() << " = " << prediction.item<float>() << std::endl;

		}
		if( epoch % 100 == 0 ) {
			loss_mean /= 4.0;
			std::cout << "count : " << epoch << ", loss : " << loss_mean << std::endl;
			std::cout << "----------------------------" << std::endl;
			if( loss_mean < 0.0001 ) break;
		}
	}
}
```

## ニューラルネットワークの定義と予測計算
ニューラルネットワークを構成する要素は下図の線で描いた部分（赤枠で示した部分）と、丸で描いた部分に分かれる。
![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161139/a016079e-4d04-1765-a16a-030e34773264.png)
丸で描いた部分が層（レイヤー）になるのだが、この部分をテンソルで表現する。入力層は2個であるので2次元ベクトル、中間層は3個であるので3次元ベクトル、出力層は1個であるので1次元ベクトル。

### 線形関数
ニューラルネットワークでは図で描いた線1本ごとに重みを持っている。さらに線がつながる箇所で重みを持っている。
![](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/161139/25315ace-4075-30df-5b9c-59da9082df70.png)
計算式としては以下のようになる。
$y_1 = x_1w_{11}+x_2w_{21} + b_1$
$y_2 = x_1w_{12}+x_2w_{22} + b_2$
$y_3 = x_1w_{13}+x_2w_{23} + b_3$
これを線形関数と呼ぶ。ちなみに、ニューラルネットワークでは・・・と書いたがニューラルネットワークでは線形関数以外も利用する。

線形関数を一般化すると下記の式になる。
$y = xW^T + B$
$x$: 入力
$y$: 出力
$W$: 重み ($W^T$は$W$の転置行列)
$B$: バイアス

### 線形関数の重みとバイアスの定義部分
3層ネットワークなので、接続部分は２つある。したがって、重みとバイアスは２つ定義する必要がある。
``` cpp
auto weight1 = torch::rand( { 3, 2 }, op );
auto bias1 = torch::rand( { 3 }, op );
auto weight2 = torch::rand( { 1, 3 }, op );
auto bias2 = torch::rand( { 1 }, op );
```
乱数で初期値を決めている。
weight1,2は重み。行列の数は「出力層の数 $\times$ 入力層の数」。
bias1,2はバイアス。ベクトルの数は「出力層の数」。

### ニューラルネットワークの計算部分
線形関数の重みとバイアスを定義したので今度は計算部分。
``` cpp
auto hidden1 = batch.first.matmul( weight1.t() ).add( bias1 );
auto hidden2 = torch::sigmoid( hidden1 );
auto prediction = hidden2.matmul( weight2.t() ).add( bias2 );
```
batch.firstは入力値を示している。
weight1.t()はweight1の転置行列を求めるメソッド。
matmul()は行列の掛け算。
add()はベクトルの足し算。

ちなみに下記のコードでは動かなかった。理由は、*(アスタリスク)では行列の掛け算ができなかったため。
``` cpp
auto hidden1 = batch.first * weight1.t() + bias1;
```
2つの線形関数の間にシグモイド関数(torch::sigmoid())を挟み込んでいるが、これは活性化関数と呼ばれる関数。活性化関数にはシグモイド関数以外にもいろいろな関数がある。

## 損失の計算
損失は、ニューラルネットワークで予測された値と正解の乖離を示す値となる。この値が０になると入力と出力が完全一致していることを意味するう。
損失の計算には平均２乗誤差と呼ばれる関数を利用する。

``` cpp
auto loss = torch::pow( ( prediction - batch.second ), 2 ).mean();
```

## 誤差逆伝播
勾配を計算して、重みとバイアスを修正する。

### 勾配計算
勾配計算はlibtorchが自動で行ってくれる（自動微分）。コードは下記の部分。

``` cpp
loss.backward();
```

### 重みとバイアスの修正
重みとバイアスの修正を行う。勾配値はgrad()で出力されるため、その値に学習率(0.01)を掛けて重みとバイアスから引く。
このときには自動微分機能を切る必要があるので、torch::NoGradGuradという構造体が存在するスコープで行う。
最後に、勾配値をゼロクリア(zero_())しておく必要がある。これを行わないと、勾配計算した結果が累積されてしまう。
``` cpp
{
	torch::NoGradGuard noguard;
	weight1 -= 0.01 * weight1.grad();
	bias1 -= 0.01 * bias1.grad();
	weight2 -= 0.01 * weight2.grad();
	bias2 -= 0.01 * bias2.grad();
	weight1.grad().zero_();
	bias1.grad().zero_();
	weight2.grad().zero_();
	bias2.grad().zero_();
}
```

## 結果
予測、損失計算、誤差逆伝播を繰り返すことにより、学習していく。何回か繰り返したところ、８０００～９０００回くらいで収束する。

```
0 0 = 0.00851332
0 1 = 0.986627
1 0 = 0.987648
1 1 = 0.0180201
count : 8900, loss : 0.000182157
----------------------------
0 0 = 0.00710408
0 1 = 0.988854
1 0 = 0.989711
1 1 = 0.0150119
count : 9000, loss : 0.000126481
----------------------------
0 0 = 0.00592092
0 1 = 0.990719
1 0 = 0.991436
1 1 = 0.0124931
count : 9100, loss : 8.76508e-05
----------------------------
```

# ライブラリ利用編
自力でやった部分をライブラリで行っていく。

## ソースコード
``` cpp:example-app.cpp
class Xor : public torch::data::datasets::Dataset<Xor> {
public:
	Xor() {
		auto op = torch::TensorOptions().requires_grad( true ).dtype( torch::ScalarType::Float );
		_data[0] = { torch::tensor( { 0.0, 0.0 }, op ), torch::tensor( { 0.0 }, op ) };
		_data[1] = { torch::tensor( { 0.0, 1.0 }, op ), torch::tensor( { 1.0 }, op ) };
		_data[2] = { torch::tensor( { 1.0, 0.0 }, op ), torch::tensor( { 1.0 }, op ) };
		_data[3] = { torch::tensor( { 1.0, 1.0 }, op ), torch::tensor( { 0.0 }, op ) };
	}

	virtual torch::optional<size_t> size() const override {
		return torch::optional<size_t>( 4 );
	}

	virtual ExampleType get( size_t index ) override {
		return _data[index];
	}

private:
	ExampleType _data[4];
};

struct Net : torch::nn::Module {
	Net() {
		fc1 = register_module( "fc1", torch::nn::Linear( 2, 3 ) );
		fc2 = register_module( "fc2", torch::nn::Linear( 3, 1 ) );
	}

	torch::Tensor forward( torch::Tensor x ) {
		x = fc1( x );
		x = torch::sigmoid( x );
		x = fc2( x );
		return x;
	}

	torch::nn::Linear fc1 = nullptr;
	torch::nn::Linear fc2 = nullptr;
};

int main() {
	auto net = std::make_shared<Net>();
	auto data_loader = torch::data::make_data_loader( 
		Xor().map( torch::data::transforms::Stack<>() ) );
	torch::optim::SGD optimizer( net->parameters(), 0.01 );

	for( size_t epoch = 0; epoch < 10001; epoch++ ) {
		float loss_mean = 0;

		for( auto& batch : *data_loader ) {
			optimizer.zero_grad();
			auto prediction = net->forward( batch.data );
			auto loss = torch::mse_loss( prediction, batch.target );
			loss.backward();
			optimizer.step();

			if( epoch % 100 == 0 )
				std::cout << batch.data[0][0].item<float>() << " " << batch.data[0][1].item<float>() << " = " << prediction[0][0].item<float>() << std::endl;
			loss_mean += loss.item<float>();
		}
		if( epoch % 100 == 0 ) {
			loss_mean /= 4.0;
			std::cout << "count : " << epoch << ", loss : " << loss_mean << std::endl;
			std::cout << "----------------------------" << std::endl;
			if( loss_mean < 0.0001 ) break;
		}
	}
}
```

## ニューラルネットワークの定義と予測計算
ニューラルネットワークの構築部分と予測計算を１つのクラス(サンプルでは構造体)にまとめる。

``` cpp
struct Net : torch::nn::Module {
	Net() {
		fc1 = register_module( "fc1", torch::nn::Linear( 2, 3 ) );
		fc2 = register_module( "fc2", torch::nn::Linear( 3, 1 ) );
	}

	torch::Tensor forward( torch::Tensor x ) {
		x = fc1( x );
		x = torch::sigmoid( x );
		x = fc2( x );
		return x;
	}

	torch::nn::Linear fc1 = nullptr;
	torch::nn::Linear fc2 = nullptr;
};
```
コンストラクタでニューラルネットワークで使用する線形関数を定義している。
torch::nn::Linearは線形関数の定義部分。第１引数が入力層の数。第２引数が出力層の数。
register_module()はtorch::nn::Moduleに定義されているメソッドで、線形関数に名前をつけている箇所。これがないと、勾配計算ができない模様なので必ずつける必要がある。
自力版と比較すると重みとバイアスが隠蔽されている。
``` cpp
fc1 = register_module( "fc1", torch::nn::Linear( 2, 3 ) );	// auto weight1 = torch::rand( { 3, 2 }, op );
															// auto bias1 = torch::rand( { 3 }, op );
fc2 = register_module( "fc2", torch::nn::Linear( 3, 1 ) );	// auto weight2 = torch::rand( { 1, 3 }, op );
															// auto bias2 = torch::rand( { 1 }, op );
```

forward()は予測計算を行う部分。ちなみに仮想関数ではない。
自力で行った部分とほぼ同一。計算部分は隠蔽されている。
``` cpp
x = fc1( x );				// auto hidden1 = batch.first.matmul( weight1.t() ).add( bias1 );
x = torch::sigmoid( x );	// auto hidden2 = torch::sigmoid( hidden1 );
x = fc2( x );				// auto output = hidden2.matmul( weight2.t() ).add( bias2 );
```

## 損失の計算
損失の計算に自力版では平均２乗法をそのまま記述したが、libtorchに関数が用意されているのでそれを利用する。
``` cpp
auto loss = torch::mse_loss( prediction, batch.target );	// auto loss = torch::pow( ( prediction - batch.second ), 2 ).mean();
```

## 誤差逆伝播
勾配を計算して、重みとバイアスを修正する。勾配計算は完全に同一。

### 重みとバイアスの修正
自力版では重みとバイアスの修正を自力で行っていたが、最適化関数であるSGD(確率的勾配降下法)を利用する。
最初に変更するパラメータと学習率を指定して、あとはstep()で勾配の修正、zero_grad()で勾配のゼロクリアである。
サンプルソースではzero_grad()は本家MNISTのサンプルに合わせてループの最初に移動させてある。

``` cpp
torch::optim::SGD optimizer( net->parameters(), 0.01 );
															// torch::NoGradGuard noguard;
optimizer.step();											// weight1 -= 0.01 * weight1.grad();
															// bias1 -= 0.01 * bias1.grad();
															// weight2 -= 0.01 * weight2.grad();
															// bias2 -= 0.01 * bias2.grad();
optimizer.zero_grad();										// weight1.grad().zero_();
															// bias1.grad().zero_();
															// weight2.grad().zero_();
															// bias2.grad().zero_();
```

## データの定義
データセットを作成する。MNISTデータセットはデフォルトで用意されているので、自力で組む必要がないが、XORのように用意されていないものは自力で作成する必要がある。
``` cpp
class Xor : public torch::data::datasets::Dataset<Xor> {
public:
	Xor() {
		auto op = torch::TensorOptions().requires_grad( true ).dtype( torch::ScalarType::Float );
		_data[0] = { torch::tensor( { 0.0, 0.0 }, op ), torch::tensor( { 0.0 }, op ) };
		_data[1] = { torch::tensor( { 0.0, 1.0 }, op ), torch::tensor( { 1.0 }, op ) };
		_data[2] = { torch::tensor( { 1.0, 0.0 }, op ), torch::tensor( { 1.0 }, op ) };
		_data[3] = { torch::tensor( { 1.0, 1.0 }, op ), torch::tensor( { 0.0 }, op ) };
	}

	virtual torch::optional<size_t> size() const override {
		return torch::optional<size_t>( 4 );
	}

	virtual ExampleType get( size_t index ) override {
		return _data[index];
	}

private:
	ExampleType _data[4];
};
```
実装しなければならない仮想関数はsize()とget()の２つ。まぁ、名前そのまんまの関数なので特に悩むところはないと思う。
getの返り値はExampleTypeは学習データの入力値(data)と結果値(target)の２つが入っている。

## 結果
収束回数・速度ともに自力版とほぼ同じ。異なるのは、データローダーを用いているためか、入力値の並び順がことなるところ。

```
1 1 = 0.0225483
0 1 = 0.984644
1 0 = 0.983162
0 0 = 0.010332
count : 8400, loss : 0.00028362
----------------------------
0 0 = 0.00779504
1 0 = 0.986217
0 1 = 0.988565
1 1 = 0.0191447
count : 8500, loss : 0.000187003
----------------------------
0 0 = 0.00645444
1 1 = 0.0146134
0 1 = 0.98963
1 0 = 0.988667
count : 8600, loss : 0.000122794
----------------------------
1 0 = 0.991126
0 0 = 0.00551018
1 1 = 0.0121528
0 1 = 0.991868
count : 8700, loss : 8.07313e-05
----------------------------
```

# 最後に
XORのサンプルをlibtorchで作成してみたがいかがだっただろうか。ここでは記載していないが、[libtorch本家のMNISTサンプル](https://pytorch.org/cppdocs/frontend.html)と比較してほしい。ニューラルネットワークの構造が異なるがほぼ同一のコードで実行できていることが確認できると思う。

# 参考文献
- 赤石雅典 . PyTorch & 深層学習プログラミング . 日経BP
- 巣籠悠輔 . 詳解　ディープラーニング 第2版 . マイナビ