【C++】【Eigen】MATLABのdiff()関数（っぽいもの）を実装する

2021-02-092024-04-16C++,Eigen

先日C++でプログラムを作成していたところ、MATLABのdiff()のようなものが必要になったので、Eigenを使って同様の関数を実装してみました。

更新履歴

22.03.25　C++のソースコードを修正。
22.01.24　文体を修正。
21.02.15　Eigenによる実装部分の文章を加筆・修正。合わせてC++のソースコードも修正。
21.02.09　初稿

1. MATLABのdiff()関数とは
2. Eigenによる実装

MATLABのdiff()関数とは

MATLABのdiff()関数は隣り合う要素の差分を返す関数です。
行列に対しては、行間の差分を取ります。m行n列の行列に対してdiff()を適用すると(m-1)行n列の行列が返されます。
（列間の差分も取れますが）

差分と近似微分 – MATLAB diff – MathWorks 日本

また、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.bottomRows(mr.rows() - 1) - mr.topRows(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.block(0, 1, rows, col_counts) - m.leftCols(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.block(1, 0, row_counts, columns) - m.topRows(row_counts);
        }

        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

動作を確認した環境

OS
1. Windows 10
  1. Microsoft C/C++ Optimizing Compiler 19.00.24210 for x86
  2. MinGW-W64 GCC version 8.1.0
2. Linux
  1. GCC version 9.3.0
Eigen version 3.3.9

C++,Eigen

Posted by izadori

qpdfviewでパスワード付きPDFが印刷できない