線形代数ライブラリの速度比較
C++ Advent Calender 2014 13日目の記事です．
最初はマイペースに作成している自作ライブラリ内で使用しているメタプログラミングに関する技術の解説を行おうと考えていましたが、この分野は他に詳しい方々が詳細なスライドや記事を作成してくださっているのと、時間の都合等々で断念しました．代わりに、Advent Calender 3日目の5mingame2さんの紹介にもあったEigenやboost.ublasといった線形代数ライブラリのコンテナ毎の速度比較結果をまとめます．
はじめに
vectorやmatrixを利用した線形代数の処理を行いたい場合、内積や行列の積といった基本的な処理は自分で容易に記述することができますが、LU分解やSVD等の自分で実装するのは少し気が引ける処理を利用したい際に、どこかに良い線形代数ライブラリが無いかな～と探すことと思います．この時、C++が好きな人は何よりもまず処理速度が気になることでしょう．そこで、今後このようなライブラリを利用したいという人の参考になれば幸いだと思い、稚拙ながら本記事を執筆します．
今回、調査対象のライブラリはC++では恐らく最も有名なEigenとublasという2つのライブラリです．これらは、有名なExpression Templateと呼ばれるメタプログラミング技法を利用して実装されており、これは計算した値が必要になる場面まで（計算式の）評価を遅延することで、普通に逐次評価する際に生じる一時変数を可能な限り無くして高速化をはかっています．ただし、コードの書き方によっては逆に速度が低下してしまうこともあるので注意が必要です（参考記事:Boost.uBLAS の遅延評価について）．また、両ライブラリ共にヘッダオンリーで手軽に利用できる点もオススメできる理由です．
実験設定
ベクトルを表すコンテナ型、および行列を表すコンテナ型に対して、ある処理タスクを開始して終了するまでの処理時間を計測します．
処理時間の計測に関しては、1度の計測だけでは偶然の影響で本来のパフォーマンスが測れない場合を考慮し、タスク毎に100回の計測を行ってその平均と標準偏差を求めました．
また、ベクトルや行列の疎密による変化も調べるため、全体における非ゼロ要素の割合が10％の場合(スパース)と90％の場合(非スパース)とに分けて測定を行いました．
実行環境を以下に示します．
- コンパイラ： Visual Studio 2015 Preview (v140), GCC（g++4.9.2）
- ライブラリ： boost 1.57.0, Eigen 3.2.3
Visual StudioはNDEBUGを指定してReleaseビルド、GCCは-DNDEBUG -O3を指定してコンパイルを実行．
ベクトル型の比較
計測を行うコンテナは以下の8種類です．相対的な理解をしやすくするために、ベースラインとしてstd::vectorと(影の薄い)std::valarrayの計測も行います．
- std::vector
- std::valarray
- Eigen::VectorXd
- Eigen::SparseVector
- ublas::vector
- ublas::mapped_vector ：内部実装に利用しているSTLのhash mapの性能に依存
- ublas::compressed_vector ：ランダムアクセスに強く、構築と追加・削除に弱い
- ublas::coordinate_vector ：追加・削除に強く、ランダムアクセスに弱い
この中で、Eigen::VectorXdとublas::vectorは特に密ベクトルと呼ばれ、一般的な配列と同じ構造をしています．一方、Eigen::SparseVector、ublas::mapped_vector、ublas::compressed_vector、ublas::coordinate_vectorは疎ベクトルと呼ばれ、値が存在する要素（非ゼロ要素）のみを保持した構造となっています．このため、大規模で欠損値の多いデータセット（推薦分野での評価値行列、グラフの隣接行列など）に対してはメモリ消費を大幅に抑えることができます．
実行するタスクは次の3種類です．
行列型の比較
計測を行うコンテナは以下の6種類です．
- Eigen::MatrixXd;
- Eigen::SparseMatrix
; - ublas::matrix
; - ublas::mapped_matrix
; - ublas::compressed_matrix
; - ublas::coordinate_matrix
;
ベクトルの場合と同様に、Eigen::MatrixXdとublas::matrixが密行列、残りが疎行列となっています．
実行するタスクは次の3種類です．
- ランダムアクセス ：任意の行・列へのアクセス（1万回）
- 積 ：行列同士の積（100×100行列）
- LU分解 ：与えられた行列を2つの三角行列に分解．逆行列や連立方程式の解を求める際に利用される手法の一種
ただし、LU分解に関してはublasの疎行列には対応していないようでした．また、Visual Studio環境ではEigenのコンパイルが上手くいかず断念しました．
実験結果
以下の表の見方としては、各セルの値は [非スパース] / [スパース] を表しており、左側の数値が処理時間の平均値、()内が100回の計測値の標準偏差を表しています．表1,3がVisual Studioでの結果、表2,4にGCCでの結果です．（ただし、GCCの計測は仮想マシン(Ubuntu)上で行っているため、Visual Studioとの単純な性能比較はあまり意味が無いことに留意）
[表1：result of vector (Visual Studio)]
| ランダムアクセス (μs) | イテレーション (ns) | 内積 (ns) | |
|---|---|---|---|
| std::vector | 518.8(10.27) / 517.7(9.60) | 74.33(110.1) / 78.29(131.8) | 72.84(137.1) / 56.28(124.3) |
| std::valarray | 505.4(18.39) / 508.2(12.65) | 32.41(90.18) / 33.15(94.30) | 33.37(114.9) / 36.41(103.7) |
| Eigen::VectorXd | 506.4(10.92) / 506.9(7.66) | 28.54(80.06) / 29.7(87.85) | 26.18(84.49) / 28.4(93.5) |
| Eigen::SparseVector | 1208(27.7) / 1070(11.39) | 29.93(87.0) / 26.48(89.79) | 36.35(114.8) / 36.34(114.8) |
| ublas::vector | 506.5(12.02) / 507.5(6.21) | 59936(5049) / 59861(4759) | 32.10(98.3) / 34.86(108.6) |
| ublas::mapped_vector | 1669(14.34) / 1127(16.84) | 874166(130327) / 191085(12161) | 888374(135080) / 194833(14858) |
| ublas::compressed_vector | 507(9.16) / 506(5.78) | 53120(737.7) / 5999(216.0) | 72245(2093) / 9383(1205) |
| ublas::coordinate_vector | 1446(18.99) / 1256(11.60) | 54020(9105) / 6122(912.2) | 73304(10612) / 9367(927.4) |
| ランダムアクセス (μs) | イテレーション (ns) | 内積 (ns) | |
|---|---|---|---|
| std::vector | 158.0(45.79) / 155.2(41.64) | 94.43(26.94) / 93.47(26.20) | 282.4(133.5) / 238.0(119.8) |
| std::valarray | 164.8(55.07) / 147.0(24.50) | 26.72(3.81) / 27.66(3.29) | 159.7(69.96) / 147.4(56.49) |
| Eigen::VectorXd | 168.5(66.26) / 148.5(28.11) | 28.6(4.16) / 30.44(4.52) | 160.3(69.66) / 168.3(75.14) |
| Eigen::SparseVector | 1180(378.7) / 856.0(167.0) | 33.46(2.76) / 30.08(2.42) | 160.2(69.31) / 147.4(55.19) |
| ublas::vector | 161.4(58.08) / 148.4(31.79) | 26.99(5.78) / 27.16(2.20) | 159.8(70.19) / 144.9(49.50) |
| ublas::mapped_vector | 1612(889.5) / 735.3(126.9) | 890425(252637) / 93699(19464) | 1872770(695327) / 700327(224015) |
| ublas::compressed_vector | 1187(403.1) / 875.3(196.1) | 181.6(54.17) / 143.2(47.26) | 173815(69871) / 16222(7674) |
| ublas::coordinate_vector | 1201(366.0) / 874.4(139.4) | 1055(9226) / 187.2(987.9) | 160998(70554) / 17240(13653) |
ベクトルの結果を見ると、Eigen::VectorXdとstd::valarrayが全体的に良く、std::vectorもそれに続く性能があり、ublas::vectorはイテレーションが必要な処理が遅いことが分かります．ubalsに関しては、疎ベクトルも含め最適化が上手く働いていないせい？か非常に遅い結果となりました（もし上述したコンパイラオプション以外にも何か最適化に影響するものを知っている方がいればご教授ください）．また、内積の結果を見ると標準偏差が他のタスクに比べて大きく、同じタスクでも処理時間のバラつきが大きいことが見て取れます．
[表3：result of matrix (Visual Studio)]
| ランダムアクセス (μs) | 積 (μs) | |
|---|---|---|
| Eigen::MatrixXd | 976.1(43.87) / 940.3(8.08) | 195.0(17.01) / 187.8(13.18) |
| Eigen::SparseMatrix | 1479(70.52) / 1244(12.47) | 1634(44.96) / 144.7(8.88) |
| ublas::matrix | 964.2(47.21) /930.1(12.97) | 994.1(23.09) / 978.8(18.10) |
| ublas::mapped_matrix | 2303(110.7) / 1706(10.47) | 158780(867.0) / 63955(456.3) |
| ublas::compressed_matrix | 964.2(47.49) / 927.7(7.55) | 46028(392.8) / 18441(202.6) |
| ublas::coordinate_matrix | 2581(119.0) / 1982(8.77) | 157322(1168) / 87679(657.5) |
[表4：result of matrix (GCC)]
| ランダムアクセス (μs) | 積 (μs) | LU分解 (μs) | |
|---|---|---|---|
| Eigen::MatrixXd | 452.3(127.3) / 443.7(72.73) | 198.7(17.19) / 196.8(9.77) | 151(10.86) / 150.5(5.86) |
| Eigen::SparseMatrix | 841.6(19.21) / 675.5(9.89) | 1625(53.3) / 156.4(11.94) | 578.0(18.53) / 484.7(10.59) |
| ublas::matrix | 441.4(19.02) / 458.5(148.1) | 862.9(33.57) / 854.9(37.72) | 216.6(8.11) / 212.8(4.83) |
| ublas::mapped_matrix | 1642(126.9) / 1042(161.1) | 141289(1145) / 61621(553.4) | - |
| ublas::compressed_matrix | 927.1(40.22) / 713.9(15.07) | 31976(680.0) / 15285(264.8) | - |
| ublas::coordinate_matrix | 1750(40.09) / 1346(28.33) | 102352(1322) / 15285(264.8) | - |
行列の結果もベクトルの結果と同様にEigen::matrixが全体的に優れた結果となりました.Eigen::SparseMatrixに関しては、行列の積の結果を見るとスパースなデータに対しては密なEigen::MatrixXdを越える性能が見られました．
追記
@TAIYAKI926さんのコメントで「#define BOOST_UBLAS_NDEBUG」というublasの最適化オプションがあるようなので、これを適用したものに修正しました．また、Eigenのバージョンも3.2.3に更新して再度実験を行いました．
変更後の顕著な変化としましては、mapped_matrixのランダムアクセス(非スパース)がVisualStudioで2267800μs→2303μs、GCCで1459870μs→1642μsに大幅改善されました(スパースも同様に改善)．
まとめ
- 基本的にEigenを使えば問題ない
- 疎ベクトル・疎行列に特化したコンテナは、密ベクトル・密行列に比べ、データが密な場合は速度が2倍程落ち、疎な場合は同等以上の速度である
- 影が薄いstd::valarrayが予想以上にデキる子だった
今回、実験に使用したコードはgithubで公開しています．
明日14日目はlapis_twさんです．よろしくお願いいたします．
