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ライブラリとしてラップした。
これがテスト計画だ。
- n個のされた範囲(0,\(2^{12}\))の整数値のデータファイルを乱数で生成する。固定値のseedで乱数生成する。
- テストは、まずデータをすべて読み込む
- テストは、RDTSC命令でサイクル数を計測する
- テストは、少なくとも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 おめでとう。
