コンパイラーを負かす

roguelazer's website: beating the compiler

なかなか面白かったので翻訳して紹介する。

たとえば、97%の場合において、僅かな効率など忘れるべきである。。早すぎる最適化は諸悪の根源である。とはいえ、残りの重要な3%の機会を逃すべからず。

-- Donald Knuth

計測せよ。計測するまで速度の最適化を施してはならぬ。たとえ計測したにせよ、一部のコードが残りを圧倒するまではまだ最適化してはならぬ。

Rob Pike

最新のWebサービスを主体とした技術の業界に長年浸かった我々は、パフォーマンスの問題を忘れがちである。SQLAlchemy ORMの中で行うリクエスト一つが8,9秒かかる中で、関数呼び出しひとつを3ミリ秒最適化したところで何になるというのか。とはいえ、時にはそのような最適化スキルを養っておくのもいいことだ。今回は、ある簡単な課題を最適化して、やっつけ仕事よりどれだけよくできるのかを見ていくことにする。普通はやらないだろうが。

課題

今回の課題はとても簡単なものだ。n個の整数を合計せよ。この課題を公平にするために、整数値は32bit符号付きintに収まるものとし、その合計も32bit符号付きintに収まるものとする。

もちろん、我々は最適化をするのであるから、パフォーマンスの計測には注意を払わねばならぬ。通常なら、我々は単なる時刻を使う（Pythonのtimeitモジュールなどで計測するものだ）。ただし、time.clock()をPythonからただ呼び出すだけで、30マイクロものコストがかかる。おやおや、これはあまりよくないぞ。言い換えてみよう。time.clock()をPythonで呼び出すと、最近のCPUでは60,000サイクル消費するのだ。かわりに、筆者はCPUサイクルをカウントする手法を用いるものとする。これには問題もある（特に、周波数スケーリングが有効になっていたり、複数のコアがあってOSのスケジューラーがリスケジュールを行いたくなった場合など）。だが、これはCPUのパフォーマンスの計測として、最も正確でオーバーヘッドの最小な単位である。Intel x86では、RDTSC命令で実装できる（より具体的に言うと、シリアライズされるRDTSCP命令だ）。筆者はこれを使いやすいようにPythonライブラリとしてラップした。

これがテスト計画だ。

1. n個のされた範囲(0,\(2^{12}\))の整数値のデータファイルを乱数で生成する。固定値のseedで乱数生成する。
2. テストは、まずデータをすべて読み込む
3. テストは、RDTSC命令でサイクル数を計測する
4. テストは、少なくとも100回実行する

この記事のテストは、Intel Core i7-4850HQを搭載した2013年の15インチのrMBPで行われた。また、Intel Core i5-4308Uを搭載した2014年の13インチ rMBPでも実行して、結果がまともであることを確かめた。

最初の試み：Python

Pythonは好きだ。簡単でダックタイピングな言語で、オブジェクト指向と関数型の便利な機能がたくさんあるプログラミング言語だ。ファイルからデータを読み込むのは簡単だ。

import array

data = array.array("i")
with open(sys.argv[1], "r") as f:
    for line in f:
        data.append(int(line.strip()))

最もやっつけな解決法がこれだ。

def sum_naive_python():
    result = 0
    for i in data:
        result += i
    return result

コードはあまりに自明だ。美しい。このコードは、筆者のコンピューターでは、50,000,000 ± 3%サイクルかかる。これは約24ミリ秒だ。だが、筆者は絶対時間についてはそれほど気にしていない。問題なのはスケールするかどうかだ。そこで筆者が使う指標は、1サイクルあたり何個処理できるかだ。理想的な単一命令実行の非スーパースカラーで完璧な分岐予測を持つCPUでは、サイクルあたり1.0個処理できるはずだ。この指標では、最初の試みは悲惨だ。0.01個/サイクルだ。整数の個数辺りのスケールは以下の通り。

とても簡単な改良方法がある。Pythonにはネイティブのsum関数があり、これはCで実装されていて、もっと速いはずだ。次のテスト関数はこうだ。

def sum_native_python():
    return sum(data)

これはマシになった。12,000,000 ± 2%サイクル、0.04個/サイクルだ。まだ悲惨だが、4倍速い。

Python用のイケてる数値ライブラリを使うとどうなるだろうか。numpyはどうだ。このライブラリは最低聞かされたFORTRANコードを呼んでくれる。すげー速いはずだ。

ダメだ。NumPyのオーバーヘッドで相殺されるようだ。0.004個/サイクルしかでない（ネイティブsum関数を組み込みのarray型に適用する場合の10分の1だ）

より高みへ：C

Pythonでそれ以上速くできないのであれば、Cに変えろというのは、Pythonでよくあるパターンだ。Cで同様のテスト環境を用意した。ファイルからヒープに確保されたintの配列に読み込み（少々長いのでわざわざコードはみせない）。この数値群を足し合わせるC関数の最初の実験がこれだ。

int sum_simple(int* vec, size_t vecsize)
{
    int res = 0;
    for (int i = 0; i < vecsize; ++i) {
        res += vec[i];
    }
    return res;
}

最初のPythonコードと同じくらい綺麗だ。さて結果は？

3,423,072サイクル(~1.5ミリ秒)、0.162個/サイクルだ。これはPythonより10倍速い。

注意：Apple LLVM 6.0 / clang-600.0.56 -O0でコンパイルした。

もっと速くできるだろうか。

そういえば昔、ループの手動展開は高パフォーマンスであると聞いたことがある。ループを展開すると、ループの最後のジャンプのオーバーヘッドを隠して、分岐予測の為事を代わりに引き受けたことになる。4回展開版が以下だ。

int sum_unroll_4(int* vec, size_t vecsize)
{
    int res = 0;
    int i;
    for (i = 0; i < vecsize - 3; i+=4) {
        res += vec[i];
        res += vec[i+1];
        res += vec[i+2];
        res += vec[i+3];
    }
    for (; i < vecsize; ++i) {
        res += vec[i];
    }
    return res;
}

展開されたループのパフォーマンスはどうか。ほぼ同じだ（0.190個/サイクル）

ひょっとしたら展開数が足りないのではないか。では、8回ループ展開をしてみよう。

int sum_unroll_8(int* vec, size_t vecsize)
{
    int res = 0;
    int i;
    for (i = 0; i < vecsize; i+=8) {
        res += vec[i];
        res += vec[i+1];
        res += vec[i+2];
        res += vec[i+3];
        res += vec[i+4];
        res += vec[i+5];
        res += vec[i+6];
        res += vec[i+7];
    }
    for (; i < vecsize; ++i) {
        res += vec[i];
    }
    return res;
}

やれやれ、0.192 PPC(Points Per Cycles)だ

どうやら、最近のプロセッサーの分岐予測というのはめちゃめちゃ優秀らしい。ループ展開は、とても短いループか、あるいは何らかの理由で分岐予測が最悪の結果になるループには有効だが、まず時間の無駄というものだ。

最適化コンパイラー

Cを書いた経験がある読者は、ここで疑問に思っていることだろう。なぜ筆者は-O0でコンパイルしているのだ？、と。最新のコンパイラーは何十年もの最適化手法が積み上げられており、筆者ごときが書いた程度のものは余裕でこすことができるだろう、と。

実にそのとおりだ。-O3を使うと、精製されたコードは他の方法を消し炭すら残らないほどの煙にしてふっ飛ばしてしまう。4.462個/サイクルだ。

待てよ・・・1個/サイクル以上だと？　どうやったらそんなことが可能なのだ。

さて、読者よ。プロセッサーの歴史のお勉強の時間だ。

複数命令発行、スーパースカラー、アウトオブオーダー、ベクトル計算機

今は昔、コンピュータープロセッサーというものは、とても簡単だった。プロセッサーは命令を一つづつ読み込み、1ないしは数サイクルを使って処理し、結果を吐き出すものだった。

これは、1965年にCray 6600がリリースされる以前の話で、デスクトップでも1995年にP6 Pentium Proが発売されたことで過去の話になった。この2つのプロセッサーは、スーパースカラーアーキテクチャの代表例だ。

簡単に言うと、スーパースカラーアーキテクチャーとは、複数の実行ユニット（大抵目的別に特化されている。たとえば2つの汎用数値演算ユニットといくつかのアドレス計算ユニットとひとつのSIMDユニットとひとつのFPUユニットなど）があり、複数の命令を並列に実行するために別々の実行ユニットに振り分けることができる。

筆者の持っているIntel Haswell CPUのアーキテクチャーは、以下のようになっている。

つまり、複数の命令を同時に実行できるということで、1.0加算/サイクルというのは、遅すぎるということだ。

最後の要素は、謎めいた略語であるSIMD(Single Instruction, Multiple Data)だ。SIMD命令はベクトル命令とも呼ばれていて、演算を複数のデータに対して一度に適用できる。例えば、IntelのSSEユニットでは、演算を4つの32bit整数値か、2つの64bit浮動小数点数に対して適用できる。

x86には多数のSIMDユニットがある。MMX拡張から始まり、SSE, SSE2, SSE3, SSE4.1, AVX, AVX2、そしてまだリリースされていないAVX-512。ここではその技術的詳細については述べないが、詳しく知りたければ、ここにすべて載っている。

Haswellプロセッサーでは、2つのベクトルInt ALUがある、理論上、2つのAVX2処理を並列に行える。それぞれが8個の32bit整数値を同時に処理できる。Intel 64 and IA32 Architectures Optimization Reference Manualによれば、それぞれの命令は1サイクルのレイテンシーがあり、つまり理論上のPPCは16.0になる。だいぶ高い。実際、これは不可能なほどに高い。1サイクルごとに16個32bit整数値をCPUに転送するには、64バイト/サイクルのメモリースループットが必要であり、これは、筆者のCPUに当てはめれば、147GB/sのメモリースループットになる。筆者の環境のPC3-12800 RAMは、最大一本あたり12.8GB/s（トータルで25.6GB/s）であり、これはたったの11バイト/サイクルだ。つまり、メモリの最大のスループットを考えれば、2.75PPCをすこし上回るほどだ。これ以上高いスコアは、筆者のCPU上の128MBものeDRAM L4キャッシュに残っていたものである。

コンパイラーを負かす

さて、コンパイラーは負かせるか？　4.462 PPCより上を叩き出せるか。もちろん無理だが、データが広帯域のL3やL4キャッシュに乗るほど小さければ、まあ、可能ではある。いかがその例だ。PPCは4.808だ。

int sum_avx2_unroll4(int* vec, size_t vecsize)
{
    unsigned int vec_out[8] __attribute__ ((aligned(64)));
    size_t processed = 0;
    asm(
        "XORQ %%rax,%%rax\n\t"
        "VPXOR %%ymm0,%%ymm0,%%ymm0\n\t"
        "VPXOR %%ymm5,%%ymm5,%%ymm5\n\t"
        "VPXOR %%ymm6,%%ymm6,%%ymm6\n\t"
        "VPXOR %%ymm7,%%ymm7,%%ymm7\n\t"
        "VPXOR %%ymm8,%%ymm8,%%ymm8\n\t"
        "top%=:\n\t"
        "VMOVDQA (%1, %%rax, 4),%%ymm1\n\t"
        "VMOVDQA 32(%1, %%rax, 4),%%ymm2\n\t"
        "VMOVDQA 64(%1, %%rax, 4),%%ymm3\n\t"
        "VMOVDQA 96(%1, %%rax, 4),%%ymm4\n\t"
        "ADDQ $32,%%rax\n\t"
        "VPADDD %%ymm1,%%ymm5,%%ymm5\n\t"
        "VPADDD %%ymm2,%%ymm6,%%ymm6\n\t"
        "VPADDD %%ymm3,%%ymm7,%%ymm7\n\t"
        "VPADDD %%ymm4,%%ymm8,%%ymm8\n\t"
        "CMP %2,%%rax\n\t"
        "JL top%=\n\t"
        "VPADDD %%ymm5,%%ymm0,%%ymm0\n\t"
        "VPADDD %%ymm6,%%ymm0,%%ymm0\n\t"
        "VPADDD %%ymm7,%%ymm0,%%ymm0\n\t"
        "VPADDD %%ymm8,%%ymm0,%%ymm0\n\t"
        "VMOVDQA %%ymm0, 0(%3)\n\t"
        "MOVQ %%rax, %0\n\t"
        : "=r"(processed)
        : "r"(vec), "r"(vecsize), "r"(vec_out)
        : "rax", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "cc"
    );
    int res = 0;
    for (int i = 0; i < 8 ; ++i) {
        res += vec_out[i];
    }
    for (; processed < vecsize; ++processed) {
        res += vec[processed];
    }
    return res;
}

これは4回展開したAVX2ベクトル化版だ。また、独立した複数の加算レジスターに合計させることにより、IPCを僅かに稼いでいる。

本当にコンパイラーを打ち負かしたのだろうか。PPCの違いを比べてみる。

技術的に言えば、筆者のコードは負けることはない。参考に、これがコンパイラーの生成するアセンブリだ。

## BB#2:                                ## %vector.body.preheader
    pxor    %xmm0, %xmm0
    xorl    %ecx, %ecx
    pxor    %xmm1, %xmm1
    .align  4, 0x90
LBB8_3:                                 ## %vector.body
                                        ## =>This Inner Loop Header: Depth=1
    movdqa  %xmm1, %xmm2
    movdqa  %xmm0, %xmm3
    movdqu  (%rdi,%rcx,4), %xmm0
    movdqu  16(%rdi,%rcx,4), %xmm1
    paddd   %xmm3, %xmm0
    paddd   %xmm2, %xmm1
    addq    $8, %rcx
    cmpq    %rcx, %rax
    jne LBB8_3
## BB#4:
    movq    %rax, %rcx
LBB8_5:                                 ## %middle.block
    paddd   %xmm1, %xmm0
    movdqa  %xmm0, %xmm1
    movhlps %xmm1, %xmm1            ## xmm1 = xmm1[1,1]
    paddd   %xmm0, %xmm1
    phaddd  %xmm1, %xmm1
    movd    %xmm1, %eax
    cmpq    %rsi, %rcx
    je  LBB8_8
## BB#6:                                ## %.lr.ph.preheader13
    leaq    (%rdi,%rcx,4), %rdx
    subq    %rcx, %rsi
    .align  4, 0x90
LBB8_7:                                 ## %.lr.ph
                                        ## =>This Inner Loop Header: Depth=1
    addl    (%rdx), %eax
    addq    $4, %rdx
    decq    %rsi
    jne LBB8_7

コンパイラーは2回展開のSSE2レジスターを使う選択をしている。いえい。

結論

Pythonや、その他の動的型付け言語は、十分に速いかもしれないが、速いことは常にいいことだ。Cのような低級言語で速度を稼げるかもしれないし、アセンブリを手書きすることでもう少し上を行けるかもしれないが、コンパイラーを信頼したほうがいい。しかし、もし大量の12bitの数値を限りなく高速に足し合わせたいのであれば、まあ、これが答えだ。

Algorithm	Cycles for 500k Points	Points per Cycle
sum_numpy	110681288 ± 0% Cycles	0.005 ± 0% PPC
sum_naive_python	48255863 ± 1% Cycles	0.010 ± 1% PPC
sum_native_python	12063768 ± 2% Cycles	0.041 ± 2% PPC
-O0 simple	3095247 ± 1% Cycles	0.162 ± 1% PPC
-O2 simple	2889994 ± 1% Cycles	0.173 ± 1% PPC
-O0 unroll_4	2630305 ± 2% Cycles	0.190 ± 2% PPC
-O0 unroll_8	2597596 ± 0% Cycles	0.192 ± 0% PPC
-O2 unroll_8	2461759 ± 1% Cycles	0.203 ± 1% PPC
-O2 unroll_4	2419625 ± 1% Cycles	0.207 ± 1% PPC
-O2 unroll_4_asscom	1134117 ± 1% Cycles	0.441 ± 1% PPC
-O0 unroll_4_asscom	1091086 ± 2% Cycles	0.458 ± 2% PPC
-O3 unroll_4	393973 ± 4% Cycles	1.269 ± 4% PPC
-O3 unroll_8	334532 ± 0% Cycles	1.495 ± 0% PPC
-O0 asm_unroll_8	266224 ± 3% Cycles	1.878 ± 3% PPC
-O3 unroll_4_asscom	254489 ± 1% Cycles	1.965 ± 1% PPC
-O2 asm_unroll_8	227831 ± 1% Cycles	2.195 ± 1% PPC
-O3 asm_unroll_8	227390 ± 0% Cycles	2.199 ± 0% PPC
-O0 avx2_and_scalar_unroll	149051 ± 0% Cycles	3.355 ± 0% PPC
-O3 avx2_and_scalar_unroll	146261 ± 0% Cycles	3.419 ± 0% PPC
-O2 avx2_and_scalar_unroll	144545 ± 0% Cycles	3.459 ± 0% PPC
-O3 sse2	125397 ± 2% Cycles	3.987 ± 2% PPC
-O0 sse2	123398 ± 0% Cycles	4.052 ± 0% PPC
-O2 sse2	119917 ± 0% Cycles	4.170 ± 0% PPC
-O3 simple	112045 ± 0% Cycles	4.462 ± 0% PPC
-O0 avx2	111489 ± 0% Cycles	4.485 ± 0% PPC
-O2 avx2	110661 ± 0% Cycles	4.518 ± 0% PPC
-O3 avx2_unroll4_mem_address	110342 ± 0% Cycles	4.531 ± 0% PPC
-O3 avx2_unroll4	109213 ± 1% Cycles	4.578 ± 1% PPC
-O3 avx2	107011 ± 0% Cycles	4.672 ± 0% PPC
-O2 avx2_unroll4	106253 ± 0% Cycles	4.706 ± 0% PPC
-O2 avx2_unroll4_mem_address	105458 ± 1% Cycles	4.741 ± 1% PPC
-O0 avx2_unroll4	103996 ± 0% Cycles	4.808 ± 0% PPC
-O0 avx2_unroll4_mem_address	101976 ± 0% Cycles	4.903 ± 0% PPC

おめでとう。