Implementation of MATLAB-like diff() in C++ using Eigen

2024-06-21C++,Eigen

Recently, while creating a program in C++, I needed something similar to MATLAB’s diff() function, so I implemented a similar function using Eigen.

What is MATLAB’s diff() function?

The MATLAB diff() function returns the difference between adjacent elements. For matrices, it computes the difference between rows. When diff() is applied to an $m \times n$ matrix, it returns an $(m-1) \times n$ matrix.

Additionally, by writing diff(X, k), the function can be applied recursively multiple times.

Y = diff(eye(5), 2)

1 -2  1  0  0
0  1 -2  0  0
0  0  1 -2  1

Similar functions are available in Python’s NumPy and R. Python’s diff(), however, behaves differently from MATLAB’s, as it returns the differences between adjacent columns, so caution is needed. On the other hand, R’s diff() works the same way as MATLAB’s, returning the differences between adjacent rows. Both can also be applied recursively multiple times.

Implementation in C++ using Eigen

Using Eigen, I’ll implement a function similar to MATLAB’s diff() in C++.

Eigen utilizes expression templates. When creating a function that takes a Matrix class object as an argument, it is better to use MatrixBase<Derived> & in a template function rather than Matrix<...> &. This allows the function to directly accept expressions like Matrix<>::Identity() or operations on matrices.

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.
（Reference for the MatrixBase class）

Eigen’s Matrix class has methods to extract submatrices, which we can use to improve efficiency.

Initially, I used .topRows() and .bottomRows() to create submatrices. However, this required calling .eval() inside the for loop, which involves object creation and copying every time, potentially reducing efficiency. Additionally, when passing VectorXd as an argument, the return type becomes MatrixXd.

In the revised version, I’ll address these issues. The updated implementation will also handle matrices like MatrixXi and allow for computing column-wise differences.

Initial Version

#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;
}

Revised Version

For dense matrices, methods to create submatrices can be used on lvalues, allowing us to place the results of calculations back into the original matrix. This approach helps to avoid object creation and copying.

When computing row-wise differences, .topRows() is used directly, but instead of using .bottomRows(), .block() is used to vary the number of rows extracted. Similarly, for column-wise differences, .leftCols() and .block() are used.

To support types other than MatrixXd, Derived::PlainObject is utilized.

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();

    // column-wise differences when 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);
    }
}

Unfortunately, it seems that you cannot use methods to create submatrices on the left-hand side for SparseMatrix. Therefore, for the SparseMatrix version, I have chosen to follow the same approach as the initial version.

The source code implementing the SparseMatrix version is available as diff.h on GitHub. It operates with just header file inclusion, similar to Eigen.

Usage Example

#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;
}

Results

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

Verified Environment

  • 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