组合Eigen和CppAD

Question

我想在本征线性代数中使用CppAD提供的自动微分机制。一个示例类型是Eigen::Matrix< CppAD::AD，-1，-1>。因为CppAD::AD是一个自定义的数字类型，所以必须提供此类型的NumTraits。CppAD在文件cppad/example/cppad_eigen.hpp中提供了这些。这将使以下最小示例进行编译：

#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);
}

一旦使用了一些更复杂的表达式，例如前面提到的

y = m * x;

代码编译失败：

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...

如果我手动将double Matrix转换为AD，它可以工作。然而，这不是一个解决方案，因为类型提升实际上在Eigen中随处可见。

在我看来，CppAD提供的NumTraits对于这种情况是不够的。以下错误消息支持这一点：

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)
    ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

其他用例会导致错误消息，如：

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> > >’

有谁能告诉我正确的方向吗？有可能NumTraits是针对旧的特征版本的。我使用的是当前主分支中的3.3.2和CppAD。

Answer 1

如果要将Matrix<CppAD::AD<double>, ...>与Matrix<double, ...>相乘，还需要专门化相应的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

这需要实现CppAD::AD<X>() * X()。

或者，您需要将矩阵m转换为AD

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