c++の線形代数ライブラリのEigenをよく使うのですが,前々から小さなサイズの動的行列の逆行列を求める処理が重いことが気になっていました.
以下のページを見て,てっきり動的行列でもサイズが小さければ公式を使って逆行列を求めるのかと思っていたのですが,よくよく見ると固定サイズでないと公式を使った計算は行わないと書いていて,実装を見てみると動的行列では必ず(1×1を除いて) partialPivLu を使って逆行列を求めているようでした.
http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
LinearSolverのsolve()を使えとか,サイズが小さい時は固定サイズを使えということなんでしょうが,やっぱりプログラム上どうしても動的行列の逆行列が欲しい場面はあるわけで,そのようなケースのために次のようなコードを書いてみました.
// fast_inverse template<typename T> void fast_inverse(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& m, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& dst) { switch (m.rows()){ case 2: dst = Eigen::internal::inverse_impl<Eigen::Matrix<T, 2, 2>>(m); break; case 3: dst = Eigen::internal::inverse_impl<Eigen::Matrix<T, 3, 3>>(m); break; case 4: dst = Eigen::internal::inverse_impl<Eigen::Matrix<T, 4, 4>>(m); break; default: dst = m.inverse(); break; } }
行列のサイズを見て,サイズが4以下なら固定サイズ用の逆行列計算ルーチンを呼び,そうでなければデフォルトのinverseメソッドを呼びます.
本当はEigenのExpressionTemplateに従って式を返すような関数にしたかったのですが,EigenのReturnByValueクラスまわりがよくわからなかったので,とりあえず出力行列を引数にとってそれに代入するようにしています.(だれかいい感じに式を返す方法を知っていれば教えて下さい)
で,こんな感じのコードで普通のEigen::Matrix4f::inverse()とEigen::MatrixXf::inverse()と速度比較を行いました.実行環境はi7-2640Mで,VisualStudio2013のRelseaseの最適化O2,拡張命令はAVXまで有効にしています.
int main(int argc, char** argv) { Eigen::VectorXf vec = Eigen::VectorXf::Random(4); Eigen::MatrixXf mat = vec * vec.transpose(); Eigen::MatrixXf inv; for (int i = 0; i < 10000; i++) { for (int j = 0; j < 10000; j++) { // mat.inverse(); fast_inverse( mat, dst ); } } }
結果は以下のとおり
Matrix4f::inverse : 2266[msec]
MatrixXf::inverse : > 30000[msec]
fast_inverse : 2941[msec]
やっぱり,MatrixXf::inverse() はいちいちLU分解しているので,小さな行列に対してはかなり遅いみたいです.それに対して fast_inverse() では固定サイズに近い程度に高速化できました.
もちろん,MatrixXf::inverse()に対してかなりアンフェアな比較で,サイズが4以上の行列に対してはオーバーヘッドが増えるだけなのですが,予め行列のサイズが小さいと見込める場合には有効な手段なんじゃないかと思います.