Combining Eigen and CppAD

Question

I want to use automatic differentiation mechanism provided by CppAD inside Eigen linear algebra. An example type is Eigen::Matrix< CppAD::AD,-1,-1>. As CppAD::AD is a custom numeric type the NumTraits for this type have to be provided. CppAD provides those in the file cppad/example/cppad_eigen.hpp. This makes the following minimal example compile:

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

int main() {
   typedef double Scalar;
   typedef CppAD::AD<Scalar> AD;

   // independent variable vector
   Eigen::Matrix<AD,Eigen::Dynamic,1> x(4);
   CppAD::Independent(x);

   // dependent variable vector 
   Eigen::Matrix<AD,Eigen::Dynamic,1> y(4);

   Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> m(4,4);
   m.setIdentity();

   y = 1. * x;
   // THIS DOES NOT WORK
   // y = m * x;

   CppAD::ADFun<Scalar> fun(x, y);
}

As soon as some more complex expression is used, e.g. the mentioned

y = m * x;

the code fails to compile:

PATH/Eigen/latest/include/Eigen/src/Core/Product.h:29:116: error: 
      no type named 'ReturnType' in 'Eigen::ScalarBinaryOpTraits<double, CppAD::AD<double>,
      Eigen::internal::scalar_product_op<double, CppAD::AD<double> > >'
  ...typename ScalarBinaryOpTraits<typename traits<LhsCleaned>::Scalar, typename traits<RhsCleaned>::Scalar>::ReturnType...

If i manually cast the double Matrix to AD, it works. However this is not a solution because the type promotion is virtually used everywhere in Eigen.

It looks to me as if the NumTraits provided by CppAD are not sufficient for this case. This is supported by a followup error message:

PATH/Eigen/latest/include/Eigen/src/Core/Product.h:155:5: error: 
      no type named 'CoeffReturnType' in
      'Eigen::internal::dense_product_base<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
      Eigen::Matrix<CppAD::AD<double>, -1, 1, 0, -1, 1>, 0, 7>'
    EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
    ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Other use cases lead to error messages like:

PATH/Eigen/src/Core/functors/Binary
Functors.h:78:92: error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<dou
ble, CppAD::AD<double>, Eigen::internal::scalar_product_op<double, CppAD::AD<double> > >’

Can anybody point me in the correct direction? It is possible that the NumTraits are for old Eigen Versions. I'm using 3.3.2 and CppAD from the current master branch.

Plimsoll

Why you don't want m to have type Eigen::Matrix<AD,Eigen::Dynamic,Eigen::Dynamic>? – Tympanites

Answer 1

If you want to multiply a Matrix<CppAD::AD<double>, ...> by a Matrix<double, ...> you also need to specialize the corresponding ScalarBinaryOpTraits:

namespace Eigen {
template<typename X, typename BinOp>
struct ScalarBinaryOpTraits<CppAD::AD<X>,X,BinOp>
{
  typedef CppAD::AD<X> ReturnType;
};

template<typename X, typename BinOp>
struct ScalarBinaryOpTraits<X,CppAD::AD<X>,BinOp>
{
  typedef CppAD::AD<X> ReturnType;
};
} // namespace Eigen

This requires that CppAD::AD<X>() * X() is implemented.

Alternatively, you need to cast your matrix m to AD:

// should work:
y = m.cast<CppAD::AD<double> >() * x;

Suspensory