ディープじゃないディープラーニングをperlのPDLのみで超簡単実装してみた

きっかけ

Matchlabさんが書かれていたディープじゃないディープラーニングをNumPyのみで超簡単実装してみたが非常に面白かったのと、遠い異国の地でperlでのMachine Learningが少ないことを嘆き、自らMXNetのperlインタフェースをつくってCPANに上げてしまったというツワモノSergeyさんの記事Machine learning in Perlで、彼の嘆きを知ったことがきっかけです。

とりあえず、MXNetの方は別途使ってみるとして、今回はMatchlabさんの記事を参考にPerlで同じ実装をしてみました。
pythonのnumpyと同じく、perlにもPDL(Perl Data Language)という言語があり画像まわりのライブラリなどもあるのですがイマイチ流行ってないようです。perl好きとしては非常にくやしいので、Matchlabさんのpythonのコードを参考にほぼ、そのままPDLを使ってperlで実装しなおしてみました。

#PDLとは
本家のサイトとはPDLサイトにあります。そのまま、pythonのnumpyとscipyが合体したようなライブラリです。画像関連やチャートを出力するライブラリもあります。ターミナルからpdlと打つと、ipythonなんかと同じくコマンドラインでのインタラクティブな処理もできます。

$ pdl
perlDL shell v1.357
（中略）
pdl> $a = $b = pdl [[1,2],[3,4]]

pdl> print $a * $b

[
 [ 1  4]
 [ 9 16]
]

pdl> print inner $a, $b
[5 25]

pdl> print $a x $b

[
 [ 7 10]
 [15 22]
]

#インストール
自分は開発用のUbuntuにはCPANのコマンドで一発で入れられました。が、MacのNoteには権限の問題でコマンドラインではうまく入らかったため、サイトのInstall PDLから辿れる "SciPDL" binary.のリンクで、バイナリを落としてインストールしました。

ちなみに、釈迦に説法かと思いますがLinuxでのCPANのコマンドは

$ perl -MCPAN -e shell
$ cpan[1]>install PDL

と入れるだけ、です。権限をグダグダ要求されて、うまくいかない場合はforce install PDLで入ると思います。（会社のサーバーなどには使わず、自分専用のPCで自己責任でお願いします）

日本語の参考になるサイトが本当に少ないので、利用するにはPDLサイトのマニュアルを読んだ方が無難です。

#コード
Matchlabさんのコードを関数をベースに、perlで書き直しました。ただし、全く同じ関数がPDLにはなかったりした場合には同じ結果になるような形で別な方法で記載してます。
perlの思想の基本は「There's more than one way to do it」です!やり方はいくらでもあります。

インプットのirisファイルは、あらかじめGithubのirisデータから#がついたコメント行や空行を削除してCSVにしておきました。

iris.pl

#!/usr/bin/env perl
use PDL;
use PDL::Ufunc;
use PDL::NiceSlice;

$TRAIN_DATA_SIZE   = 50;    # 150個のデータのうちTRAIN_DATA_SIZE個を訓練データとして使用。残りは教師データとして使用。
$HIDDEN_LAYER_SIZE = 6;     # 中間層(隠れ層)のサイズ(今回は中間層は1層なのでスカラー)
$LEARNING_RATE     = 0.1;   # 学習率
$ITERS_NUM         = 1000;  # 繰り返し回数
$DELTA             = 0.01;

# データを読み込み
# デフォルトで'#'の行をを飛ばすようにはなってないので、一旦外の処理でコメント行を削除してから読み込み
$x     = rcols("iris.csv",[0..3], { DEFTYPE=>float, COLSEP=>"," })->transpose;
$raw_t = rcols("iris.csv",[4], { DEFTYPE=>float, COLSEP=>"," });
$onehot_t = float zeroes(3,150);

for (0 .. shape($onehot_t)->at(1)-1){
    $onehot_t($raw_t->at($_,0) , $_) .= float 1;
}

$train_x = $x(:, :$TRAIN_DATA_SIZE-1);
$train_t = $onehot_t(:, :$TRAIN_DATA_SIZE-1);
$test_x  = $x(:, $TRAIN_DATA_SIZE:);
$test_t  = $onehot_t(:, $TRAIN_DATA_SIZE:);

# 重みとバイアスの初期化
$W1 = grandom($HIDDEN_LAYER_SIZE, 4) x sqrt(2/4);
$W2 = grandom(3, $HIDDEN_LAYER_SIZE) x sqrt(2/$HIDDEN_LAYER_SIZE);
$b1 = zeroes($HIDDEN_LAYER_SIZE);
$b2 = zeroes(3);

# テストデータの結果
$test_y = forward($test_x);
print float($test_y->maximum_ind == $test_t->maximum_ind)->sum, '/', 150- $TRAIN_DATA_SIZE, "\n";

for (1 .. $ITERS_NUM){

    # 順伝搬withデータ保存
    $y1 = $train_x x $W1 + $b1;
    $y2 = relu($y1);
    $train_y = softmax($y2 x $W2 + $b2);

    # 損失関数計算
    $L = cross_entropy_error($train_y, $train_t);

    print $L,"\n" if $_ % 100 == 0;

    # 勾配計算
    # 計算グラフで求めた式を使用
    $a1 = ($train_y - $train_t) / $TRAIN_DATA_SIZE;
    $b2_gradient = dsumover $a1->transpose;
    $W2_gradient = $y2->transpose x $a1;
    $a2 = $a1 x $W2->transpose;
    $a2 = $a2 * (pdl [$y1 > 0]);

    $b1_gradient = (dsumover $a2->transpose)->clump(-1);
    $W1_gradient = $train_x->transpose x $a2->clump(1..2);

    # パラメータ更新
    $W1 = $W1  - $LEARNING_RATE * $W1_gradient;
    $W2 = $W2  - $LEARNING_RATE * $W2_gradient;
    $b1 = $b1  - $LEARNING_RATE * $b1_gradient;
    $b2 = $b2  - $LEARNING_RATE * $b2_gradient;

}

# 最終訓練データのL値
$L = cross_entropy_error(forward($train_x), $train_t);
print $L,"\n";

# テストデータの結果
$test_y = forward($test_x);
print pdl($test_y->maximum_ind == $test_t->maximum_ind)->sum, '/', 150- $TRAIN_DATA_SIZE, "\n";


# 以下、計算に必要な各種関数
sub relu{
    my $x = shift @_; 
    my $idx = $x >= 0.0;
    return $x * $idx;
}

sub softmax{
    my $x = shift @_;
    # 後ろの除算のoverflow対策で最大値を引いておく
    $x = exp $x - max($x);

    if($x->getndims == 1){
        $sum = dsumover $x;
        return $x / $sum;
    }elsif($x->getndims == 2){
        $sum = dsumover $x;
        return $x / $sum->transpose;
    }else{
        die "Dimensions are not identical!";
    }
}

sub cross_entropy_error{
    my ($y, $t) = @_;
    my $z = (shape $y)->at(1);

    if ( $y->getndims != $t->getndims ){
        die "Dimensions are not identical!";
    }elsif( $y->getndims == 1 ){
        return -($t * log($y))->sum;
    }elsif( $y->getndims == 2 ){
        return -($t * log($y))->sum / $z;
    }else{
        die "Dimension exceeds 2";
    }

}

sub forward{
    my $x = shift @_;

    my $y1 = relu($x x $W1 + $b1);
    my $y2 = softmax($y1 x $W2 + $b2);
 
    return $y2;
}

numpyとPDLの違い

いくつかの点で大きな違いがあります。まず、配列の指定の仕方の際の行列が逆、ということです。

例えば、４行３列のゼロ埋めした配列を用意する場合の例です。

• numpy: np.zeros(4,3)
• pdl: zeroes(3,4)

pythonのコードをコピーすると、全てが逆になるので、結構頭が混乱するのですが機械にとっては縦も横もないわけなので慣れてしまえば平気です。

あと、testデータとtrainデータの分割時の:を使った指定の際、From-ToのToの扱いが違うので、うまく合うように調整してやる必要があります。

• numpy 0:N_SIZE (N_SIZEは含まない)
• pdl 0:N_SIZE-1 (N_SIZEは含まないので、マイナス１を加える)

実行結果

$ perl iris.pl

58/100
0.347918856815158
0.224271306476562
0.110087931759781
0.0828852730545221
0.0444535190411007
0.0356200821472588
0.0296333456682152
0.0252544894186067
0.0219212824222726
0.0193079334160957
0.0192846699135605
96/100

ちゃんと学習できています！

実行時間の比較

自分のサーバースペックは以下です。処理時間は、そんなに差はありませんがトータルだと若干pdlの方が２割くらい長いです。ただ、OSの処理時間はperlの方が短め。適当な関数がなかったので少し強引な処理を入れたことも関係しているかもしれませんので、厳密な意味でA2Aの比較ではないです。

• Core™ i7-7700 プロセッサー
• NVIDIA® GeForce® GTX 1070
• DDR4 16GB (8GBx2)

$ time perl iris.pl
（学習結果は省略）
real 0m0.273s
user 0m0.260s
sys 0m0.008s

$ time python iris.py
（学習結果は省略）
real 0m0.227s
user 0m0.180s
sys 0m0.012s

次は、elgoogさんのRubyでディープラーニングを参考に、perlでCNNを実装してみたいと思います。