Matrix manipulation via nullary-expressions

The main purpose of the class CwiseNullaryOp is to define procedural matrices such as constant or random matrices as returned by the Ones(), Zero(), Constant(), Identity() and Random() methods. Nevertheless, with some imagination it is possible to accomplish very sophisticated matrix manipulation with minimal efforts such that implementing new expression is rarely needed.

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:

\[ \begin{bmatrix} 1 & 8 & 4 & 2 \\ 2 & 1 & 8 & 4 \\ 4 & 2 & 1 & 8 \\ 8 & 4 & 2 & 1 \end{bmatrix} \]

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>
typedef Matrix<typename ArgType::Scalar,
ArgType::SizeAtCompileTime,
ArgType::SizeAtCompileTime,
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>
{
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>
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;
std::cout << mat << std::endl;
}

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

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:

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>
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 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 EIGEN_COMP_CXXVER >= 11
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:

col
m col(1)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:66
B
Definition: test_numpy_dtypes.cpp:299
makeCirculant
CwiseNullaryOp< circulant_functor< ArgType >, typename circulant_helper< ArgType >::MatrixType > makeCirculant(const Eigen::MatrixBase< ArgType > &arg)
Definition: make_circulant2.cpp:36
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
main
int main(int argc, char **argv)
Definition: cmake/example_cmake_find_gtsam/main.cpp:63
ceres::Matrix
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > Matrix
Definition: gtsam/3rdparty/ceres/eigen.h:42
A
Definition: test_numpy_dtypes.cpp:298
operator()
internal::enable_if< internal::valid_indexed_view_overload< RowIndices, ColIndices >::value &&internal::traits< typename EIGEN_INDEXED_VIEW_METHOD_TYPE< RowIndices, ColIndices >::type >::ReturnAsIndexedView, typename EIGEN_INDEXED_VIEW_METHOD_TYPE< RowIndices, ColIndices >::type >::type operator()(const RowIndices &rowIndices, const ColIndices &colIndices) EIGEN_INDEXED_VIEW_METHOD_CONST
Definition: IndexedViewMethods.h:73
unaryExpr
const EIGEN_DEVICE_FUNC CwiseUnaryOp< CustomUnaryOp, const Derived > unaryExpr(const CustomUnaryOp &func=CustomUnaryOp()) const
Apply a unary operator coefficient-wise.
Definition: CommonCwiseUnaryOps.h:135
arg
Definition: cast.h:1277
row
m row(1)
mat_indexing
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)
Definition: nullary_indexing.cpp:33
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
indexing_functor
Definition: nullary_indexing.cpp:8
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:319
circulant_functor
Definition: make_circulant2.cpp:8
circulant_helper
Definition: make_circulant2.cpp:23
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Wed May 15 2024 15:27:19