先日C++でプログラムを作成していたところ、MATLABのdiff()のようなものが必要になったので、Eigenを使って同様の関数を実装してみた。
更新履歴
MATLABのdiff()関数とは
MATLABのdiff()関数は隣り合う要素の差分を返す関数。
行列に対しては、行間の差分を取る。m行n列の行列に対してdiff()を適用すると(m-1)行n列の行列が返される。
（列間の差分も取れるが）
また、diff(X, k)のように書くことで再帰的に複数回適用できる。
Y = diff(eye(5), 2) 1 -2 1 0 0 0 1 -2 0 0 0 0 1 -2 1
似たような関数はpythonのnumpyやRにもある。pythonのdiff()ではMATLABと動作が違い、隣り合う列の差を取ったものを返す。
一方、Rのdiff()はMATLABと同じで隣り合う行の差を取ったものを返すようだ。いずれも再帰的に複数回適用することも可能。
Eigenによる実装
Eigenを使い、C++でこれと似たような関数を実装する。
Eigenは式テンプレートという仕組みを使用していることから、Matrixクラスのオブジェクトを引数に取る関数を使う場合は、Matrix<...> &とするよりは、MatrixBase<Derived> &として、テンプレート関数にしたほうがよい。このようにすると、Matrix<>::Identity()や行列を演算したものなどを直接引数に取れる。
When writing a function taking Eigen objects as argument, if you want your function to take as argument any matrix, vector, or expression, just let it take a MatrixBase argument.
（訳）Eigenオブジェクトを引数として取る関数を作成するときに、関数に任意の行列、ベクトル、または式を引数として使用させたい場合は、MatrixBase引数をとるようにする。
（MatrixBaseクラスのリファレンスより）
実装にあたり、EigenのMatrixクラスには、行列の一部を取り出して部分行列を作るメソッドがあるので、これを使う。
初稿では、下記のように.topRows()と.bottomRows()を使うバージョンだったが、forループの中で.eval()が必要であり、毎回オブジェクトの生成とコピーが発生していたことから、効率が良くなかった（と思われる）。また、これだとVectorXdなども引数に取れるが、戻り値がMatrixXdになってしまう。
改訂版では、これを修正することにする。ついでにMatrixXiなどMatrixXd以外にも対応し、列間の差分も取れるようにした。
初稿版
#include <Eigen/Dense> template <typename Derived> const Eigen::MatrixXd Diff(const Eigen::MatrixBase<Derived> & m, const unsigned int n = 1) { Eigen::MatrixXd mr = m; for(unsigned int i = 0; i < n; i++) { mr = (mr.topRows(mr.rows() - 1) - mr.bottomRows(mr.rows() - 1)).eval(); } return mr; }
改訂版
実装方針として、密行列に対しては部分行列を作るメソッドは左辺値にも使うことができるので、部分行列で計算した結果を元の行列に戻すことで、オブジェクトの生成とコピーの発生を抑制する。
行間の差分を取る際に.topRows()はそのまま使うが、.bottomRows()ではなく、.block()で取り出す行数を変えるようにする。列間の差分を取る場合は、同様に.leftCols()と.block()を使用する。
なお、MatrixXd以外に対応させるために、Derived::PlainObjectを使う。
template <typename Derived> const typename Derived::PlainObject Diff(const Eigen::MatrixBase<Derived> & m0, const int n = 1, const int dim = 0) { typename Derived::PlainObject m = m0; int rows = m.rows(), columns = m.cols(); // dim == 1で列間の差分を取る if(dim == 1) { for(int i = 0; i < n ; i++) { int col_counts = columns - i - 1; m.leftCols(col_counts) = m.leftCols(col_counts) - m.block(0, 1, rows, col_counts); } return m.leftCols(columns - n); } else { for(int i = 0; i < n ; i++) { int row_counts = rows - i - 1; m.topRows(row_counts) = m.topRows(row_counts) - m.block(1, 0, row_counts, columns); } return m.topRows(rows - n); } }
ところで、残念なことにSparseMatrixに対しては左辺値に使えない様子なので、SparseMatrix版の実装は初稿版と同じ方針を取ることにした。
（なにか良い方法はないものか）
SparseMatrix版まで実装したソースコードはdiff.hとしてGitHubに公開した。Eigenと同じようにヘッダファイルのインクルードのみで動作する。
使用例
#include <iostream> #include "diff.h" int main() { Eigen::MatrixXf mf(4, 4); Eigen::MatrixXi mi = Eigen::MatrixXi::Identity(5, 5); Eigen::VectorXd vd(5); Eigen::RowVectorXd rvd(5); Eigen::SparseMatrix<double> sm; mf << -2, 1, 0, 3, 3, 2, -2, 3, 5, 0, -1, -3, 1, -1, 4, 1; vd << 0, 3, 2, 2, 4; rvd << 5, 3, 2, 0, 1; std::cout << "Diff(Eigen::MatrixXd::Identity(5, 5), 2) = " << std::endl; std::cout << Diff(Eigen::MatrixXd::Identity(5, 5), 2) << std::endl << std::endl; std::cout << "Diff(mf, 2, 1) = " << std::endl; std::cout << Diff(mf, 2, 1) << std::endl << std::endl; std::cout << "Diff(mi - Eigen::MatrixXi::Ones(5, 5), 3) = " << std::endl; std::cout << Diff(mi - Eigen::MatrixXi::Ones(5, 5), 3) << std::endl << std::endl; std::cout << "Diff(vd) = " << std::endl; std::cout << Diff(vd) << std::endl << std::endl; std::cout << "Diff(rvd, 1, 1) = " << std::endl; std::cout << Diff(rvd, 1, 1) << std::endl << std::endl; sm = Eigen::VectorXd::Ones(5).asDiagonal(); std::cout << "Diff(sm, 2) = " << std::endl; std::cout << Diff(sm, 2) << std::endl; return 0; }
出力結果
Diff(Eigen::MatrixXd::Identity(5, 5), 2) = 1 -2 1 0 0 0 1 -2 1 0 0 0 1 -2 1 Diff(mf, 2, 1) = -4 4 -3 9 4 -1 7 -8 Diff(mi - Eigen::MatrixXi::Ones(5, 5), 3) = 1 -3 3 -1 0 0 1 -3 3 -1 Diff(vd) = -3 1 0 -2 Diff(rvd, 1, 1) = 2 1 2 -1 Diff(sm, 2) = Nonzero entries: (1,0) (-2,0) (1,1) (1,0) (-2,1) (1,2) (1,1) (-2,2) (1,2) Outer pointers: 0 1 3 6 8 $ 1 -2 1 0 0 0 1 -2 1 0 0 0 1 -2 1