いくぷろ〜育児プログラマーの日記

プログラミングをしながら子供を見守る日々(記載事項に間違いがあればご指摘ください)

ニューラルネットワークの層を抽象化しよう[Arrayfire]

by ·

六花です。

以前、三層のニューラルネットワークをC++のライブラリであるArrayfireとEigenを使って実装しましたが、今回はそれを改善しようと思います。

3層のニューラルネットワークの実装
https://komochiduki.net/ikupro/2019/12/30/195/

モデルを構築する際、各層をコンテナとして扱えるようにするとソースコードがすっきりしてわかりやすくなります。
※ソースコード全文は最後に載せます。

§環境
windows 10
Microsoft Visual Studio Community 2022 (64 ビット) Version 17.3.4
ArrayFire v3.8.2


struct Layer
{
    virtual void init() {}          // 勾配の初期化に用いる
    virtual void forward() = 0;     // 順伝播
    virtual void backward() = 0;    // 逆伝播
    virtual void SGD(){}            // 重みの更新
};

まず、コンテナとして扱うために抽象クラスを作ります。(構造体として定義していますが、メンバのデフォルトアクセス権がpublicであるだけなので実質クラスです。)
学習状態の保存と読み込みなど、機能が足りないかもしれませんが、とりあえずこれで最後まで作ります。
これがstd::vectorなりに入って、順番にforward()やbackward()が呼ばれることになります。
SGD()を始めとした最適化アルゴリズムについては、別途抽象クラスを作って引数として渡すことで処理の一元化ができますが、今回は割愛します。
全結合層は以下のように記述します。


struct FC_layer : public Layer
{
    const af::array& x;
    af::array& y;
    af::array& dx;
    const af::array& dy;

    const int size_batch;

    af::array W;
    af::array B;

    af::array dW;
    af::array dB;

    FC_layer(const af::array& x, af::array& y, af::array& dx, const af::array& dy)
        : x(x)
        , y(y)
        , dx(dx)
        , dy(dy)
        , size_batch(x.dims(1))
    {
        // dimsでその次元方向の要素数を得られる
        // Heの初期値
        W = af::randn(y.dims(0), x.dims(0), 1, dtype_t) * std::sqrt(var_t(2) / var_t(x.dims(0)));
        B = af::randn(y.dims(0), 1, 1, dtype_t) * std::sqrt(var_t(2) / var_t(x.dims(0)));

        dW = W * 0.0;
        dB = B * 0.0;
    }

    virtual void init()
    {
        dW = 0.0;
        dB = 0.0;
    }

    virtual void forward()
    {
        // af::tileは各次元の方向に整数倍コピーする
        y = af::matmul(W, x) + af::tile(B, 1, size_batch);
        y.eval();
    }

    virtual void backward()
    {
        dy.eval();

        // af::matmulNTは2番目の引数の行列を転置して積を求める
        // ミニバッチ分の勾配が全て集計されるので、ミニバッチ数が大きくなるとdW・dBに入る値も大きくなる
        // 今か後か、ミニバッチ数で割ることが必要
        dW += af::matmulNT(dy, x) / (var_t)size_batch;
        dB += af::sum(dy, 1) / (var_t)size_batch;

        dx += af::matmulTN(W, dy);

        dW.eval();
        dB.eval();
    }

    virtual void SGD()
    {
        constexpr var_t alpha = (var_t)0.01;
        W -= alpha * dW;
        B -= alpha * dB;

        W.eval();
        B.eval();
    }
};

最初のメンバであるx,y,dx,dyは順伝播された評価値と逆伝播された勾配です。
これらは各層で共有されるので、クラス外で定義されて、層のコンストラクタで「登録」します。
y方向(二番目の次元)は「1」で表現されることが多いですが、行列で同時に計算する都合上、その次元にはミニバッチ数を割り当てます。
eval()はその時点で行列を評価するための関数なので、必要ではないと思います。(一応その後行列の更新が必要なくなったタイミングで入れています。)
forwardで代入(=)しない場合、勾配と同じように初期化が必要です。
backwardで加算(+=)しているのは、Shortcut_Layerを作るときに必要になるからです。


struct ReLU_layer : public Layer
{
    const af::array& x;
    af::array& y;
    af::array& dx;
    const af::array& dy;

    ReLU_layer(const af::array& x, af::array& y, af::array& dx, const af::array& dy)
        : x(x)
        , y(y)
        , dx(dx)
        , dy(dy)
    {
    }

    virtual void forward()
    {
        y = af::select((x >= 0.0), x, 0.0); // 要素ごとに、x>=0ならばx,そうでないなら0.0
        y.eval();
    }

    virtual void backward()
    {
        dy.eval();
        dx += af::select((x >= 0.0), dy, 0.0); // 要素ごとに、x>=0ならばdy,そうでないなら0.0
    }
};

活性化関数のReLUを実行する層です。
ReLUはGPUの苦手とする場合分けを必要としますが、af::selectで分岐処理を実装できます。


int main()
{
    af::info();

    constexpr int size_input = 2;
    constexpr int size_hidden = 300;
    constexpr int size_output = 1;
    constexpr int size_batch = 128;

    constexpr int size_data = 100000;
    constexpr int size_data_train = (int)(size_data * 0.7);

    // Sphere function
    af::array input = af::randu(size_input, size_data, dtype_t) * (var_t)2 - (var_t)1;
    af::array output = af::sum((input * input), 0);


    // モデル構築

    std::vector<af::array> data; // 評価値
    data.push_back(af::constant(0.0, size_input, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_output, size_batch, dtype_t));

    std::vector<af::array> grad; // 誤差(傾き)
    grad.push_back(af::constant(0.0, size_input, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_output, size_batch, dtype_t));

    std::vector<std::shared_ptr<Layer>> layer; // 処理層
    layer.push_back(std::make_shared<FC_layer>(data.at(0), data.at(1), grad.at(0), grad.at(1)));
    layer.push_back(std::make_shared<ReLU_layer>(data.at(1), data.at(2), grad.at(1), grad.at(2)));
    layer.push_back(std::make_shared<FC_layer>(data.at(2), data.at(3), grad.at(2), grad.at(3)));

    // 続く

main関数に入ります。
今回の練習用の関数として、Sphere functionを選びました。f(x1, x2, x3, …)におけるすべてのxnの二乗を足し合わせたものです。
モデル構築について、今回もいわゆる三層ニューラルネットワークですが、活性化関数も一層と数えているので実装上は四層になります。
抽象クラスをそろえたことによって、FC_layerもReLU_layerも同じstd::vectorに入れることができました。
これで、後はfor文やfor_each文で順番に呼べば処理をきちんと行ってくれます。


    // 続きから

    // gradient checking
    {
        // 入力データの設定まで
        auto init_phase = [&]()
        {
            // 誤差(傾き)の初期化
            for (auto& itm : grad) { itm = 0.0; }
            for (auto& itm : layer) { itm->init(); }

            // 入力データの設定
            data.front() = input(af::span, af::seq(0, size_batch - 1));
        };

        // 順伝播と逆伝播
        auto forback = [&]()
        {
            // 順伝播
            for (auto& itm : layer) { itm->forward(); }

            // 誤差の計算
            grad.back() = data.back() - output(af::span, af::seq(0, size_batch - 1)); // y - t

            //逆伝播
            std::for_each(layer.rbegin(), layer.rend(), [](auto& itm) { itm->backward(); });
        };

        var_t grad_test = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 順伝播と逆伝播
            forback();

            grad_test = af::sum<var_t>(grad.front()(0, 0));
        }

        var_t grad_plus = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 誤差を加える
            data.front()(0, 0) += eps;

            // 順伝播と逆伝播
            forback();

            grad_plus = af::sum<var_t>(af::pow(grad.back(), 2.0)(0, 0) / (var_t)2); // mse
        }

        var_t grad_minus = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 誤差を加える
            data.front()(0, 0) -= eps;

            // 順伝播と逆伝播
            forback();

            grad_minus = af::sum<var_t>(af::pow(grad.back(), 2.0)(0, 0) / (var_t)2); // 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;
    }

    // 続く

gradient checkingで処理が正しいか確認しています。
inputの一つに着目して、普通に順伝播/逆伝播したときの勾配と、誤差を加えて順伝播したときの評価値の数値微分を比較して、両者が誤差の範囲内であるかどうかを調べます。
比較している値が違う場所なので混乱しますが、加えた「誤差」の影響が誤差の範囲で収まっていればいいようです。
都合、大体同じ処理を三度書く必要があるので、ラムダ式で括っています。
ちゃんと調べるならば、以下のように内部的に倍精度浮動小数点数を使った方がいいです。


using var_t = double;
constexpr auto dtype_t = af::dtype::f64;

なお、pow(…, 2.0)のところをpow2(…)と書いた時、計算が合わなかったのでpowを使っています。


    // 続きから

    // 学習
    int epoch = 0;
    constexpr int size_data_in_epoch = size_data_train / size_batch;
    while (true)
    {
        ++epoch;

        // ランダム選出
        // シャッフルされたインデックスを作成する
        af::array idx_data;
        {
            af::array vals_data;
            af::array sort_data = af::randu(size_data_train, 1, dtype_t);
            af::sort(vals_data, idx_data, sort_data, 0);
        }

        for (int step = 0; step < size_data_in_epoch; ++step)
        {
            // 誤差(傾き)の初期化
            for (auto& itm : grad) { itm = 0.0; }
            for (auto& itm : layer) { itm->init(); }

            // 今回のステップの学習対象
            // af::seq は範囲を指定するためのもの
            // コンストラクタによって指定される範囲が特殊なので以下を参照のこと
            // https://arrayfire.org/docs/classaf_1_1seq.htm
            af::array idx_target = idx_data(af::seq((step + 0) * size_batch, (step + 1) * size_batch - 1));

            // 入力値を設定
            data.front() = input(af::span, idx_target);

            // 順伝播
            for (auto& itm : layer) { itm->forward(); }

            // 誤差の計算
            grad.back() = data.back() - output(af::span, idx_target); // y - t

            //逆伝播
            std::for_each(layer.rbegin(), layer.rend(), [](auto& itm) { itm->backward(); });

            // 勾配の更新
            for (auto& itm : layer) { itm->SGD(); }
        }

        // 一定期間ごとにログを取る
        if (epoch % 100 == 1)
        {
            auto diff = af::mean<var_t>(af::abs(grad.back()));
            auto norm = af::norm(grad.back());
            std::cout << "epoch : " << epoch << "\t" << "diff : " << diff << "\t" << "norm : " << norm << std::endl;
        }
        else
        {
            std::cout << "epoch : " << epoch << "\r";
        }
    }

    return 0;
}

学習部分です。
問題がなければdiffやnormの値が少しずつ減っていきます。
自作の定義としてepochとstepがあります。
epochはテストデータ一周分を指します。
stepはミニバッチ一回分を指します。
epochの最初に、テストデータ分のシャッフルされたインデックスを作成します。(乱数の入った行列をソートすることでシャッフルしています。)
stepごとに、それをミニバッチ分取り出し、inputに代入します。
100epochごとにログを出していますが、それはお好みでいいと思います。

実行した結果が以下の通りです。

ArrayFire v3.8.2 (CUDA, 64-bit Windows, build 5752f2dc)
Platform: CUDA Runtime 11.4, Driver: 11070
[0] NVIDIA GeForce RTX 3080, 10240 MB, CUDA Compute 8.6
-1- NVIDIA GeForce GTX 1080 Ti, 11264 MB, CUDA Compute 6.1
--------
     3.57628
     3.17111 = 0
    0.405172
     3.57628
     3.98145
     0.79647 = 0
-------
epoch : 1       diff : 0.0608721        norm : 0.92664
epoch : 101     diff : 0.00852684       norm : 0.120807
epoch : 201     diff : 0.00659542       norm : 0.0925958
epoch : 301     diff : 0.00570944       norm : 0.0812211
epoch : 401     diff : 0.00526673       norm : 0.0735091
epoch : 501     diff : 0.00464438       norm : 0.0688159
epoch : 601     diff : 0.00434553       norm : 0.0645664
epoch : 701     diff : 0.00448307       norm : 0.06588
epoch : 801     diff : 0.00417014       norm : 0.0585569
epoch : 901     diff : 0.00465379       norm : 0.0646378
epoch : 1001    diff : 0.00367416       norm : 0.0524162
epoch : 1101    diff : 0.00412977       norm : 0.0569032
epoch : 1201    diff : 0.00342694       norm : 0.0488249
epoch : 1301    diff : 0.00338654       norm : 0.0469121
epoch : 1401    diff : 0.00352478       norm : 0.0477821
epoch : 1501    diff : 0.00339974       norm : 0.0463401
epoch : 1601    diff : 0.00346398       norm : 0.0488337
epoch : 1701    diff : 0.00274355       norm : 0.0385687
epoch : 1801    diff : 0.00294569       norm : 0.0437747

ちゃんとdiffやnormが次第に減っていっています。(難しい問題だと構築したモデルの限界で減らないこともあります)
gradient checkingが変なことになっていますが、精度に対して小さすぎる誤差を加えることによって発生します。(これは32bitに対して1e-7)
32bitならば1e-4くらいだとちゃんとした結果になります。
64bitで1e-7にした結果が以下のものです。

--------
    -1.17166
 7.57211e-09 = 0
     1.17166
     1.17166
     2.34332
 3.23136e-09 = 0
-------

これで、自力で3層ニューラルネットワークが組めるようになりました。
あとはこれの拡張で、いろんな層に変更したり、追加して層を深くしていくだけです。

最後に、全文を載せておきます。


#include <arrayfire.h>
#undef max
#undef min

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <array>
#include <vector>
#include <memory>
#include <random>

using var_t = double;
constexpr auto dtype_t = af::dtype::f64;

constexpr var_t eps = 1e-7; // floatのときは1e-4くらい、doubleのときは1e-7くらい

struct Layer
{
    virtual void init() {}          // 勾配の初期化に用いる
    virtual void forward() = 0;     // 順伝播
    virtual void backward() = 0;    // 逆伝播
    virtual void SGD(){}            // 重みの更新
};

struct FC_layer : public Layer
{
    const af::array& x;
    af::array& y;
    af::array& dx;
    const af::array& dy;

    const int size_batch;

    af::array W;
    af::array B;

    af::array dW;
    af::array dB;

    FC_layer(const af::array& x, af::array& y, af::array& dx, const af::array& dy)
        : x(x)
        , y(y)
        , dx(dx)
        , dy(dy)
        , size_batch(x.dims(1))
    {
        // dimsでその次元方向の要素数を得られる
        // Heの初期値
        W = af::randn(y.dims(0), x.dims(0), 1, dtype_t) * std::sqrt(var_t(2) / var_t(x.dims(0)));
        B = af::randn(y.dims(0), 1, 1, dtype_t) * std::sqrt(var_t(2) / var_t(x.dims(0)));

        dW = W * 0.0;
        dB = B * 0.0;
    }

    virtual void init()
    {
        dW = 0.0;
        dB = 0.0;
    }

    virtual void forward()
    {
        // af::tileは各次元の方向に整数倍コピーする
        y = af::matmul(W, x) + af::tile(B, 1, size_batch);
        y.eval();
    }

    virtual void backward()
    {
        dy.eval();

        // af::matmulNTは2番目の引数の行列を転置して積を求める
        // ミニバッチ分の勾配が全て集計されるので、ミニバッチ数が大きくなるとdW・dBに入る値も大きくなる
        // 今か後か、ミニバッチ数で割ることが必要
        dW += af::matmulNT(dy, x) / (var_t)size_batch;
        dB += af::sum(dy, 1) / (var_t)size_batch;

        dx += af::matmulTN(W, dy);

        dW.eval();
        dB.eval();
    }

    virtual void SGD()
    {
        constexpr var_t alpha = (var_t)0.01;
        W -= alpha * dW;
        B -= alpha * dB;

        W.eval();
        B.eval();
    }
};

struct ReLU_layer : public Layer
{
    const af::array& x;
    af::array& y;
    af::array& dx;
    const af::array& dy;

    ReLU_layer(const af::array& x, af::array& y, af::array& dx, const af::array& dy)
        : x(x)
        , y(y)
        , dx(dx)
        , dy(dy)
    {
    }

    virtual void forward()
    {
        y = af::select((x >= 0.0), x, 0.0); // 要素ごとに、x>=0ならばx,そうでないなら0.0
        y.eval();
    }

    virtual void backward()
    {
        dy.eval();
        dx += af::select((x >= 0.0), dy, 0.0); // 要素ごとに、x>=0ならばdy,そうでないなら0.0
    }
};

int main()
{
    af::info();

    constexpr int size_input = 2;
    constexpr int size_hidden = 300;
    constexpr int size_output = 1;
    constexpr int size_batch = 128;

    constexpr int size_data = 100000;
    constexpr int size_data_train = (int)(size_data * 0.7);

    // Sphere function
    af::array input = af::randu(size_input, size_data, dtype_t) * (var_t)2 - (var_t)1;
    af::array output = af::sum((input * input), 0);


    // モデル構築

    std::vector<af::array> data; // 評価値
    data.push_back(af::constant(0.0, size_input, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    data.push_back(af::constant(0.0, size_output, size_batch, dtype_t));

    std::vector<af::array> grad; // 誤差(傾き)
    grad.push_back(af::constant(0.0, size_input, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_hidden, size_batch, dtype_t));
    grad.push_back(af::constant(0.0, size_output, size_batch, dtype_t));

    std::vector<std::shared_ptr<Layer>> layer; // 処理層
    layer.push_back(std::make_shared<FC_layer>(data.at(0), data.at(1), grad.at(0), grad.at(1)));
    layer.push_back(std::make_shared<ReLU_layer>(data.at(1), data.at(2), grad.at(1), grad.at(2)));
    layer.push_back(std::make_shared<FC_layer>(data.at(2), data.at(3), grad.at(2), grad.at(3)));

    // gradient checking
    {
        // 入力データの設定まで
        auto init_phase = [&]()
        {
            // 誤差(傾き)の初期化
            for (auto& itm : grad) { itm = 0.0; }
            for (auto& itm : layer) { itm->init(); }

            // 入力データの設定
            data.front() = input(af::span, af::seq(0, size_batch - 1));
        };

        // 順伝播と逆伝播
        auto forback = [&]()
        {
            // 順伝播
            for (auto& itm : layer) { itm->forward(); }

            // 誤差の計算
            grad.back() = data.back() - output(af::span, af::seq(0, size_batch - 1)); // y - t

            //逆伝播
            std::for_each(layer.rbegin(), layer.rend(), [](auto& itm) { itm->backward(); });
        };

        var_t grad_test = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 順伝播と逆伝播
            forback();

            grad_test = af::sum<var_t>(grad.front()(0, 0));
        }

        var_t grad_plus = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 誤差を加える
            data.front()(0, 0) += eps;

            // 順伝播と逆伝播
            forback();

            grad_plus = af::sum<var_t>(af::pow(grad.back(), 2.0)(0, 0) / (var_t)2); // mse
        }

        var_t grad_minus = (var_t)123;
        {
            // 入力データの設定まで
            init_phase();

            // 誤差を加える
            data.front()(0, 0) -= eps;

            // 順伝播と逆伝播
            forback();

            grad_minus = af::sum<var_t>(af::pow(grad.back(), 2.0)(0, 0) / (var_t)2); // 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;
    }

    // 学習
    int epoch = 0;
    constexpr int size_data_in_epoch = size_data_train / size_batch;
    while (true)
    {
        ++epoch;

        // ランダム選出
        // シャッフルされたインデックスを作成する
        af::array idx_data;
        {
            af::array vals_data;
            af::array sort_data = af::randu(size_data_train, 1, dtype_t);
            af::sort(vals_data, idx_data, sort_data, 0);
        }

        for (int step = 0; step < size_data_in_epoch; ++step)
        {
            // 誤差(傾き)の初期化
            for (auto& itm : grad) { itm = 0.0; }
            for (auto& itm : layer) { itm->init(); }

            // 今回のステップの学習対象
            // af::seq は範囲を指定するためのもの
            // コンストラクタによって指定される範囲が特殊なので以下を参照のこと
            // https://arrayfire.org/docs/classaf_1_1seq.htm
            af::array idx_target = idx_data(af::seq((step + 0) * size_batch, (step + 1) * size_batch - 1));

            // 入力値を設定
            data.front() = input(af::span, idx_target);

            // 順伝播
            for (auto& itm : layer) { itm->forward(); }

            // 誤差の計算
            grad.back() = data.back() - output(af::span, idx_target); // y - t

            //逆伝播
            std::for_each(layer.rbegin(), layer.rend(), [](auto& itm) { itm->backward(); });

            // 勾配の更新
            for (auto& itm : layer) { itm->SGD(); }
        }

        // 一定期間ごとにログを取る
        if (epoch % 100 == 1)
        {
            auto diff = af::mean<var_t>(af::abs(grad.back()));
            auto norm = af::norm(grad.back());
            std::cout << "epoch : " << epoch << "\t" << "diff : " << diff << "\t" << "norm : " << norm << std::endl;
        }
        else
        {
            std::cout << "epoch : " << epoch << "\r";
        }
    }

    return 0;
}

タグ:

おすすめ

コメントを残す

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