CODEWITHSUNDEEP
c++templatesmatrixtypedefeigen
marcman
marcman

Reputation: 3383

Passing fixed-size Eigen matrices as arguments to function calling for dynamic-size matrices

I am writing a small linear algebra utility library on top of Eigen for my personal code-base. To try to make it as flexible as possible, I typedefed different Eigen matrix types for use as parameters. However, an issue I keep running into is that when I use it, I can't pass a fixed-size (i.e. set at compile-time) matrix as the argument to a function that has a dynamically-sized (set at runtime) matrix typedef as a parameter. I could understand the inverse--not being able to pass a dynamic-sized matrix as fixed due to compile-time checks, but it seems this ought to work.

A test-able example is the pdist2 function below (which really ought to have a native implementation in the Eigen API).

#include <Eigen/Core>

namespace Eigen
{
    template <typename T>
    using MatrixXT = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
}

// X is M x N
// Y is M x K
// Output is N x K
template <typename T>
inline Eigen::MatrixXT<T> pdist2(const Eigen::MatrixXT<T> &X, const Eigen::MatrixXT<T> &Y)
{
    // ASSERT(X.rows() == Y.rows(), "Input observations must have same number of rows (" + 
    //  std::to_string(X.rows()) + "!=" + std::to_string(Y.rows()) + ")");

    Eigen::MatrixXT<T> dists = X.colwise().squaredNorm().transpose() * Eigen::MatrixXT<T>::Ones(1, Y.cols()) +
        Eigen::MatrixXT<T>::Ones(X.cols(), 1) * Y.colwise().squaredNorm() -
        2 * X.transpose() * Y;

    return dists;
}

This code does not compile:

Eigen::Matrix<double, 3, 5> X;
X << 8.147236863931790, 9.133758561390193, 2.784982188670484, 9.648885351992766, 9.571669482429456,
     9.057919370756192, 6.323592462254095, 5.468815192049838, 1.576130816775483, 4.853756487228412,
     1.269868162935061, 0.975404049994095, 9.575068354342976, 9.705927817606156, 8.002804688888002;

Eigen::Matrix<double, 3, 4> Y;
Y << 1.418863386272153, 7.922073295595544, 0.357116785741896, 6.787351548577734,
     4.217612826262750, 9.594924263929030, 8.491293058687772, 7.577401305783335,
     9.157355251890671, 6.557406991565868, 9.339932477575505, 7.431324681249162;

Eigen::Matrix<double, 5, 4> D = pdist2(X, Y);

The above function has been unit tested and evaluates correctly, but it will only work if X and Y are Eigen::MatrixXd types. It seems it must be my template typedef causing the issue, but it's just a dynamically (i.e. at runtime) sized matrix with a templated type.

The error reads:

error: no matching function for call to ‘pdist2(Eigen::Matrix<double, 3, 5>&, Eigen::Matrix<double, 3, 4>&)’
  Eigen::Matrix<double, 5, 4> D = Util::Math::pdist2(X, Y);
                                                              ^
note: candidate: template<class T> Eigen::MatrixXT<T> Util::Math::pdist2(Eigen::MatrixXT<T>&, Eigen::MatrixXT<T>&)
   inline Eigen::MatrixXT<T> pdist2(const Eigen::MatrixXT<T> &X, const Eigen::MatrixXT<T> &Y)
                         ^
note:   template argument deduction/substitution failed:
note:   template argument ‘3’ does not match ‘#‘integer_cst’ not supported by dump_decl#<declaration error>’
  Eigen::Matrix<double, 5, 4> D_est = Util::Math::pdist2(X, Y);
                                                              ^
note:   ‘Eigen::Matrix<double, 3, 5>’ is not derived from ‘Eigen::MatrixXT<T>’

Why is this not working? Or perhaps more specifically, how do I use a templated typedef to ensure that my fixed-size matrices do derive from Eigen::MatrixXT<T>?

Note: This is all using Eigen 3.3.3.

Upvotes: 5

Views: 3028

Answers (1)

ggael
ggael

Reputation: 29205

The problem is that Matrix<double, 3, 5> and MatrixXT<double> are not the same type, therefore, the only way to pass the former to pdist2 is to convert it to a MatrixXT<double>. This would be done automatically by the compiler if pdist2 were not a template function:

MatrixXT<double> pdist2(const MatrixXT<double>&,const MatrixXT<double>&);

but since pdist2 is templated and that MatrixXT<double> is not a base class of Matrix<double, 3, 5>, in C++ the compiler is not allowed to automatically deduced the template parameter to instantiate pdist2. The solution for you is to generalize even more your function to take any MatrixBase<> as inputs:

template<typename D1,typename D2>
Matrix<typename D1::Scalar,D1::ColsAtCompileTime,D2::ColsAtCompileTime>
pdist2(const MatrixBase<D1>& _X,const MatrixBase<D2>& _Y);

Since X and Y are going to be used multiple times, and that they can now be arbitrary expressions (including a costly matrix product), you might want to let Eigen evaluate the arguments if needed, to this end, you can use Ref:

Ref<const typename D1::PlainObject> X(_X);
Ref<const typename D2::PlainObject> Y(_Y);

This way X will be evaluated if it cannot be represented as a pointer+stride to actual values.

Upvotes: 4

Related Questions