OpenBLASが速かった

ちょっと行列演算をしたかったので適当にコードを書き散らしていたのですが、速度が必要な演算であるために既成のライブラリを使うこととしました。今回は行列乗算ができれば良いのでBLASを生で使ってもそんなに問題ありません。BLASの中でも高速な実装であるOpenBLASを使い、適当に書いたコードと速度の比較をしてみます。

結論から言うと、OpenBLASの速さを舐めてました。まさか50倍の速度差がでるとは。

比較に用いたマシンは、MacBook Pro 2013 Late 13 inch (2.4 GHz Intel Core i5) です。

インストール

Macの場合はbrewで。LinuxでもOpenBLASのパッケージはたいていのディストリビューションにはあると思います。

$ brew install homebrew/science/openblas

Matrix クラス

1024 × 1024の大きさの実正方行列の掛け算をします。以下N = 1024とします。

const int N = 1024;

適当にMatrixクラスを定義します。精度はそこそこで良いので今回はfloatを使います。

class Matrix {
public:
    Matrix(int row_size, int col_size) :
        row_size_(row_size),
        col_size_(col_size),
        data_(new float[row_size * col_size]) {}

    float& operator()(int row, int col) {　return data_[row * col_size_ + col]; }
    const float& operator()(int row, int col) const { return data_[row * col_size_ + col]; }

    float* data() { return data_.get(); }
    const float* data() const { return data_.get(); }

private:
    int row_size_;
    int col_size_;
    unique_ptr<float[]> data_;
};

operator()を定義していて、行列Aのi行目j列にはA(i, j)でアクセスできるようにしておきます。

1024 × 1024 の行列 A, B を掛け算し、C = ABを計算するものとします。C++の変数としては小文字a, b, cを用いて次のように表すとします。

Matrix a(N, N);
Matrix b(N, N);
Matrix c(N, N);

ナイーブに掛け算をした場合

ナイーブに実装すると次のようになるでしょう。

for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
        float s = 0;
        for (int k = 0; k < N; ++k) {
            s += a(i, k) * b(k, j);
        }
        c(i, j) = s;
    }
}

この場合、b(k, j)の部分がメモリを飛び飛びにアクセスすることになるので速度はでません。clang -O3で試すと、おおよそ6〜7秒程度の時間がかかります。

乗算順序を変更

i, j, kの順序でアクセスするのではなく、i, k, jの順でアクセスするとメモリアクセスが飛び飛びにならないのでいくらか速度が稼げます。 ただしc(i, j)は0で初期化しておく必要があります。

for (int i = 0; i < N; ++i) {
    for (int k = 0; k < N; ++k) {
        for (int j = 0; j < N; ++j) {
            c(i, j) += a(i, k) * b(k, j);
        }
    }
}

これでおよそ1.1〜1.2秒ぐらいです。演算回数としては10243なので1G回の演算があるわけですが、2.4GHzのCPUで1.2秒となると、まあちょっと遅いですねという気持ちになります。

OpenBLAS の場合

関数cblas_sgemmを使います。cblas_sgemmは、行列A, B, Cとスカラーα, βでαAB + βCを計算します。計算結果はCに格納されます。単に掛け算ABを計算したければ、βを0とすれば良いでしょう。

cblas_sgemmの引数はちょっと見ただけではわかりづらいです。

行列の要素はメモリ上では行ごとに並んでいるので、CblasRowMajorを指定する必要があります。FORTRANなどではメモリ上は列ごとにならべるのが普通ですがC/C++では行ごとにならべるのが普通です。また、cblas_sgemmでは転置してから行列を掛け算することもできます。今回は転置などせずに掛け算するので、CblasNoTransを指定します。

cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
    N, N, N,
    1, a.data(), N,
    b.data(), N,
    0, c.data(), N);

これの計算時間はわずか0.03秒でした。

ちなみにEigenでも同じような計算をしてみたのですが、こちらは0.11秒程度でした。それでも速いですね。

2016-01-17
活動日記