Eigen » Extending/Customizing Eigen » Matrix manipulation via nullary-expressions

Example 1: circulant matrix

To explore these possibilities let us start with the circulant example of the implementing new expression topic. Let us recall that a circulant matrix is a matrix where each column is the same as the column to the left, except that it is cyclically shifted downwards. For example, here is a 4-by-4 circulant matrix:

A circulant matrix is uniquely determined by its first column. We wish to write a function makeCirculant which, given the first column, returns an expression representing the circulant matrix.

For this exercise, the return type of makeCirculant will be a CwiseNullaryOp that we need to instantiate with: 1 - a proper circulant_functor storing the input vector and implementing the adequate coefficient accessor operator(i,j) 2 - a template instantiation of class Matrix conveying compile-time information such as the scalar type, sizes, and preferred storage layout.

Calling ArgType the type of the input vector, we can construct the equivalent squared Matrix type as follows:

template<class ArgType>
struct circulant_helper {
  typedef Matrix<typename ArgType::Scalar,
                 ArgType::SizeAtCompileTime,
                 ArgType::SizeAtCompileTime,
                 ColMajor,
                 ArgType::MaxSizeAtCompileTime,
                 ArgType::MaxSizeAtCompileTime> MatrixType;
};

This little helper structure will help us to implement our makeCirculant function as follows:

template <class ArgType>
CwiseNullaryOp<circulant_functor<ArgType>, typename circulant_helper<ArgType>::MatrixType>
makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
  typedef typename circulant_helper<ArgType>::MatrixType MatrixType;
  return MatrixType::NullaryExpr(arg.size(), arg.size(), circulant_functor<ArgType>(arg.derived()));
}

As usual, our function takes as argument a MatrixBase (see this page for more details). Then, the CwiseNullaryOp object is constructed through the DenseBase::NullaryExpr static method with the adequate runtime sizes.

Then, we need to implement our circulant_functor, which is a straightforward exercise:

template<class ArgType>
class circulant_functor {
  const ArgType &m_vec;
public:
  circulant_functor(const ArgType& arg) : m_vec(arg) {}

  const typename ArgType::Scalar& operator() (Index row, Index col) const {
    Index index = row - col;
    if (index < 0) index += m_vec.size();
    return m_vec(index);
  }
};

We are now all set to try our new feature:

int main()
{
  Eigen::VectorXd vec(4);
  vec << 1, 2, 4, 8;
  Eigen::MatrixXd mat;
  mat = makeCirculant(vec);
  std::cout << mat << std::endl;
}

If all the fragments are combined, the following output is produced, showing that the program works as expected:

1 8 4 2
2 1 8 4
4 2 1 8
8 4 2 1

This implementation of makeCirculant is much simpler than defining a new expression from scratch.

Example 2: indexing rows and columns

The goal here is to mimic MatLab's ability to index a matrix through two vectors of indices referencing the rows and columns to be picked respectively, like this:

A =
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9

A([1 2 1], [3 2 1 0 0 2]) =
 0  1 -6 -2 -2  1
 9  0 -3  6  6  0
 0  1 -6 -2 -2  1

To this end, let us first write a nullary-functor storing references to the input matrix and to the two arrays of indices, and implementing the required operator()(i,j):

template<class ArgType, class RowIndexType, class ColIndexType>
class indexing_functor {
  const ArgType &m_arg;
  const RowIndexType &m_rowIndices;
  const ColIndexType &m_colIndices;
public:
  typedef Matrix<typename ArgType::Scalar,
                 RowIndexType::SizeAtCompileTime,
                 ColIndexType::SizeAtCompileTime,
                 ArgType::Flags&RowMajorBit?RowMajor:ColMajor,
                 RowIndexType::MaxSizeAtCompileTime,
                 ColIndexType::MaxSizeAtCompileTime> MatrixType;

  indexing_functor(const ArgType& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
    : m_arg(arg), m_rowIndices(row_indices), m_colIndices(col_indices)
  {}

  const typename ArgType::Scalar& operator() (Index row, Index col) const {
    return m_arg(m_rowIndices[row], m_colIndices[col]);
  }
};

Then, let's create an indexing(A,rows,cols) function creating the nullary expression:

template <class ArgType, class RowIndexType, class ColIndexType>
CwiseNullaryOp<indexing_functor<ArgType,RowIndexType,ColIndexType>, typename indexing_functor<ArgType,RowIndexType,ColIndexType>::MatrixType>
mat_indexing(const Eigen::MatrixBase<ArgType>& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
{
  typedef indexing_functor<ArgType,RowIndexType,ColIndexType> Func;
  typedef typename Func::MatrixType MatrixType;
  return MatrixType::NullaryExpr(row_indices.size(), col_indices.size(), Func(arg.derived(), row_indices, col_indices));
}

Finally, here is an example of how this function can be used:

  Eigen::MatrixXi A = Eigen::MatrixXi::Random(4,4);
  Array3i ri(1,2,1);
  ArrayXi ci(6); ci << 3,2,1,0,0,2;
  Eigen::MatrixXi B = mat_indexing(A, ri, ci);
  std::cout << "A =" << std::endl;
  std::cout << A << std::endl << std::endl;
  std::cout << "A([" << ri.transpose() << "], [" << ci.transpose() << "]) =" << std::endl;
  std::cout << B << std::endl;

This straightforward implementation is already quite powerful as the row or column index arrays can also be expressions to perform offsetting, modulo, striding, reverse, etc.

  B =  mat_indexing(A, ri+1, ci);
  std::cout << "A(ri+1,ci) =" << std::endl;
  std::cout << B << std::endl << std::endl;
#if __cplusplus >= 201103L
  B =  mat_indexing(A, ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3));
  std::cout << "A(ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3)) =" << std::endl;
  std::cout << B << std::endl << std::endl;
#endif

and the output is:

A(ri+1,ci) =
 9  0 -3  6  6  0
 9  3  6  6  6  3
 9  0 -3  6  6  0

A(ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3)) =
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3