不確定特異点

広く深く、ところどころ超深く

Python (NumPy) と Common Lisp (LLA) で行列積の実行速度を比較する

近年、機械学習や Deep Learning などのデータサイエンス分野を筆頭に数値計算の需要が高まってますが、その中でよく使われているのが Python の NumPy というライブラリです。NumPy を使うことで、動的で柔軟なスクリプト言語上で比較的高速に数値計算が可能になり、NumPy を使った SciPy や matplotlib などのライブラリも充実していることから、近頃の Python 人気が高まっているのだと思います。

個人的には Common Lisp が好きなので、NumPy との速度差が気になるところです。そこで単精度浮動小数点数を要素とする100行100列の行列積を対象に、PythonCommon Lisp、C でベンチマークをとって比較してみました。(あ、ちなみに僕は Python も好きですよ。)

実行環境は以下の通りです。

$ uname -a
Linux vagrant-ubuntu-trusty-64 3.13.0-62-generic #102-Ubuntu SMP Tue Aug 11 14:29:36 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux

Python (NumPy)

まずは Python (NumPy) です。バージョンは 2.7.6 です。

# coding: utf-8
import numpy as np

N = 100

def simple_gemm(ma, mb):
    """NumPyの行列計算を使わずに実装した版。とてつもなく遅い。"""
    mc = np.zeros((ma.shape[0], mb.shape[0]))
    for i in range(ma.shape[0]):
        for j in range(mb.shape[1]):
            for k in range(ma.shape[1]):
                mc[i][j] += ma[i][k] * mb[k][j]
    return mc

def numpy_gemm(ma, mb):
    """NumPyによる高速な行列計算(BLASライブラリを利用)"""
    return np.dot(ma, mb)

# N行N列の行列を生成して、要素を[0, 1.0]の範囲の乱数で初期化
ma = np.random.rand(N, N)
mb = np.random.rand(N, N)

def run_simple(n = 10000):
    for i in range(n):
        simple_gemm(ma, mb)

def run_gemm(n = 10000):
    for i in range(n):
        numpy_gemm(ma, mb)

simple_gemm は比較用に NumPy を使わずに実装したもので、numpy_gemm が NumPy の機能を使ったものです。ipython 上で実行時間を計測してみます。

In [1]: run test.py

In [2]: time run_simple (10)
CPU times: user 7.5 s, sys: 8.29 ms, total: 7.51 s
Wall time: 7.68 s

In [3]: time run_gemm (10000)
CPU times: user 10.3 s, sys: 12.3 ms, total: 10.3 s
Wall time: 10.6 s

わかりづらいですが、simple_gemm を 10 回、numpy_gemm を 10000 回実行してます。simple_gemm の方は遅すぎて 10 回しか計算しませんでした。一方、 numpy_gemm の方は 70 倍くらい高速化されてます。この高速化の秘密を追って NumPy のソースコードを調べてみると、BLAS ライブラリに行き着きました。

BLAS (Basic Linear Algebra Subprograms)

BLAS とはもともと Fortran で書かれた線形代数演算ルーチンライブラリのデファクトスタンダードであり、現在はオリジナルのリファレンス BLAS の他に Intel MKL, ATLAS, OpenBLAS(GotoBLAS) など色々あるようです。どのライブラリがインストールされていてリンクされるかは次のコマンドで確認できます。

$ update-alternatives --display libblas.so.3
libblas.so.3 - manual mode
  link currently points to /usr/lib/libblas/libblas.so.3
/usr/lib/libblas/libblas.so.3 - priority 10
  slave libblas.so.3gf: /usr/lib/libblas/libblas.so.3
/usr/lib/openblas-base/libblas.so.3 - priority 40
  slave libblas.so.3gf: /usr/lib/openblas-base/libblas.so.3
Current 'best' version is '/usr/lib/openblas-base/libblas.so.3'.

使うライブラリを変更するには、--display の代わりに --config を指定します。

$ update-alternatives --config libblas.so.3
There are 2 choices for the alternative libblas.so.3 (providing /usr/lib/libblas.so.3).

  Selection    Path                                 Priority   Status
------------------------------------------------------------
  0            /usr/lib/openblas-base/libblas.so.3   40        auto mode
* 1            /usr/lib/libblas/libblas.so.3         10        manual mode
  2            /usr/lib/openblas-base/libblas.so.3   40        manual mode

Press enter to keep the current choice[*], or type selection number:

この OS には OpenBLAS とリファレンス BLAS がインストールされていて、リファレンス BLAS が使われていることがわかりました。

BLASループアンローリングやキャッシュヒット率の向上など等 CPU の性能を最大限に発揮できるような最適化が施されているため、これに Common Lisp Only で立ち向かうのは大変そうです。というわけで Common Lisp でも BLAS を使えるライブラリがないか調べてみました。へたれですみません。

Common Lisp (LLA)

LLA というライブラリがありました。

github.com

BLASLAPACK といった有名な数値計算ライブラリを Common Lisp から使えるようにするためのラッパーライブラリであり、NumPy のようなリッチな機能は搭載されてませんが、シンプルな分、今回の目的には十分です。

LLA が使うライブラリはデフォルトで libblas.so が指定されてますが、別のパスを設定する場合は README にあるように、*lla-configuration* という CL-USER パッケージのシンボルにライブラリへのパスを .sbclrc に設定しておきます。

(defvar *lla-configuration*
        '(:libraries ("/usr/lib/libblas/libblas.so.3.0")))

通常は libblas.so から適当なライブラリをシンボリックリンクしているため、上記の設定は不要だと思いますが、システム非標準のパスにあるライブラリを指定する場合などに使うことになるのかもしれません。

以下のようなベンチマーク用プログラムを作りました。Common Lisp の処理系は SBCL 1.2.15 です。

(defparameter *N* 100)

(defun make-matrix (rows cols)
  (make-array (list rows cols) :initial-element 0.0 :element-type 'single-float))

(defun randomize (m)
  "行列の要素を[0, 1]の範囲の乱数で初期化する"
  (destructuring-bind (rows cols)
      (array-dimensions m)
    (dotimes (row rows)
      (dotimes (col cols)
        (setf (aref m row col) (random 1.0)))))
  m)

(defun simple-gemm (ma mb)
  "Common Lispのみを使った行列積の計算"
  (declare (optimize (speed 3) (debug 0) (safety 0)))
  (declare (type (simple-array single-float (* *)) ma mb))
  (let ((rows (array-dimension ma 0))
        (cols (array-dimension mb 1)))
    (declare (type fixnum rows cols))
    (let ((result (make-matrix rows cols)))
      (declare (type (simple-array single-float (* *)) result))
      (dotimes (row rows)
        (dotimes (col cols)
          (dotimes (k cols)
            (incf (aref result row col)
                  (* (aref ma row k) (aref mb k col))))))
      result)))

(defun run (gemm &optional (n 10000))
  "gemmをn回計算して実行時間を測定する。行列の初期化は計測しない。"
  (let ((ma (randomize (make-matrix *N* *N*)))
        (mb (randomize (make-matrix *N* *N*))))
    (time
     (loop repeat n
        do (funcall gemm ma mb)
        finally (return nil)))))

(defun run-simple-gemm (&optional (n 10000))
  "simple-gemmのベンチマーク"
  (run #'simple-gemm n))

(defun run-lla-mm (&optional (n 10000))
  "llaのmm関数のベンチマーク"
  (run #'lla:mm n))

Python版と同様に、BLAS を使わない simple-gemmBLAS を使った lla-mm との比較をします。単精度浮動小数点数を要素とする100行100列の行列積を 10000 回実行した場合の実行時間です。

CL-USER> (run-simple-gemm 10000)
Evaluation took:
  33.395 seconds of real time
  33.236077 seconds of total run time (33.148582 user, 0.087495 system)
  [ Run times consist of 0.032 seconds GC time, and 33.205 seconds non-GC time. ]
  99.52% CPU
  80,349,715,978 processor cycles
  401,280,592 bytes consed
  
NIL
CL-USER> (run-lla-mm 10000)
Evaluation took:
  8.212 seconds of real time
  8.167335 seconds of total run time (8.137797 user, 0.029538 system)
  [ Run times consist of 0.027 seconds GC time, and 8.141 seconds non-GC time. ]
  99.45% CPU
  19,756,699,566 processor cycles
  406,440,112 bytes consed
  
NIL

Common Lisp の simple 版は optimize により speed と型に最適化を施してます。これにより、Python では時間がかかりすぎていた simple 版も最適化した Common Lisp ならぼーっとしている間に計算できます。一方、BLAS を使った LLA 版は simple 版に比べて 4 倍程度高速化され、また NumPy のそれよりも 20% 弱高速化されていることがわかります。PythonCommon Lisp も同じライブラリを使っているので、この差はライブラリ呼び出しに伴うオーバーヘッドの差と考えられます。

では、そもそもライブラリを直接実行した場合の時間はどのくらいなのか?というのが気になります。というわけで、真打ちの C 言語版です。

C

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>

#define N 100

/* BLASを使わずに実装した行列計算 */
void simple_gemm(const float *ma, const float *mb, float *mc)
{
    int i, j, k;
    float sum;
    
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            sum = 0;
            for (k = 0; k < N; k++) {
                sum += ma[i*N+k]*mb[k*N+j];
            }
            mc[i*N+j] = sum;
        }
    }
}

/* BLASを使った高速な行列計算 */
void cblas_gemm(const float *ma, const float *mb, float *mc)
{
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                N, N, N,
                1.0, ma, N, mb, N,
                0.0, mc, N);
}

/* 行列の要素を[0, 1]の範囲の乱数で初期化 */
void randomize(float *m)
{
    int i;
    
    for (i = 0; i < N; i++)
        m[i] = rand() / (float)RAND_MAX;
}

int main(int argc, char *argv[])
{
    float ma[N*N];
    float mb[N*N];
    float *mc;

    /* arg[2] で計算回数を指定 */
    long times = strtol(argv[2], NULL, 10);
    
    randomize(ma);
    randomize(mb);

    /* arg[1] で実行時間を測定する行列計算の関数を指定 */
    if (strcmp(argv[1], "simple") == 0) {
        while (times-- > 0) {
            mc = (float*)calloc(sizeof(float), N*N);
            simple_gemm(ma, mb, mc);
            free(mc);
        }
    }
    else if (strcmp(argv[1], "cblas") == 0) {
        while (times-- > 0) {
            mc = (float*)calloc(sizeof(float), N*N);
            cblas_gemm(ma, mb, mc);
            free(mc);
        }
    }
    else {
        printf("unknown function\n");
    }
    
    return 0;
}

これまで同様、BLAS を使わない simple_gemmBLAS を使った cblas_gemm との比較です。

$ gcc -Wall -O2 test.c -lblas && time ./a.out simple 10000

real    0m9.999s
user    0m9.770s
sys     0m0.021s
$ gcc -Wall -O2 test.c -lblas && time ./a.out cblas 10000

real    0m0.417s
user    0m0.408s
sys     0m0.000s

simple 版に対して BLAS 版は約 24 倍高速化されました。これは Python BLAS 版 の 25 倍、Common Lisp BLAS 版の 20 倍高速です。記憶域の malloc/free をベンチマークに含めるかどうか悩みましたが、PythonCommon Lisp も計算結果を返却するためのオブジェクトを返す API になっているため、C も malloc/free を入れました。出力結果 mc をヒープではなくスタックから取るようにすれば、malloc/free のオーバーヘッドがなくなるため、さらに倍近く高速化されましたが、現実的な使い方ではないため比較する意味は薄いと思います。

まとめ

上記の実行結果をまとめます。なお、今回は簡易的な実験であり、入力する行列サイズが固定だったり、ライブラリによる差も比較できていないため、ベンチマークとしては不十分です。例えば行列サイズを変えて試すと、スケールする部分とスケールしない部分が見えてくるため、言語間での速度比は変わると思います。

100 行 100 列の行列積を 10000 回実行した結果

言語 ライブラリ SIMPLE版実行時間 BLAS版実行時間 備考
Python NumPy 7680[s] 10.6[s] SIMPLE版は推定値
Common Lisp LLA 33.395[s] 8.212[s] SIMPLE版はoptimize適用後
C(heap) - 9.999[s] 0.417[s] malloc/freeあり
C(stack) - 9.512[s] 0.243[s] malloc/freeなし

動的なプログラミング言語の特徴は、アルゴリズム検討などの実験的なコードを書く場合にも向いていると思いますが、このようなベンチマークをとってみると、数値計算を多用するデータサイエンス分野で Python が注目されているのと同じように、それよりも高速なコードを書ける Common Lisp も十分活躍できる気がしてきます。