Combining Eigen and CppAD
I want to use the automatic differentiation mechanism provided by CppAD inside its own linear algebra. An example of the type Eigen :: Matrix <CppAD :: AD, -1, -1>. Since CppAD :: AD is a custom numeric type, NumTraits must be supplied for this type. CppAD provides the ones found in the cppad / example / cppad_eigen.hpp file . This does it after compiling the minimal example:
#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 a more complex expression is used like the one mentioned
y = m * x;
the code won't 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 drop the double matrix in AD it works. However, this is not a solution because the type promotion is actually used throughout Eigen.
It seems to me that the NumTraits provided by CppAD is not sufficient for this case. This is confirmed by the 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 anyone point me in the right direction? Possibly NumTraits for old native versions. I am using 3.3.2 and CppAD from the current master branch.
source to share
If you want to propagate Matrix<CppAD::AD<double>, ...>
by Matrix<double, ...>
, you will also need to specialize the appropriate one 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 CppAD::AD<X>() * X()
.
Alternatively, you need to cast the matrix down m
to AD
:
// should work:
y = m.cast<CppAD::AD<double> >() * x;
source to share