こんにちは。労働者です。とあるプログラムで学生さんの課題を添削していたら面白い話に出会いました。
僕は今、主に学部生向けのインターン研修的なプログラムでメンターなるものをやっています。メンターとしての仕事は、学生さんの課題へフィードバックを返し、Office Hourというセッションを毎週設けて質問受けやCSに関するトークを行うといった内容になっています。今回話題に取り上げるのはその中の課題の1つ、「行列積のプログラムを書いて時間を計測せよ」という何気ない話で、続く課題たちのいわば前座のようなものです。こういったところに沼は隠されているものですね。
担当している学生さんたちが細かい実験を行ってくれて以下のような疑問が提示されました。
「行列積の計算が N = 1024のときだけ N = 1023, 1025のときに比べて3倍遅いのはなぜ？」
配列のサイズが2のべき乗になるのは避けるべきといった話は聞いたことがあるような(ないような)気がしたのですが、自分にも答えが出ませんでした。というわけで調査を行いましたので、その結果を書きます。
※正確には「行列積の計算が N = 2^nのときだけ N = 2^n+1のときに比べてかなり遅いのはなぜ？」といった結果を提出してくれました。
結論
多分キャッシュライン衝突が原因です。
行列を縦方向にアクセスしていくステップで、Nが2べき乗の際に同じキャッシュタグを踏みまくり、キャッシュの中に収まりきらなくなってキャッシュロードの時間がかかる、ということだと思います。以下で詳細を語ります。
今回は行列積を例にして考えましたが、一番の重要ポイントは「縦方向のアクセスをする機会があるか」というところであり、行列積以外の問題にも適用できる話だと思います。みなさん、2のべき乗サイズの多次元配列は避けましょう。
前提知識とスコープ
ここで考える行列はint型の要素を持つN=1024付近の行列で、mallocなどで連続するメモリを確保したとします。行列積の計算は特に工夫せず、以下のようなシンプルな実装を仮定します。
for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) for (int k = 0; k < N; k++) Answer[i][j] += Matrix1[i][k] * Matrix2[k][j]
メモリ上で行列の要素がどう保管されているかというと、１行目のデータが左から順に入り、次に２行目、３行目…と並んで入っています。
ところでCPUには一度データを探しにメモリへアクセスすると、そのメモリから連続する部分をキャッシュにロードする機能があります。 キャッシュにロードされていると、2度目以降のアクセスが高速化されます。一度にロードされるデータ量はCPUによって異なりますが、今回は僕のPCの仕様である64 byteを採用します。ちなみに、この値はキャッシュラインサイズと呼ばれます。
※今回の話はvector型では再現できません。vectorでは連続するメモリに行列のデータを保存しておらず、各行データを持つ配列のポインタを配列として持っており、この各行データはバラバラな位置に配置されている可能性があるためです。
※行列積の例では「Matrix2を転置すればよい」などの解決策がありますがOut of Scopeです。この記事の主目的は現象の解明・検証であり、サブ目的は「縦方向のアクセスをする際2のべき乗サイズの配列が危ない」というより大きい問題についてかみ砕いて説明することでそのために行列積を例として使用しています。もし行列積だと気になる場合は、横方向でも縦方向でもアクセスする必要があるような問題（画像処理など）に置き換えて考えてもらえれば幸いです。
行列計算
では、実際に行列計算のどの部分が重くなるのかを考えていきましょう。行列のある要素2行目2列目を計算するには以下のような計算を行う必要があります。
Matrix1では連続したメモリを読んでいけばいいので、1列目のデータを読み込めば2列目以降64 byteまではキャッシュを読むだけで調べられます。一方でMatrix2では1行目のデータを読み込んだとしても2行目はまだキャッシュにないデータであり、3行目、4行目でも同様に新たなキャッシュをロードしに行かないといけません。これをN回行う必要があります。
では次の行列の要素2行目3列目を計算する時を考えましょう。
先ほどと同様Matrix1は連続したメモリを読んでいきキャッシュからデータをとれるので問題はありません。一方でMatrix2については、少し事情が異なります。先ほどの計算で一度1行目からN行目まですべての行の最初の64 byteをキャッシュにロードしているので、先の計算時のキャッシュが残っていれば、今回の計算で再度メモリを探しに行く必要はなくなります。
つまり、1024個分のデータをキャッシュに保存することができれば、次のMatrix2を縦方向に辿るステップで再度データをロードする必要がなくなり早くなるということです。
ここまでをまとまとめると、何個分のキャッシュラインを保存することができるかできないかがボトルネックになりそうですね。
キャッシュライン衝突
それではキャッシュに保存できるのかを見ていきましょう。キャッシュサイズなどはCPU依存ですが、今回は僕のPCスペックを参照して、以下のようにします。
- L1D キャッシュ: 32KB/core, 8-way set associative, 64 sets, 64 byte line size
- L2 キャッシュ: 1MB/core, 16-way set associative, 1024sets, 64 byte line size
L1キャッシュが一番アクセスしやすくてL2キャッシュは少し遅いですが、ここまではコアごとのキャッシュであり共有ではないので十分早いそうです。なので、L2に入るかどうかを考えます。もしキャッシュの容量をフルに使えれば、1024行分のデータをすべて保存できます(L1Dキャッシュでは足りません)が、キャッシュはハッシュテーブルのようなものなので、全容量使い切れないケースが多いです。というわけで、キャッシュの中身をもう少し詳細に考えます。
セットアソシアティブ方式 (set associative)：ひとつのタグに対していくつかのデータを格納できる方式です。 n-way set associative, m sets といった形で呼ばれ、この場合m個のタグを持ちそれぞれn個までデータを格納できるキャッシュとなります。
今回の仕様である「1024セットの16-way set associative」 では、キャッシュタグを1024個持ち、各キャッシュタグに対し16個分の64byteデータを格納できるということになります。そしてタグは64byte分のデータに共通するアドレス(addr / 64)をキャッシュタグの最大値で割ったものになるので
tag(addr) = (addr / 64) % 1024
のように計算できます。この値は、Matrix2を縦方向に辿っていくとどうなるでしょう？N=1024のとき、addrは 4byte * 1024 ずつ増えていのくので、tagは16回ごとに同じ値を取ります。そしてL2キャッシュは16-wayなので256(=16*16)までしか保存できないということになります。これは1024個分のデータには足りません。
一方N=1023, 1025のときは1024の時と違って1024で綺麗に割り切れないので余りは異なる値を巡回することになり、キャッシュの容量をうまく使ってデータを保存しきることができます。
ここで全節の内容を思い出すと、
1024個分のデータをキャッシュに保存することができれば、次のMatrix2を縦方向に辿るステップで再度データをロードする必要がなくなり早い。
という話がありました。今節の結果を合わせるとN=1024のケースだけ遅くなることがわかりました。
それっぽい検証
計測時間の比較
メモリアクセス時間をものすごく適当に見積もります。
1回ずつの計算ステップで行っているのは以下です。
Answer[i][j] = Answer[i][j] + Matrix1[i][k] * Matrix2[k][j]
ここでのメモリ呼び出しは Read 3回, Write 1回 となっています。この中でMatrix2[k][j]のReadだけ遅くなりがちな要素となります。
ReadもWriteもだいたい同じ時間tで、キャッシュミスの時のみ10倍程度の時間10tがかかると見積もると、N=1024のときは毎回キャッシュをロードするので13t(=t+t+t+10t)、N=1023, 1025のときはだいたいすでにキャッシュに入っているデータで計算できるので4t(=t+t+t+t)の時間がかかるといえるでしょう。
これらを比較するとN=1024ではN=1023, 1025に比べて3倍強(~13t/4t)の時間がかかるという見積もりになりました。
ここで実際に行列計算を行って計測してみると
N = 1023: 1.07754s N = 1024: 3.32129s N = 1025: 1.02125s
という結果になりました。なんかそれっぽい！
キャッシュロード回数
キャッシュロードの回数をすごく適当に見積もります。
N = 1024のとき、すべてのステップでMatrix2のキャッシュがヒットせずロードが起きているので、だいたい10^9(~1024*1024*1024)程度のオーダーになりそうです。実際にはMatrix1や解行列のロードもあるので、これより少し多くなるでしょう。
一方N = 1023, 1025の時は、各要素の計算開始時には毎回キャッシュのロードをやり直す必要がありますが、各要素の計算中はMatrix1が横方向に64byte進むたびにロードしなおし、すなわち64(~1024/(64/4: int型のbyte数))回程度ロードすることになるので、だいたい10^8(~1024*1024*64)程度のオーダーになりそうです。実際にはこれより多くなるでしょう。
以下のコマンドでキャッシュのロードの回数を確認しました。
perf stat -e LLC-loads ファイル名
このとき
N = 1023: 133,622,632 回 N = 1024: 2,149,731,970 回 N = 1025: 87,673,800 回
という結果を得られました。なんかそれっぽい！
※今回はあくまでオーダー確認が目的なので他プロセスの影響を排除するなどの作業行っていません。
以上2つの検証からそれっぽい！ということがわかりました。
議論の余地はあると思うので、訂正・改善があり次第更新します。
Special Thanks to @hikalium.
