row times column is not a scalar, with custom scalar type

user_9908011 0
Viewed 134

I am trying to use the CppAD scalar type with Eigen.

The following fails to compile:

#include <Eigen/Dense>
#include <cppad/cppad.hpp>
#include <cppad/example/cppad_eigen.hpp>

int main()
{
    using Scalar = CppAD::AD<double>;
    //using Scalar = double;

    Eigen::Matrix<Scalar, 1,4> row;
    Eigen::Matrix<Scalar, 4,1> col;
    Scalar scalar = 5;
    
    Scalar res2 = row * col + scalar; //fails

    return 0;
}

The error is basically that it cannot add a scalar to the result of the multiplication. However, the result of the multiplication is a scalar itself, so it should not be a problem. Indeed, when using double as the Scalar type, there are no issues.

Here is the compiler error:

cppad-eigen-problem.cpp:14:29: error: no match for ‘operator+’ (operand types are ‘const Eigen::Product<Eigen::Matrix<CppAD::AD<double>, 1, 4>, Eigen::Matrix<CppAD::AD<double>, 4, 1, 0, 4, 1>, 0>’ and ‘Scalar’ {aka ‘CppAD::AD<double>’})
   14 |     Scalar res2 = row * col + scalar; //fails
      |                   ~~~~~~~~~ ^ ~~~~~~
      |                       |       |
      |                       |       Scalar {aka CppAD::AD<double>}
      |                       const Eigen::Product<Eigen::Matrix<CppAD::AD<double>, 1, 4>, Eigen::Matrix<CppAD::AD<double>, 4, 1, 0, 4, 1>, 0>

There is an issue on the CppAD project, but I am not sure where the problem is:

  1. If I use another trivial custom scalar, I cannot reproduce the error...

  2. On the other hand, the CppAD Eigen traits for their scalar types looks ok to me.

Versions: Eigen 3.3.7, latest master of CppAD, g++ 9.3.0

Any clue?

Thanks

eigen linear-algebra c++
1 Answers

While it is mathematically true that a (1x1) matrix is a scalar, C++ is a different story: The return type of operator* in row * col is not a scalar and neither a (1x1) matrix but a product expression. This product expression is implicitly convertible to the scalar type of its operands, in this case CppAD::AD<double>.

The reason for the error you see is that CppAD::AD<double> is a template and its operators are thus function templates, e.g. something like

template<typename T>
CppAD::AD<T> operator+(const CppAD::AD<T>& lhs, const CppAD::AD<T>& rhs);

The problem now is that the above templated operator+ cannot be called since there is no implicit conversion in template arguments and exactly this would be needed to cast the product expression to CppAD::AD<double>.

That's also the reason why avoiding the product works (e.g. by using Eigen's .dot function).

Nevertheless you can solve this problem by e.g. defining an appropriate operator+, something like (not tested)

template<typename Derived>
CppAD::AD<double>  operator+(const MatrixBase<Derived>& lhs, const CppAD::AD<double> & rhs) {
  return lhs.derived().coeff(0, 0) + rhs;
                     //^^^^^^^^^^^ No need for an implicit conversion
}

or you could use Eigen's plugin mechanism and add an operator to MatrixBase.

user_64323220
Related