Eigenを使って疎行列の計算をする

行列計算ライブラリEigenには密行列だけではなく、疎行列の計算を行うためのルーチンも含まれている。
（Eigenの紹介についてはこちらを参照 C++行列計算ライブラリEigen入門）

疎行列とはほとんどの要素が０で、非ゼロの要素が小数個しかない行列。
全行列要素をメモリに置くかわりに、非ゼロ要素の位置と値だけを保持する。密行列表現ではメモリに乗らないような巨大な行列も取り扱うことができる。

行列要素のセット

一番シンプルなコードは以下のようになる。

sparse.cpp

#include <iostream>
#include <vector>
#include <Eigen/Sparse>

typedef Eigen::Triplet<double> T;

int main() {
  std::vector<T> tripletVec;

  tripletVec.push_back( T(0,0,0.2) );
  tripletVec.push_back( T(0,1,0.5) );
  tripletVec.push_back( T(1,0,1.0) );
  tripletVec.push_back( T(1,1,2.0) );

  Eigen::SparseMatrix<double> M(2,2);
  M.setFromTriplets(tripletVec.begin(), tripletVec.end());
  std::cout << M << std::endl;

  return 0;
}

Nonzero entries:
(0.2,0) (1,1) (0.5,0) (2,1)

Outer pointers:
0 2  $

0.2 0.5
1 2

これで $\begin{pmatrix} 0.2 & 0.5 \\ 1 & 2 \ \end{pmatrix}$ という2x2行列を表現している。（この例は全く疎行列ではないが。）

行列要素をセットしていくには Eigen::Triplet の配列を作る。Tripletのコンストラクタは３つの引数を取り（行index、列index、要素の値）を表している。
Tripletのvectorを作ったら、M.setFromTriplets()を呼んで要素に値をセットすることができる。

詳細は http://eigen.tuxfamily.org/dox/group__TutorialSparse.html を参照のこと。

このようにして作った疎行列に対しては密行列とどうように +,-,* などの演算子が適切に定義されている。だいたい直感的に使えるようになっている。
ただし、計算コストがかかる場合があるので一応注意が必要。

密行列との演算も可能。

  Eigen::MatrixXd D(2,2);
  D << 3, 2,
       1, 4;
  std::cout << "M*D: " << M*D << std::endl;

連立一次方程式を解く

疎行列に対するSolverもEigenにいくつか用意されている。
http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html

どのmoduleを使って解くにしろ、以下のようなコードを書くことになる。

#include <Eigen/RequiredModuleName>
// ...
SparseMatrix<double> A;
// fill A
VectorXd b, x;
// fill b
// solve Ax = b
SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Success) {
  // decomposition failed
  return;
}
x = solver.solve(b);
if(solver.info()!=Success) {
  // solving failed
  return;
}

具体的なコードは以下の通り。例えばQR分解を使って解く場合。

sparse_solve.cpp

#include <iostream>
#include <vector>
#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <Eigen/SparseQR>

typedef Eigen::Triplet<double> T;

int main() {
  std::vector<T> tripletVec;

  tripletVec.push_back( T(0,0,0.2) );
  tripletVec.push_back( T(0,1,0.5) );
  tripletVec.push_back( T(1,0,1.0) );
  tripletVec.push_back( T(1,1,2.0) );

  Eigen::SparseMatrix<double> M(2,2);
  M.setFromTriplets(tripletVec.begin(), tripletVec.end());
  std::cout << M << std::endl;

  Eigen::SparseMatrix<double> b(2,1);
  std::vector<T> tripletVecb;
  tripletVecb.push_back( T(0,0,1.0) );
  tripletVecb.push_back( T(1,0,2.0) );  // b = (1.0 2.0)
  b.setFromTriplets(tripletVecb.begin(), tripletVecb.end());
  std::cout << b << std::endl;

  Eigen::VectorXd x;  // 解のベクトル

  Eigen::SparseQR< Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > solver;  // solverオブジェクトを構築する。
  solver.compute(M);

  if( solver.info() != Eigen::Success ) {
    std::cerr << "decomposition failed" << std::endl;
  }
  x = solver.solve(b);
  if( solver.info() != Eigen::Success ) {
    std::cerr << "solving failed" << std::endl;
  }
  std::cout << x << std::endl;

  return 0;
}