3層のニューラルネットワークの実装[Eigen vs Arrayfire]
六花です。
以前、tiny-dnnについて記事を書きましたが、自分で実装してみたくなったのでこの一年間は色々と調べていました。
その結果、まったくニューラルネットワークを知らない人向けの説明、数式を用いて詳細に書いてあるもの、「ではライブラリを使いましょう」というサイトしか見当たらなかったので、自分のような初級者向けに記事を書くことにしました。
今回は、単純な3層ニューラルネットワークを行列計算ライブラリの力を借りながら実際にc++の実装例を示したいと思います。
CPU向けのEigen、GPU向けのArrayfireについて、比較できるようになるべく同じ流れで記述したいと思います。
■環境
Visual Studio Community 2019 Version 16.4.0
Eigen 3.3.7
Arrayfire 3.6.4
■共通の流れ
Eigen、Arrayfire共通の流れは以下の通りです。
・入力と出力の定義
今回は「y = x1 + x2 + x3」を満たす3入力1出力のデータセットを乱数で用意します。
・行列群の定義
ニューラルネットの本体となる行列群を定義します。順伝播に使うdata系と逆伝播に使うgrad系があります。
重みの初期値も設定します。
・順伝播
・誤差の取得
・逆伝播
・重みの更新
・(誤差が最低値を更新したらログとして表示)
・→順伝播へ
■実際のコード
#include <array>
#include <random>
#include <iostream>
#include "../Eigen/Eigen/Core"
#include <arrayfire.h>
int main()
{
constexpr int size_input = 3;
constexpr int size_hidden = 100;
constexpr int size_output = 1;
constexpr int size_data = 10000;
// 乱数の準備
std::random_device rd;
std::array<std::seed_seq::result_type, std::mt19937_64::state_size> seed_data;
std::generate(seed_data.begin(), seed_data.end(), ref(rd));
std::seed_seq s_seq(seed_data.begin(), seed_data.end());
std::mt19937_64 mt(s_seq);
for (int i = 0; i < 100; ++i) { mt(); }
auto dist_input = std::uniform_real_distribution<float>(-1.0, +1.0);
float input_src[size_data][size_input];
float output_src[size_data];
// y = x1 + x2+ x3
for (int i1 = 0; i1 < size_data; ++i1)
{
output_src[i1] = float(0);
for(int i2 = 0; i2 < size_input; ++i2)
{
input_src[i1][i2] = dist_input(mt);
output_src[i1] += input_src[i1][i2];
}
}
// Eigen
//-------------------
{
std::cout << "Eigen version : " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << std::endl;
Eigen::MatrixXf input(size_input, size_data);
Eigen::MatrixXf output(size_output, size_data);
for (int i1 = 0; i1 < size_data; ++i1)
{
for (int i2 = 0; i2 < size_input; ++i2)
{
input.coeffRef(i2, i1) = input_src[i1][i2];
}
output.coeffRef(0, i1) = output_src[i1];
}
// z1 = W * x + B
// z2 = sigmoid(z1)
// y = V * z + C
Eigen::MatrixXf x_data(size_input, 1);
Eigen::MatrixXf x_grad(size_input, 1);
Eigen::MatrixXf W_data(size_hidden, size_input);
Eigen::MatrixXf W_grad(size_hidden, size_input);
Eigen::MatrixXf B_data(size_hidden, 1);
Eigen::MatrixXf B_grad(size_hidden, 1);
Eigen::MatrixXf z1_data(size_hidden, 1);
Eigen::MatrixXf z1_grad(size_hidden, 1);
Eigen::MatrixXf z2_data(size_hidden, 1);
Eigen::MatrixXf z2_grad(size_hidden, 1);
Eigen::MatrixXf V_data(size_output, size_hidden);
Eigen::MatrixXf V_grad(size_output, size_hidden);
Eigen::MatrixXf C_data(size_output, 1);
Eigen::MatrixXf C_grad(size_output, 1);
Eigen::MatrixXf y_data(size_output, 1);
Eigen::MatrixXf y_grad(size_output, 1);
// Xavierの初期値
auto dist_W = std::normal_distribution<float>(0.0, 1.0 / std::sqrt(float(size_input)));
auto dist_V = std::normal_distribution<float>(0.0, 1.0 / std::sqrt(float(size_hidden)));
// 重みの初期化
for (int i_row = 0; i_row < W_data.rows(); ++i_row)
{
for (int i_col = 0; i_col < W_data.cols(); ++i_col)
{
W_data.coeffRef(i_row, i_col) = dist_W(mt);
}
}
for (int i_row = 0; i_row < V_data.rows(); ++i_row)
{
for (int i_col = 0; i_col < V_data.cols(); ++i_col)
{
V_data.coeffRef(i_row, i_col) = dist_V(mt);
}
}
// バイアスの初期化
B_data.setZero();
C_data.setZero();
auto f_sigmoid = [](const float p) { return float(1) / (float(1) + exp(-p)); };
auto df_sigmoid = [&f_sigmoid](const float p) { return f_sigmoid(p) * (float(1) - f_sigmoid(p)); };
auto dist_data = std::uniform_int_distribution<int>(0, size_data - 1);
std::cout << "Eigen main loop" << std::endl;
for (int i_epoch = 0; i_epoch < int(1e+5); ++i_epoch)
{
std::cout << i_epoch << "\r";
auto target_data = dist_data(mt);
// forward
x_data = input.col(target_data);
z1_data = W_data * x_data + B_data;
z2_data = z1_data.unaryExpr(f_sigmoid);
y_data = V_data * z2_data + C_data;
// diff
auto diff = (y_data - output.col(target_data)).eval();
y_grad = diff;
// backward
C_grad = y_grad;
V_grad = y_grad * z2_data.transpose();
z2_grad = V_data.transpose() * y_grad;
z1_grad = z2_grad.cwiseProduct(z1_data.unaryExpr(df_sigmoid));
B_grad = z1_grad;
W_grad = z1_grad * x_data.transpose();
x_grad = W_data.transpose() * z1_grad;
// update
constexpr float alpha = 1e-2;
C_data -= alpha * C_grad;
V_data -= alpha * V_grad;
B_data -= alpha * B_grad;
W_data -= alpha * W_grad;
const float norm_diff = diff.cwiseAbs().norm();
static float min_diff = norm_diff;
if (norm_diff < min_diff)
{
min_diff = norm_diff;
std::cout << i_epoch << "\tnorm_diff :\t" << norm_diff << std::endl;
}
}
}
std::cout << std::endl;
// Arrayfire
//-------------------
{
af::info();
// af::constantの最初の引数はテンプレートになっているので、「0」にするとint型と判断されるのでバグの元
af::array input = af::constant(0.0, size_input, size_data);
af::array output = af::constant(0.0, size_output, size_data);
for (int i1 = 0; i1 < size_data; ++i1)
{
for (int i2 = 0; i2 < size_input; ++i2)
{
input(i2, i1) = input_src[i1][i2];
}
output(0, i1) = output_src[i1];
}
// z1 = W * x + B
// z2 = sigmoid(z1)
// y = V * z + C
af::array x_data = af::constant(0.0, size_input, 1);
af::array x_grad = af::constant(0.0, size_input, 1);
af::array W_data = af::randn(size_hidden, size_input) * (1.0 / std::sqrt(size_input)); // Xavierの初期値
af::array W_grad = af::constant(0.0, size_hidden, size_input);
af::array B_data = af::constant(0.0, size_hidden, 1);
af::array B_grad = af::constant(0.0, size_hidden, 1);
af::array z1_data = af::constant(0.0, size_hidden, 1);
af::array z1_grad = af::constant(0.0, size_hidden, 1);
af::array z2_data = af::constant(0.0, size_hidden, 1);
af::array z2_grad = af::constant(0.0, size_hidden, 1);
af::array V_data = af::randn(size_output, size_hidden) * (1.0 / std::sqrt(size_hidden)); // Xavierの初期値
af::array V_grad = af::constant(0.0, size_output, size_hidden);
af::array C_data = af::constant(0.0, size_output, 1);
af::array C_grad = af::constant(0.0, size_output, 1);
af::array y_data = af::constant(0.0, size_output, 1);
af::array y_grad = af::constant(0.0, size_output, 1);
std::cout << "Arrayfire main loop" << std::endl;
for (int i_epoch = 0; i_epoch < int(1e+5) / size_data; ++i_epoch)
{
// std::vectorでいうところのstd::shuffleの代わり
// val_dataにはソート済みの行列、idx_dataにはval_dataの要素のソート前の行列の位置が入る
// 乱数で生成した行列をソートすると、idx_dataに一意の数字が入っているshuffleされた配列が手に入る
af::array idx_data;
{
af::array val_data, src_data = af::randu(size_data, 1);
af::sort(val_data, idx_data, src_data);
}
for (int i_epoch_sub = 0; i_epoch_sub < size_data; ++i_epoch_sub)
{
std::cout << i_epoch << "\t" << i_epoch_sub << "\r";
af::array idx_target = idx_data(i_epoch_sub);
// forward
x_data = input(af::span, idx_target);
z1_data = af::matmul(W_data, x_data) + B_data;
z2_data = af::sigmoid(z1_data);
y_data = af::matmul(V_data, z2_data) + C_data;
// diff
af::array diff = (y_data - output(af::span, idx_target));
y_grad = diff;
// backward
C_grad = y_grad;
V_grad = af::matmul(y_grad, z2_data.T());
z2_grad = af::matmul(V_data.T(), y_grad);
z1_grad = z2_grad * (af::sigmoid(z1_data) * (1.0 - af::sigmoid(z1_data)));
B_grad = z1_grad;
W_grad = af::matmul(z1_grad, x_data.T());
x_grad = matmul(W_data.T(), z1_grad);
// update
constexpr float alpha = 1e-2;
C_data -= alpha * C_grad;
V_data -= alpha * V_grad;
B_data -= alpha * B_grad;
W_data -= alpha * W_grad;
const float norm_diff = af::norm(diff);
static float min_diff = norm_diff;
if (norm_diff < min_diff)
{
min_diff = norm_diff;
std::cout << i_epoch << "\t" << i_epoch_sub << "\tnorm_diff :\t" << norm_diff << std::endl;
}
}
}
}
std::cout << std::endl << std::endl << "end" << std::endl;
return 0;
}
実行結果:
Eigen version : 3.3.7
Eigen main loop
1 norm_diff : 0.632244
3 norm_diff : 0.29782
6 norm_diff : 0.185658
15 norm_diff : 0.164348
17 norm_diff : 0.0715959
21 norm_diff : 0.00959095
44 norm_diff : 0.0019688
45 norm_diff : 0.000986647
430 norm_diff : 0.000916868
433 norm_diff : 0.000192344
659 norm_diff : 0.00015834
723 norm_diff : 5.28097e-05
1474 norm_diff : 2.22623e-05
1800 norm_diff : 3.63588e-06
20933 norm_diff : 1.66893e-06
22576 norm_diff : 2.98023e-07
29815 norm_diff : 1.63913e-07
99999
ArrayFire v3.6.4 (CUDA, 64-bit Windows, build 1b8030c)
Platform: CUDA Toolkit 10.0, Driver: 10020
[0] GeForce GTX 1080 Ti, 11264 MB, CUDA Compute 6.1
-1- GeForce GTX 1080, 8192 MB, CUDA Compute 6.1
Arrayfire main loop
0 6 norm_diff : 0.0608592
0 18 norm_diff : 0.00686159
0 144 norm_diff : 0.00495607
0 166 norm_diff : 0.00368699
0 244 norm_diff : 0.00269824
0 360 norm_diff : 0.000216067
0 685 norm_diff : 0.000169039
0 1123 norm_diff : 0.000109673
0 1302 norm_diff : 7.17491e-05
0 1486 norm_diff : 9.29832e-06
0 2227 norm_diff : 7.98702e-06
0 2992 norm_diff : 3.45707e-06
1 9335 norm_diff : 1.43051e-06
2 2979 norm_diff : 8.34465e-07
2 5198 norm_diff : 5.96046e-07
4 3107 norm_diff : 1.19209e-07
9 9999
end
■雑感
・メインメモリからGPUのメモリに読み込むのは非常に遅いので、多少の無理や不効率は許容してでもArrayfireの変数や関数で完結させた方が良い。
・メインメモリからGPUのメモリに読み込むのは、実際に必要になったときに行われる。
・CPUのキャッシュに収まる程度であれば、Eigenの方がよっぽど速い。