2015-12-20

高速なビット行列乗算

前置き

この記事Competitive Programming (その2) Advent Calendar 201512月20日の分です.12月19日はhadroriさんの競技プログラマー入門者用単語集12月21日はroiti46さんです.

皆さんこんばんは.Mです.Advent Calendarを書くのは初めてなのでドキドキしていますが,どうぞよろしくおねがいします.普段あまり記事を書かないので anond を使わせてもらっています

ここでは,ちょっと役立つ小ネタとして,今年書いたコードを1つ紹介します.「ビット行列を高速に乗算するコード」です.ごく簡単なコードですが定数倍効率効果が大きいので嘘解法に使えます

ビット行列

要素がブール値(True/False)であるような行列ビット行列と呼ぶことにします.このような行列に対する演算特に乗算)はアルゴリズム理論ではよく出てきます.最も有名な例はグラフの推移包閉(各点から行ける点を全部求める)です;行列 A をグラフの隣接行列とし,和をbit-or, 積をbit-and で定義すると,行列のべき乗和「A^0 + A^1 + ... + A^{n-1}」の (i,j) 要素が True であることと,j から i に到達可能であることが同値になります.なお,このべき乗和は (I + A)^{n-1} と等しいので,高速べき乗で一発で求めることが可能です.グラフの推移包閉を求める現在理論的に)最速の手法はこのアプローチに基づいており,計算量 O(n^\omega / log n) を達成します.O(n^\omega) は行列乗算の計算量で,現時点では \omega = 2.3728639... が最良です(Francois 2014).また,log n は Method of Four Russians と呼ばれるビット演算高速化する一般的テクニックで,サイズ log n までの演算結果を全部ハッシュに突っ込んでおくものです.

さて,この Method of Four Russian というテクニックは,実際に実装してもあまり早くないことが知られていますテーブル引きが遅い・単純なループが早い).ただし「いくつかのビットブロック単位計算する」というアイデア実用的にも有用です.ブロック単位演算ビット演算実装できるとき,そのアルゴリズムは「ビットパラレルアルゴリズム」と呼ばれています編集距離などの例が有名です.

ここでは,ビット行列に対するビットパラレル行列乗算を実装してみました.

コード&実測結果

A^n を計算するプログラムを書きました.実測結果を以下に示します.MacBook Pro; 2.8GHz Intel Core i7; 16GB 1600 MHz DDR3; g++ -std=c++11 -O3.実装http://ideone.com/8AsuI2 にあります

n提案手法[s]通常手法[s]
160.0000140.000082
320.0000250.000602
640.0000740.004485
1280.0004400.036760
2560.0027570.311192
5120.0201632.847787
10240.20055624.648609
20481.567657205.503351
409613.894987---
8196124.414617---

それなりに大きな n について 120倍くらい高速化しました.これだけ差があると嘘解法が通るようになります

解説

ビット行列を 8×8 のブロックに分割し,それぞれを unsigned long long (64bit) 1つで保存します.64が平方数というのが美しいですね.全体の乗算は8×8ブロックの乗算を普通に行えばよいので,結局 unsigned long longで表現された2つの8×8行列の積を考えれば十分です.

ここでは8×8行列の積を外積形式実装します.外積形式というのは C = A B という積を C = (Aの1列目×Bの1行目) + ... + (Aのn列目×Bのn行目) という外積の和の形で表現するものです.各外積は,すべての列がAのk列目と等しい8×8行列とすべての行がBのk行目と等しい8×8行列bit andに等しいので,Aからすべての列がAのk列目と等しい行列を作る方法とBからすべての行がBのk行目に等しい行列を作る方法を考えれば十分ですが,これはビットマスクして定数乗算すれば実装できます

余談

このコードとある問題に対する嘘解法用に作成したものですが,結局普通のほうでも通るようになってしまったので,オフィシャルにお披露目する機会がありませんでした.

記事への反応(ブックマークコメント)

ログイン ユーザー登録
ようこそ ゲスト さん