その他　実装してみた関数について

　ここでは機械学習アルゴリズム自体の話ではなく、そのために実装した関数について書いてみたいと思います。FortLearnerで新たに用意した（と思われる）関数の紹介です。

Sort関数

　ソート関数といえば実用的にはまずクイックソートがあげられるでしょう。さらに実用性のことを思えば、イントロソートを実装すべきでしょう。例えばクイックソートのピボットとしてわざと性能が悪くなるような配列を入力された場合への対応などです。しかしながら、機械学習の入力データは基本的にランダムな並びであることが非常に多いため、対応はしていません。

クイックソート

　FortLearnerはそのように外部からの学習データの投入を考えていないので、単純に配列の真ん中、もしくは要素の最初か最後も加えてその中央値をとるような実装になっています。分割された各配列の要素数が指定された数をした回った場合は挿入ソートに切り替えられます。要素数で分けたのは理由を全く思い出せないのですが、多分そのほうが速かったからでしょう。

　実装の参考もとはこちらです。一部の実装を変更しています。参考もとの引数と定義は（一部変数は削除しています）、

recursive subroutine quicksort(a, first, last)
    implicit none
    real*8  a(*)
    integer first, last
    integer i, j
    ! クイックソートの処理
    if (first < i-1) call quicksort(a, first, i-1)
    if (j+1 < last)  call quicksort(a, j+1, last)
end subroutine quicksort

ですが、FortLearnerでは（こちらも一部変数は削除しています。変数の型については無視してください）、

recursive subroutine quick_sort_r8(vector, num)
    implicit none
    real(kind=8), intent(inout) :: vector(num)
    integer(kind=8), intent(in) :: num
    integer(kind=8)             :: i, j
    ! クイックソートの処理
    if (1 < i-1)    call quick_sort_r8(vector(1:i-1),  i-1)
    if (j+1 < num)  call quick_sort_r8(vector(j+1:num),num-j)
end subroutine quick_sort_r8

のように定義しています。

　このように実装した理由は、参考もとの実装では再帰的に呼ばれるクイックソートに毎回すべての配列を渡しているのは無駄ではないかと思ったらからです。かなり前に実装しましたが、その時は20%ほど速くなった記憶があります。

ブロッククイックソート

　さらにBlockQuickSortも一応実装しています。クイックソートでは配列の左右からスキャンを行い、ピボットとの比較が行われ、適宜交換が行われます。

do
    do while (vector(i) < pivot)
        i=i+1
    end do
    do while (pivot < vector(j))
        j=j-1
    end do
    if (i >= j) exit
    tmp = vector(i);  vector(i) = vector(j);  vector(j) = tmp
    i=i+1
    j=j-1
end do

　この配列の間の空いた$i$と$j$へのアクセスは条件分岐になっているため、BranchMissPredictionが発生しやすいため、これを起こらないようにしてやれば性能が向上すると論文では主張しています。オフセットを用意し、これに格納していきますが、この時If文でなく、vector(i) < pivot（またはpivot < vector(j)）を整数型の変数に格納します。こうすることで分岐予測は発生せず、それぞれの条件に合致した際にはオフセットに格納され、そうでなければ次の要素に移行します。

    if ( num_l .eq. 1 ) then
        start_l = 0
        do i=1, 31
            offset_l(num_l) = i
            one = pivot .lt. vector(l+i-1)
            num_l = num_l + one
        end do
    end if

    if ( num_r .eq. 1 ) then
        start_r = 0
        do i=1, 31
            offset_r(num_r) = i
            one = pivot .ge. vector(r-i+1)
            num_r = num_r + one
        end do
    end if

　比較対象としては C++ std::sortやSuperScalarSampleSortなどいくつかのソート関数を挙げています。発表した直後にstd::sort側に更新が入ったようで、論文タイトルが変更されていました。

BlockQuicksort: Avoiding Branch Mispredictions in Quicksort
→ BlockQuicksort: How Branch Mispredictions don't affect Quicksort

古いバージョンを見るとstd::sortの処理時間を超えてしまったようです。std::sortの変更としてはピボットの選択を少し変更したようです。

The very small difference in the implementation of choosing the second instead of the first element as part of the sample for pivot selection makes a enormous difference when sorting special permutations like decreasingly sorted arrays. This shows how important not only the size of the pivot sample but also the proper selection is. In the other benchmarks both implementations were relatively close, so we do not show both of them.
https://arxiv.org/pdf/1604.06697v1.pdf

その後再度BlockQuickSort側に更新が入ったようで、今度は処理速度がstd::sortより速くなっています。

　実装はしたのですが、FortLearner内ではほぼ使っていません。上記の一度std::sortに負けたとの記述でまあ使わなくていいかなとなってしまっていました。今回改めて調べてみると記憶と違いあれ？と思ったら、今回のようになっていました。std::sortの修正内容についてもう一度確認し、取り込めないかと検討したいと思います。

pbucket_sort

　pはpseudoの意味で、バケットソートもどきとして作成しました。クイックソートの計算量はよく知られているように、入力配列の大きさが$N$としたときに$N \log N$です。pbucket_sortは、バケットソートのように、いくつかのバケツを用意します。例えば、バケツが二つなら閾値は1個、バケツが3つなら閾値は2個のように増えていきます。例えばバケツの数をＢ個用意したとすると、各バケツに対してクイックソートを適用すれば各バケツの計算量は$N/B \log (N/B)$となります。これが、Ｂ個あるわけですから、最終的な計算量は$N \log (N/B)$と計算できます。

　対数の中がちょっと小さくなったｔ程度であまり変わらないように見えますが、例えばN＝100万、B＝1000程度の値を設定すると

\begin{align}
N \log N & \sim 13815510.558 \\
N \log (N/B) &\sim  6907755.27898 \\
\rightarrow \frac{N \log N}{N \log (N/B)} &\sim \frac{13815510.558}{6907755.27898} = 2
\end{align}

計算量オーダーをそのまま計算するのはどうかと思いますが、このように2倍程度高速化できる見込みです。当然、そもそもバケツにフル分ける処理に時間がかかるので、そのままとはいきませんが、ある程度速くなることが確認できました。
　
　ちなみに、閾値はランダムでsqrt(N)個のサンプルを抽出し、そのユニーク値を用いています。バケットも再帰的にやることも考えたのですが、2度以上やると逆に速度が低下したので1度のみにとどめています。

　下記のグラフはソート関数の処理時間の比較です。横軸が要素数、縦軸が要素ごとの処理時間です。FortLearnerで実装されているのは、グラフ下部の判例で言うと、myqsort, bqsort, bqsort w/o sorting_net, pbucket_sortです。最初に3つはクイックソート＋挿入ソート、BlockQuickSort、BlockQuickSortのソーティングネットワークなし（なんでこんなのを実装したかが思い出せません）です。ほかはFortranで利用できるものをいくつかピックアップしています。配列サイズが小さいときはpbucket_sortは遅いですが、ある領域からはずっと早いままとなってくれています。時間計測のコードに関してはこちらを参考にさせていただきました。このように単一の機能の処理時間計測はあまり出てくる数値が信頼できないことが多い気がする（特に要素数が小さい領域で）ので、もっと良い計測方法をご存じの方がいらっしゃいましたらご教授ください。（2年以上前の結果のため、今はかなり違うかもしれません）

　この結果から、配列のサイズが小さい領域では、myqsort(通常のクイックソート＋挿入ソート)を利用するようにしています。BlockQuickSortもこう見るとかなり速いですね。置き換えてやってもよい気がします。

　pbucket_sortは、実はskasortと似たようなアルゴリズムのはずです。思いついた（と思ったときは）ときはうれしかったのですが、あれ？となってちゃんと調べてみると同じような実装でした。

　注意として、私の実装は入力が一様分布の時かつ、ユニークな値ばかりのときに最も性能を発揮します。（ほぼ）ソート済み配列であったり、同じ値が多い場合などは性能が発揮できないはずです（クイックソートを呼ぶようにしています）。

　今はほかにもいろいろなソート関数が出ています。とりあえず今の処理速度で満足してしまったので、今は下記を調べるのはずっと先になると思います。

get_minmax

　例えばExtraTree内では、1次元配列の中の最大値最小値が必要なタイミングがあります。その時単純にmaxval, minvalを適用してしまうと同じ配列を2度も全スキャンしなければなりません。これを避けるために、get_minmaxサブルーチンを定義しています。Fortranでは下記のようなイメージです。

min_val = huge(0d0)
max_val = -huge(0d0)
do i=1, size(vector), 1
    val = vector(i)
    min_val = minval([val, min_val])
    max_val = maxval([val, max_val])
end do

　実際は、LoopUnrollingをしていたり、Cでインラインアセンブラを実装していたりなどをしています。このコンパイル時の分岐が良くなくて、私の現在のPCでしか動かないので、早めに直したいです。

get_matrix_minmax

　さらに上記を拡張したものとして、get_matrix_minmaxを定義しています。これは、2次元配列の各列の最大値最小値の配列を同時に求めるものです。引数を見てみます。
get_matrix_minmax_r8(min_vals, max_vals, mat_t, indices, n_indices, n_rows, n_cols)

最小値の配列。要素数はn_cols
最大値の配列。要素数はn_cols
入力する2次元行列、転置している必要がある。形状は(n_rows, n_cols)
アクセス対象の行番号。要素数はn_indices
４のサイズ。n_indices<=n_rows
入力する2次元行列の行数
入力する2次元行列の列数

　転置している理由は後述します。

get_matrix_count_and_sum_up

　こちらもget_matrix_minmaxと似たような感じです。まず引数を見ます。
get_matrix_count_and_sum_up(sum_vals, cnt_vals, thr_vals, mat_t, y, indices, n_indices, n_rows, n_cols)

列ごとの閾値以下のyの合計。要素数はn_cols
列ごとの閾値以下のサンプルの数。要素数はn_cols
閾値。要素数はn_cols
入力する2次元行列、転置している必要がある。形状は(n_rows, n_cols)
回帰対象の値。2値分類であれば流用可能。
アクセス対象の行番号。要素数はn_indices
４のサイズ。n_indices<=n_rows
入力する2次元行列の行数
入力する2次元行列の列数

入力する2次元行列が転置している理由

　決定木では、現在のノードにあるデータへのアクセスは元データの先頭のポインタと、データのインデックスを用いて行うと説明したと思います。単純にデータを読み込むと、行列の最も内側の添え字は行番号、次が列番号になるように読み込みます。ノードの保持するインデックス番号（＝行番号）が常に連続であることは根ノードでしか保証されません。深くなればなるほどとびとびのメモリアクセスが強制され、キャッシュヒットミスが増加します。

　ここでExtraTreeのコードを眺めてみるとほしいのは

各列の最大値最小値 get_matrix_minmax
閾値以下のデータの数、回帰対象の値の和 get_matrix_count_and_sum_up

のみであることがわかります（閾値は求めた最大値最小値の間から一様乱数で取得します）。つまり、律儀に列ごとにデータへのアクセスをする必要がなく、ある程度は連続であることが保証されている列方向（列番号が増える方向。行方向と列方向の定義が良くわかっていません。使わないようにします。）にアクセスすればキャッシュヒットミスを削減できることが期待できます。実際高速化ができ、列数が増えるほど効果が増加している傾向が見られました。

　当然FortranはColumnMajorなので、"行列の最も内側の添え字は行番号、次が列番号"のように読み込んでいては意味がないので、入力の行列は転置しています。

　この処理はインラインアセンブラで実装してしまったので、このメモリアクセス方法自体がどれだけ効果的かはわかっていません。

read2bin, read_bin

　CSVを読み込みバイナリファイルとして出力する関数と、その出力されたバイナリファイルを読み込む関数です。

　基本的にきれいにしたデータしか読み込むことを想定していません。きれいにしたとは、特徴量計算が済んだ状態です。特徴量計算はPythonでなどで実施しています。文字列が入っていたり、NaNがあったりは全く考えていないため、かなり利用範囲は狭いです。

read2binの引数

　引数を見てみます。
read2bin_2d_32bit(input_file_name, output_file_name, n_samples, n_columns, skip_header, input_dtype, output_dtype)

入力ファイル名
出力ファイル名
行数
列数（read2bin_1d_32bitの場合にはこの引数はない）
ヘッダーを飛ばすか否か
入力データの型
出力データの型

　行数も列数を指定する必要がある部分が一番気に食わない部分です。n_samplesが64bit整数なら保存されるバイナリファイルも中身は64bitの整数か浮動小数点で保存されます。

read_bin_2dの引数

read_bin_2d_r4(file_name, matrix)

読み込み対象のファイル名
読み込み結果の行列

出力されたバイナリファイルの先頭に保存されているデータの形状、型なども入っているので、二つ目の引数と不一致であればエラー文が出ます。

非常に使いにくいですが、毎回CSVからデータを読み込むのはアルゴリズムの動作検証までの時間短縮を思えば一度で済むので、妥協しています。

以上