Writing Functions Taking Eigen Types as Parameters

Eigen's use of expression templates results in potentially every expression being of a different type. If you pass such an expression to a function taking a parameter of type Matrix, your expression will implicitly be evaluated into a temporary Matrix, which will then be passed to the function. This means that you lose the benefit of expression templates. Concretely, this has two drawbacks:

Fortunately, all this myriad of expression types have in common that they all inherit a few common, templated base classes. By letting your function take templated parameters of these base types, you can let them play nicely with Eigen's expression templates.

Table of contents

Some First Examples

This section will provide simple examples for different types of objects Eigen is offering. Before starting with the actual examples, we need to recapitulate which base objects we can work with (see also The class hierarchy).

EigenBase Example

Prints the dimensions of the most generic object present in Eigen. It coulde be any matrix expressions, any dense or sparse matrix and any array.

Example:Output: 
#include <iostream>
#include <Eigen/Core>
using namespace Eigen;

template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
  std::cout << "size (rows, cols): " << b.size() << " (" << b.rows()
            << ", " << b.cols() << ")" << std::endl;
}

int main()
{
    Vector3f v;
    print_size(v);
    // v.asDiagonal() returns a 3x3 diagonal matrix pseudo-expression
    print_size(v.asDiagonal());
}

 




 DenseBase Example 



 Prints a sub-block of the dense expression. Accepts any dense matrix or array expression, but no sparse objects and no special matrix classes such as DiagonalMatrix. 


template <typename Derived>
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
{
  std::cout << "block: " << b.block(x,y,r,c) << std::endl;
}

  ArrayBase Example 



 Prints the maximum coefficient of the array or array-expression. 


template <typename Derived>
void print_max(const ArrayBase<Derived>& a, const ArrayBase<Derived>& b)
{
  std::cout << "max: " << (a.max(b)).maxCoeff() << std::endl;
}

  MatrixBase Example 



 Prints the inverse condition number of the given matrix or matrix-expression. 


template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
  const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType&
    sing_vals = a.jacobiSvd().singularValues();
  std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}

These examples are just intended to give the reader a first impression of how functions can be written which take a plain and constant Matrix or Array argument. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold.



In which cases do functions taking plain Matrix or Array arguments work?


Let's assume one wants to write a function computing the covariance matrix of two input matrices where each row is an observation. The implementation of this function might look like this 


MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

 and contrary to what one might think at first, this implementation is fine unless you require a genric implementation that works with double matrices too and unless you do not care about temporary objects. Why is that the case? Where are temporaries involved? How can code as given below compile? 


MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);

 In this special case, the example is fine and will be working because both parameters are declared as const references. The compiler creates a temporary and evaluates the expression x+z into this temporary. Once the function is processed, the temporary is released and the result is assigned to C.


Note: Functions taking const references to Matrix (or Array) can process expressions at the cost of temporaries.



In which cases do functions taking a plain Matrix or Array argument fail?


Here, we consider a slightly modified version of the function given above. This time, we do not want to return the result but pass an additional non-const paramter which allows us to store the result. A first naive implementation might look as follows. 


// Note: This code is flawed!
void cov(const MatrixXf& x, const MatrixXf& y, MatrixXf& C)
{
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

 When trying to execute the following code 


MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));

 the compiler will fail, because it is not possible to convert the expression returned by MatrixXf::block() into a non-const MatrixXf&. This is the case because the compiler wants to protect you from writing your result to a temporary object. In this special case this protection is not intended -- we want to write to a temporary object. So how can we overcome this problem?


The solution which is preferred at the moment is based on a little hack. One needs to pass a const reference to the matrix and internally the constness needs to be cast away. The correct implementation for C98 compliant compilers would be 


template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;

  const Scalar num_observations = static_cast<Scalar>(x.rows());

  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;

  const_cast< MatrixBase<OtherDerived>& >(C) =
    (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

 The implementation above does now not only work with temporary expressions but it also allows to use the function with matrices of arbitrary floating point scalar types.


Note: The const cast hack will only work with templated functions. It will not work with the MatrixXf implementation because it is not possible to cast a Block expression to a Matrix reference!



How to resize matrices in generic implementations?


One might think we are done now, right? This is not completely true because in order for our covariance function to be generically applicable, we want the follwing code to work 


MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);

 This is not the case anymore, when we are using an implementation taking MatrixBase as a parameter. In general, Eigen supports automatic resizing but it is not possible to do so on expressions. Why should resizing of a matrix Block be allowed? It is a reference to a sub-matrix and we definitely don't want to resize that. So how can we incorporate resizing if we cannot resize on MatrixBase? The solution is to resize the derived object as in this implementation. 


template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C_)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;

  const Scalar num_observations = static_cast<Scalar>(x.rows());

  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;

  MatrixBase<OtherDerived>& C = const_cast< MatrixBase<OtherDerived>& >(C_);
  
  C.derived().resize(x.cols(),x.cols()); // resize the derived object
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

 This implementation is now working for parameters being expressions and for parameters being matrices and having the wrong size. Resizing the expressions does not do any harm in this case unless they actually require resizing. That means, passing an expression with the wrong dimensions will result in a run-time error (in debug mode only) while passing expressions of the correct size will just work fine.


Note: In the above discussion the terms Matrix and Array and MatrixBase and ArrayBase can be exchanged and all arguments still hold.



Summary











libicr

Author(s): Robert Krug

autogenerated on Mon Jan 6 2014 11:34:08