ニューラルネットワークの検算、GradientCheckingの実装[Eigen vs Arrayfire]
六花です。
既製のライブラリを使うときには気にすることはありませんが、自作するときに確認しないといけないことがあります。
それは、「そもそもその計算が正しいのか」ということです。
何度やっても誤差や重みがinf/nanになってしてしまう場合も勿論、見かけ上は誤差が小さくなっていっていても計算が正しくない場合があります。
また、計算が正しくないと判ったとき、バグの箇所を絞り込めないと多大な労力を強いられます。
そのために、gradient checkingと呼ばれる作業を行います。
簡潔に言うと、「誤差逆伝播法で求めた勾配」と、「数値微分の結果」がほぼ等しいことを確認することです。
数値微分、ひいては浮動小数点数は常に誤差を含むので完全に一致はしませんが、一般に微小な値を1e-4～1e-7として、差分がそれより小さな値であれば良いようです。
注意点として、「微小な値」を小さくし過ぎると丸め誤差の影響で数値微分の結果が0になってしまい、正常に機能しません。
私の理解のために、以下のサイト等を参考にしました。
https://towardsdatascience.com/how-to-debug-a-neural-network-with-gradient-checking-41deec0357a9
ニューラルネットワークを「y=f(x)」とすれば、「出力層 = f(入力層)」なので、「入力層への勾配」と「入力層の一要素に微小な変化を加えた数値微分の結果」を比べればいいことになります。
ちなみに、一回の検査では一か所の勾配しか確かめられないので、完全に検査したいならば全要素数分だけ繰り返す必要があります。(その計算を一回に圧縮したのが誤差逆伝播法)
サンプルコードでは、入力の一か所の勾配のみ確かめています。
また、説明を単純化したため関数fの引数を入力層の要素だけとしましたが、「出力層 = f(入力層, 重み)」と考えることもできるので重みの勾配も同様に可能です。(前掲のリンク先では重みの勾配に微小な値を加えられています。)
誤差逆伝播法によって求めた勾配は、コスト関数から算出した誤差に対するものであることを忘れてはいけません。
コスト関数の微分が「出力層 – 教師信号」です。
ですから、数値微分の方も出力層の値そのものではなく、出力層と教師信号とコスト関数から算出した誤差を使います。
下記のサンプルコードのコスト関数はMSE(二乗誤差)です。
なお、サンプルコードの行列の宣言において、floatとdoubleを一括で変更できるように記述を変更しました。
■環境
Visual Studio Community 2019 Version 16.4.0
Eigen 3.3.7
Arrayfire 3.6.4
■サンプルコード
#include <array>
#include <random>
#include <iostream>
#include <iomanip>
#include "../Eigen/Eigen/Core"
#include <arrayfire.h>
using var_t = double;
auto dtype_t = af::dtype::f64;
int main()
{
constexpr int size_input = 3;
constexpr int size_hidden = 100;
constexpr int size_output = 1;
constexpr int size_data = 1;
// 乱数の準備
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<var_t>(-1.0, +1.0);
var_t input_src[size_data][size_input];
var_t output_src[size_data];
// y = x1 + x2+ x3
for (int i1 = 0; i1 < size_data; ++i1)
{
output_src[i1] = var_t(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::Matrix<var_t, -1, -1> input(size_input, size_data);
Eigen::Matrix<var_t, -1, -1> 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::Matrix<var_t, -1, -1> x_data(size_input, 1);
Eigen::Matrix<var_t, -1, -1> x_grad(size_input, 1);
Eigen::Matrix<var_t, -1, -1> W_data(size_hidden, size_input);
Eigen::Matrix<var_t, -1, -1> W_grad(size_hidden, size_input);
Eigen::Matrix<var_t, -1, -1> B_data(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> B_grad(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> z1_data(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> z1_grad(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> z2_data(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> z2_grad(size_hidden, 1);
Eigen::Matrix<var_t, -1, -1> V_data(size_output, size_hidden);
Eigen::Matrix<var_t, -1, -1> V_grad(size_output, size_hidden);
Eigen::Matrix<var_t, -1, -1> C_data(size_output, 1);
Eigen::Matrix<var_t, -1, -1> C_grad(size_output, 1);
Eigen::Matrix<var_t, -1, -1> y_data(size_output, 1);
Eigen::Matrix<var_t, -1, -1> y_grad(size_output, 1);
// Xavierの初期値
auto dist_W = std::normal_distribution<var_t>(0.0, 1.0 / std::sqrt(var_t(size_input)));
auto dist_V = std::normal_distribution<var_t>(0.0, 1.0 / std::sqrt(var_t(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 var_t p) { return var_t(1) / (var_t(1) + exp(-p)); };
auto df_sigmoid = [&f_sigmoid](const var_t p) { return f_sigmoid(p) * (var_t(1) - f_sigmoid(p)); };
auto dist_data = std::uniform_int_distribution<int>(0, size_data - 1);
constexpr var_t eps = 1e-4;
var_t data_base = input.coeff(0, 0);
{
// forward
x_data = input.col(0);
x_data.coeffRef(0, 0) = data_base;
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(0);
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;
}
var_t grad_test = x_grad.coeff(0, 0);
decltype(y_data) diff;
{
// forward
x_data = input.col(0);
x_data.coeffRef(0, 0) = data_base + eps;
z1_data = W_data * x_data + B_data;
z2_data = z1_data.unaryExpr(f_sigmoid);
y_data = V_data * z2_data + C_data;
// diff
diff = y_data - output.col(0);
}
auto grad_plus = (diff.cwiseAbs2() * (var_t(1) / var_t(2))).sum(); // mse
{
// forward
x_data = input.col(0);
x_data.coeffRef(0, 0) = data_base - eps;
z1_data = W_data * x_data + B_data;
z2_data = z1_data.unaryExpr(f_sigmoid);
y_data = V_data * z2_data + C_data;
// diff
diff = y_data - output.col(0);
}
auto grad_minus = (diff.cwiseAbs2() * (var_t(1) / var_t(2))).sum(); // mse
auto gradapprox = (grad_plus - grad_minus) / (2.0 * eps); // 数値微分
auto numerator = sqrt(pow(grad_test - gradapprox, 2));
auto norm_grad = sqrt(grad_test * grad_test);
auto norm_gradapprox = sqrt(gradapprox * gradapprox);
auto denominator = norm_grad + norm_gradapprox;
auto difference = numerator / denominator; // 検証すべき差分
std::cout
<< "--------" << std::endl
<< std::setw(12) << gradapprox << "" << std::endl
<< std::setw(12) << numerator << " = 0" << std::endl
<< std::setw(12) << norm_grad << std::endl
<< std::setw(12) << norm_gradapprox << "" << std::endl
<< std::setw(12) << denominator << std::endl
<< std::setw(12) << difference << " = 0 " << std::endl;
std::cout << "-------" << std::endl;
}
std::cout << std::endl;
// Arrayfire
//-------------------
{
af::info();
af::array input = af::constant(0.0, size_input, size_data, dtype_t);
af::array output = af::constant(0.0, size_output, size_data, dtype_t);
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, dtype_t);
af::array x_grad = af::constant(0.0, size_input, 1, dtype_t);
af::array W_data = af::randn(size_hidden, size_input, dtype_t) * (1.0 / std::sqrt(size_input)); // Xavierの初期値
af::array W_grad = af::constant(0.0, size_hidden, size_input, dtype_t);
af::array B_data = af::constant(0.0, size_hidden, 1, dtype_t);
af::array B_grad = af::constant(0.0, size_hidden, 1, dtype_t);
af::array z1_data = af::constant(0.0, size_hidden, 1, dtype_t);
af::array z1_grad = af::constant(0.0, size_hidden, 1, dtype_t);
af::array z2_data = af::constant(0.0, size_hidden, 1, dtype_t);
af::array z2_grad = af::constant(0.0, size_hidden, 1, dtype_t);
af::array V_data = af::randn(size_output, size_hidden, dtype_t) * (1.0 / std::sqrt(size_hidden)); // Xavierの初期値
af::array V_grad = af::constant(0.0, size_output, size_hidden, dtype_t);
af::array C_data = af::constant(0.0, size_output, 1, dtype_t);
af::array C_grad = af::constant(0.0, size_output, 1, dtype_t);
af::array y_data = af::constant(0.0, size_output, 1, dtype_t);
af::array y_grad = af::constant(0.0, size_output, 1, dtype_t);
constexpr var_t eps = 1e-4;
var_t data_base = af::sum<var_t>(input(0, 0)); // 基準値を保存する
{
// forward
x_data = input(af::span, 0);
x_data(0, 0) = data_base; // 変化を書き換え(下記の同様の分との比較のために一応記述してある)
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, 0);
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 = af::matmul(W_data.T(), z1_grad);
}
var_t grad_test = af::sum<var_t>(x_grad(0, 0));
af::array diff;
{
// forward
x_data = input(af::span, 0);
x_data(0, 0) = data_base + eps; // 変化を書き換え
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
diff = y_data - output(af::span, 0);
}
var_t grad_plus = af::sum<var_t>((diff * diff * (1.0 / 2.0))); // mse
{
// forward
x_data = input(af::span, 0);
x_data(0, 0) = data_base - eps; // 変化を書き換え
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
diff = y_data - output(af::span, 0);
}
var_t grad_minus = af::sum<var_t>((diff * diff * (1.0 / 2.0))); // mse
auto gradapprox = (grad_plus - grad_minus) / (2.0 * eps); // 数値微分
auto numerator = sqrt(pow(grad_test - gradapprox, 2));
auto norm_grad = sqrt(grad_test * grad_test);
auto norm_gradapprox = sqrt(gradapprox * gradapprox);
auto denominator = norm_grad + norm_gradapprox;
auto difference = numerator / denominator; // 検証すべき差分
std::cout
<< "--------" << std::endl
<< std::setw(12) << gradapprox << "" << std::endl
<< std::setw(12) << numerator << " = 0" << std::endl
<< std::setw(12) << norm_grad << std::endl
<< std::setw(12) << norm_gradapprox << "" << std::endl
<< std::setw(12) << denominator << std::endl
<< std::setw(12) << difference << " = 0 " << std::endl;
std::cout << "-------" << std::endl;
}
std::cout << std::endl << std::endl << "end" << std::endl;
return 0;
}
■実行結果
Eigen version : 3.3.7
--------
-0.0499762
6.71696e-11 = 0
0.0499762
0.0499762
0.0999523
6.72016e-10 = 0
-------
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
--------
-0.0754641
3.93744e-11 = 0
0.0754641
0.0754641
0.150928
2.60882e-10 = 0
-------
end