The documentation for Writing Functions Taking Eigen Types as Parameters (http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html) nicely describes how to write functions that take generic Eigen types as input parameters and output parameters. For example `cov(x,y,C)` computes the covariance matrix of the inputs `x` and `y` into the output `C`.

I'd like to extend upon this example to consider generic Eigen types as return arguments. Then I'd like to be able to do things like `C = cov(x,y)` and `C.block(0,0,3,3) = cov(x,x+y);` etc.

Here's a little example program that contains a fully templated input-output implementation of cov (e.g. `cov(x,y,C)`) and then a broken wrapper with `C` as a return argument. In the `main` function I included some example calls to check compilation. Obviously I'd like to implement `cov` only once, so hopefully the return-arg version is just a thin wrapper on top of the full output-arg version. Ideally, I'd like things like `C.block(0,0,3,3) = cov(x,y);` to work too.

Code: Select all

using namespace Eigen;
template <typename Derivedx, typename Derivedy, typename OtherDerived>
void cov(const MatrixBase<Derivedx>& x, const MatrixBase<Derivedy>& y, MatrixBase<OtherDerived> const & C_)
{
  typedef typename Derivedx::Scalar Scalar;
  typedef typename internal::plain_row_type<Derivedx>::type RowVectorTypex;
  typedef typename internal::plain_row_type<Derivedy>::type RowVectorTypey;
  const Scalar num_observations = static_cast<Scalar>(x.rows());
  const RowVectorTypex x_mean = x.colwise().sum() / num_observations;
  const RowVectorTypey 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;
}

template <typename Derivedx, typename Derivedy>
WhatShouldThisReturnTypeBe?
cov(const MatrixBase<Derivedx>& x, const MatrixBase<Derivedy>& y)
{
  WhatShouldThisReturnTypeBe? C;
  cov(x,y,C);
  return C;
}

int main()
{
  MatrixXf x = MatrixXf::Random(100,3);
  MatrixXf y = MatrixXf::Random(100,3);
  {
    MatrixXf C;
    cov(x,y,C);
  }
  {
    MatrixXf C;
    cov(x,x+y,C);
  }
  {
    MatrixXf C = MatrixXf::Zero(3,6);
    cov(x,y, C.block(0,0,3,3));
  }
  {
    MatrixXf C = cov(x,x+y);
  }
  {
    MatrixXf C = MatrixXf::Zero(3,6);
    C.block(0,0,3,3) = cov(x,y);
  }
}


The documentation for this is on dev version:
http://eigen.tuxfamily.org/dox-devel/To ... nType.html

Even known it's in dev, I think that the latest release support this api already, so you are good to go.

Instead of writing a full expression, one can exploit Eigen::ReturnByValue to reduce a bit the amount of work. In both cases, the general idea is to let cov(x,y) returns a proxy object, let's call it CovProxy storing references to x and y, and from which Eigen will call CovProxy::evalTo(C) into which you call cov(this->x, this->y, C);

However, even using Eigen::ReturnByValue, this still require some amount of boiler-plate code. Perhaps, this could be wrapped into a more generic proxy simply storing a pointer to the arguments, and a function pointer.

For the fun, I've started a proof of concept in C++11:

File "rvo_plugin.h":

Code: Select all

template<typename T>
Matrix& operator=(const std::function<void(T)>& func) {
  func(derived());
  return derived();
}


User code:

Code: Select all

#define EIGEN_MATRIX_PLUGIN "rvo_plugin.h"

#include <Eigen/Core>
#include <iostream>
#include <functional>
using namespace Eigen;

void cov_inplace(const MatrixXd &X, const MatrixXd& Y, MatrixXd &C) {
  std::cout << "Evaluate cov(X,Y) into C\n";
}

std::function<void(MatrixXd&)> cov(const MatrixXd &X, const MatrixXd &Y)  {
  return std::bind(&cov_inplace, X, Y, std::placeholders::_1);
}

int main() {
  MatrixXd A, B, C;
  A = cov(B,C);
}


This example is still very limited, but it can surely be extended to be more generic on both the arguments and destination.