ニューラルネットワークの検算、GradientCheckingの実装[Eigen vs Arrayfire]

by · 公開済み · 更新済み

　　

六花です。

既製のライブラリを使うときには気にすることはありませんが、自作するときに確認しないといけないことがあります。
それは、「そもそもその計算が正しいのか」ということです。
何度やっても誤差や重みが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

　　

■サンプルコード

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
#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

タグ:

おすすめ

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です