[eigen] generic argument declaration taking a matrix, vector or expression

[Bernhard Zeisl <bzeisl@xxxxxxxxxxxxxxxx>]

Hi all,

we are trying to implement functions, which can take a reference to a 
(templated) Eigen matrix, vector or an Eigen expression.
At http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html 
it's documented how to achieve this, however this approach is not 
compiling in our case.

Assume the following:

Eigen::MatrixXf A; // float matrix
Eigen::MatrixXi B; // integer matrix
Eigen::VectorXf c; // float vector

Particularly we would like to call a templated function foo with all of 
the following statements

foo (A)
foo (B)
foo (A.col(i));
foo (c)

such that no data is copied and results can be written back to the 
particular matrix (or column of the matrix).
If we understand the tutorial correctly, foo should look like

template <typename TypeT>
void foo (const Eigen::MatrixBase<TypeT> &mat);

where a const cast needs to be used to write data back to mat.

I attached a an example, which should make more clear what we are trying 
to achieve.
The compiler error we get for this example reads as follows:

make test_eigen_matrices
[100%] Building CXX object 
CMakeFiles/test_eigen_matrices.dir/tests/test_eigen_matrices.cpp.o
tests/test_eigen_matrices.cpp: In function ‘int main(int, char**)’:
tests/test_eigen_matrices.cpp:65:38: error: no matching function for 
call to ‘solveLeastSquares(Eigen::MatrixXf&, 
Eigen::DenseBase<Eigen::Matrix<float, -0x00000000000000001, 
-0x00000000000000001> >::ColXpr, Eigen::VectorXf&)’
tests/test_eigen_matrices.cpp:73:47: error: no matching function for 
call to ‘solveLeastSquares(Eigen::MatrixXf&, 
Eigen::DenseBase<Eigen::Matrix<float, -0x00000000000000001, 
-0x00000000000000001> >::ColXpr, Eigen::DenseBase<Eigen::Matrix<float, 
-0x00000000000000001, -0x00000000000000001> >::ColXpr)’
make[3]: *** 
[CMakeFiles/test_eigen_matrices.dir/tests/test_eigen_matrices.cpp.o] Error 1
make[2]: *** [CMakeFiles/test_eigen_matrices.dir/all] Error 2
make[1]: *** [CMakeFiles/test_eigen_matrices.dir/rule] Error 2
make: *** [test_eigen_matrices] Error 2

Can you please help us to understand how the correct declaration would 
look like.
Is such a generic declaration possible after all? If not, how can we 
avoid copying data and handle a matrix, vector or expression as possible 
argument.

Thanks for your help.

Cheers,
Bernhard

/*
 * test_eigen_matrices.cpp
 *
 *  Created on: Jan 10, 2012
 *      Author: bzeisl
 */

#include <iostream>
#include <Eigen/Dense>

enum LEAST_SQUARES_METHOD {
  SVD = 0,
  CHOLESKY,
  PSEUDO_INVERSE
};

// template <typename TypeT> void
// solveLeastSquares (const Eigen::Matrix<TypeT, Eigen::Dynamic, Eigen::Dynamic> &A, const Eigen::Matrix<TypeT, Eigen::Dynamic, Eigen::Dynamic> &b,
//                    Eigen::Matrix<TypeT, Eigen::Dynamic, Eigen::Dynamic> &x, LEAST_SQUARES_METHOD method = SVD)

template <typename TypeT> void
solveLeastSquares (const Eigen::MatrixBase<TypeT> &A, const Eigen::MatrixBase<TypeT> &b, Eigen::MatrixBase<TypeT> const &x_, LEAST_SQUARES_METHOD method = SVD)

{
  Eigen::MatrixBase<TypeT> &x = const_cast<Eigen::MatrixBase<TypeT> &> (x_);

  switch (method)
  {
    case SVD:
      x = A.jacobiSvd (Eigen::ComputeThinU | Eigen::ComputeThinV).solve (b);
      break;
    case CHOLESKY:
      x = A.ldlt ().solve (b);
      break;
    case PSEUDO_INVERSE:
      x = (A.transpose() * A).inverse() * A.transpose() * b;
      break;
    default:
      std::cout << "[solveLeastSquares] Least squares method (" << (int)method << ") not known." << std::endl;
      break;
  }
}

int
main (int argc, char* argv[])
{
  /** objective: solve the least squares system
   *  arg min_{x_i} 1/2 sum_i || b_i - A x_i ||_2^2  =  arg min_X 1/2 || B - A X ||_F^2
   */

  Eigen::MatrixXf A = Eigen::MatrixXf::Random (1000,200);
  Eigen::MatrixXf B = Eigen::MatrixXf::Random (1000,100);
  Eigen::MatrixXf X (200,100);


  // 1st test: solve for whole matrices
  solveLeastSquares (A, B, X); // compiles fine
  // expected result: X contains the solution


  // 2nd test: loop over data vector B
  Eigen::VectorXf x;
  for (int i = 0; i < B.cols(); ++i)
  {
    solveLeastSquares (A, B.col(i), x); // compile error
    // expected result: x contains the solution for the data sample b_i = B.col(i)
  }


  // 3rd test: loop over data matrix B and over result matrix X
  for (int i = 0; i < B.cols(); ++i)
  {
      solveLeastSquares (A, B.col(i), X.col(i)); // compile error
  }
  // expected result: the solution x_i for b_i = B.col(i) is written consecutively in the ith column in X .
  //                   finally X contains the result.

  return EXIT_SUCCESS;
}