Eigen vs ArrayFire LU部門 - パラレルに恋して

この記事はGPGPU Advent Calendar19日目の記事です

山田「まじでやばいです」
w_o大先生「あと四時間あるよ」

という会話をしたのが20時のこと。

ArrayFire vs Boltをやろうとしたのですが、BoltがLinuxのSDKに入っていない疑惑が漂ってちょっと探してみたのですがあんまり見つからなかったので、今日はEigen vs ArrayFireです
はじめに謝っておきます。
このGPGPUAdCの2日目に偉大なるw_o大先生の論文紹介にありましたが……GPUとCPUの性能比較はCPUでガチガチにしたあとで行われるべき、という話があります。
Eigenの内部を見まくったわけではないのですが、ガチガチに最適化されているかというと、それがわかりませんというのが率直な話です。とはいえ、ある程度最適化はされているようです。ですので、中身を良く知らないままArrayFireと勝負させるのはあまりフェアではない、というご指摘はもっともです。

とはいえ、僕の12月のテーマはAccelerEyesがんばれ、ですのでその辺はご容赦を。
……だって、Jacketが買収されたんだもん、可哀想じゃん

あと、こういったライブラリを取り上げる理由として「GPUもCPUも、簡単に書けて高性能」というのが理想である、という僕の持論があります
GPUだってコレぐらい楽に書ければもっと使い道があるかもしれない、誰かが導入の検討をするかもしれない。もしかしたら誰かの目に留まるかもしれない
これぐらい速ければ使い物になるなぁ、とか誰かが思うかもしれない
というもくろみがあります
なにもそれはGPUに限ったことではないので、比較対象としてCPUでも簡単に書けるライブラリが必要だろう。ということでEigenさんにご登場いただきました

ところでなんでLUかというと、
「速いLU分解のライブラリない？」
という話をもってこられたので、それのアンサーです。深い意味はありません

Eigenとは？

Eigen
超高速、超ポータブルな行列演算ライブラリです
山田がとやかくいう以上に皆さんご存知だと思うんでこの話はあんまりしない方がいい。多分。墓穴を掘る。

ArrayFireとは？

第一回読んでください

EigenとArrayFireの書き方について

EigenとArrayFireは非常に似通った書き方が可能なライブラリです。すばらしい。
こっちがEigen

#include <eigen3/Eigen/Core>
#include <iostream>
#include "gettimeofday.h"

#define PRINT_MAT(X) cout << #X << ":\n" << X << endl << endl

using namespace std;
using namespace Eigen;

int main()
{   
    const size_t num = 1024;

    MatrixXd A = MatrixXd::Random(num, num);
    MatrixXd B = MatrixXd::Random(num, num);

    volatile double t1 = gettimeofday_sec();

    MatrixXd C = A * B;

    volatile double t2 = gettimeofday_sec();

    std::cout << "time: " << t2 - t1 << " sec." << std::endl;

#ifdef ENABLE_PRINT
    PRINT_MAT(A);
    PRINT_MAT(B);
    PRINT_MAT(C);
#endif

    return 0;
}

こっちがArrayFire

#include <iostream>
#include <cstdio>
#include <arrayfire.h>
#include "gettimeofday.h"

int main()
{   
    const size_t num = 1024;

    af::array A = af::randu(num, num);
    af::array B = af::randu(num, num);

    volatile double t1 = gettimeofday_sec();

    af::array X = af::solve(A, B);

    volatile double t2 = gettimeofday_sec();

//      af::print(A);
//      af::print(B);

    std::cout << "time: " << t2 - t1 << " sec." << std::endl;

    return 0;
}

さて、最適化ありありでその両者のパフォーマンスはといえば

Eigen	time: 0.157702 sec.
ArrayFire	time: 0.00214005 sec.

さすがはGPUだぜ…行列だけは速いな！

さて、ここまでおさらいです

Linear algebra

こう書くとかっこいい（こなみかん
ArrayFireとEigenは、線形代数のどういったアルゴリズムをサポートしているのかについてを以下に述べます

	ArrayFire	Eigen
LU	○	○
QR	○	○
cholesky	○	○
hessenberg	○	○
eigen	○	○
SVD	○	○

つまり、Eigenでできることは、ArrayFireでもできます。ただし、Eigenがサポートしている全てのアルゴリズムがArrayFireで使えるというわけではありません。まぁむべなるかな。

LUのかきかた

LUを、ArrayFireとEigenではどうやって書くのでしょうか？
またしても一行です。すばらしい

Eigen

#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/LU>
#include <iostream>
#include "gettimeofday.h"
#include <cstdlib>

#define PRINT_MAT(X) cout << #X << ":\n" << X << endl << endl

using namespace std;
using namespace Eigen;

int main(int argc, char** argv)
{
    size_t num = 6;
    if(argc >= 2) {
        num = atoi(argv[1]);
    }

    MatrixXf A = MatrixXf::Zero(num, num);
    for(int j = 0; j < num; ++j) {
        for(int i = 0; i < num; ++i) {
            A(j,i) = static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
        }
    }

    volatile double t1 = gettimeofday_sec();

    FullPivLU< MatrixXf > flu(A);

    volatile double t2 = gettimeofday_sec();

    PartialPivLU< MatrixXf > plu(A);

    volatile double t3 = gettimeofday_sec();

    std::cout << "FullPiv time: " << t2 - t1 << " sec." << std::endl;
    std::cout << "PartPiv time: " << t3 - t2 << " sec." << std::endl;

#ifdef ENABLE_PRINT
    PRINT_MAT(A);
#endif

    return 0;
}

ArrayFire

#include <iostream>
#include <cstdio>
#include <arrayfire.h>
#include <vector>
#include <cstdlib>
#include "gettimeofday.h"

int main(int argc, char** argv)
{

    size_t num = 6;
    if(argc >= 2) {
        num = atoi(argv[1]);
    }

    std::vector<float> vec(num * num);
    for(int i = 0; i < num; ++i) {
        vec[i] = static_cast<float>(rand()) / static_cast<float>(RAND_MAX);
    }

    af::array A(num, num, &vec[0]);

    volatile double t1 = gettimeofday_sec();

    af::array lu = af::lu(A);

    volatile double t2 = gettimeofday_sec();

#ifdef ENABLE_PRINT
    af::print(A);
//      af::print(lu);                                                                                                                                          
#endif

    std::cout << "time: " << t2 - t1 << " sec." << std::endl;

    return 0;
}

にたような感じになりましたね。いやはや。初期化処理はほぼ同等のことをやっていますが、Eigenでポインタを渡して初期化ができるのかどうか怪しかったので保留
Eigenの場合、単純にLUといってもアルゴリズムによって違うんですが。部分ピボットかそうでないかがFullとPartialですね
Eigenの解説じゃないのでここではあまり言及しません

さて、双方を1024 * 1024要素でLU分解かけてみますと、結果がどうなるかというと

Eigen Full	time: 0.3918 sec.
Eigen Partial	time: 0.0314651 sec.
ArrayFire	0.0242221 sec.

部分ピボットHAEEEEEEEE!!!!
しかしそれでもArrayFireの高速っぷりはゆらがない。
ちなみに言うと、ArrayFireはサイズを大きくしていっても時間の変化がほとんどないのがすてきなところです。大規模行列にも適用可能！
上限？あんまりまわすとGPU燃えるんで勘弁してください

まとめ

Eigenは十分に高速ですが、それでもArrayFireは速いです。
AccelerEyesさんには割と本気でがんばって頂きたいので、ArrayFireを導入して誰かいいアプリケーション作ればいいと思います(こなみかん

参考:
http://www.singularpoint.org/blog/c/eigen_1/