Products of MatrixXd and MatrixXdual not working

[rath3t]

In my code the products of type double and type dual inside Eigen Vectors and Matrices seems to be not working. Do i overlook something? My MWE looks like the following

#include <iostream>
#include <Eigen/Core>
#include <autodiff/forward.hpp>
#include <autodiff/forward/eigen.hpp>

using namespace Eigen;
using namespace autodiff;
using namespace std;

dual testEnergy(const Vector3dual& DDual) {

    MatrixX3dual   C(2,3);
    Matrix2Xd A(2,4);
    MatrixX3dual B(4 ,3);
    A << 1, 2, 3, 4, 5, 6, 7, 8;

    B.setIdentity(4,3);
    B = DDual.norm()*B;
    C = A*B; //Not working
//    C = A.cast<dual>() * B; //Working but useless temporary constructed of type dual type?
return C.norm();
}

int main() {
    Vector3dual DDual;
    DDual<< 1,5,13;
    dual dualenergy;
    VectorXd g = gradient(testEnergy, wrt(DDual), at(DDual), dualenergy);
    cout<<g<<endl;
    return 0;
}

which produces the following compilation error:

Scanning dependencies of target Test_AutoDiff
[ 50%] Building CXX object CMakeFiles/Test_AutoDiff.dir/main.cpp.obj
In file included from PATH/src/eigen3/Eigen/Core:313:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h: In instantiation of 'void Eigen::internal::gebp_traits<_LhsScalar, _RhsScalar, _ConjLhs, _ConjRhs, Arch, _PacketSize>::madd(const LhsPacketType&, const RhsPacketType&, AccPacketType&, RhsPacketType&, const LaneIdType&) const [with LhsPacketType = autodiff::forward::Dual<double, double>; RhsPacketType = double; AccPacketType = autodiff::forward::Dual<double, double>; LaneIdType = Eigen::internal::FixedInt<0>; _LhsScalar = autodiff::forward::Dual<double, double>; _RhsScalar = double; bool _ConjLhs = false; bool _ConjRhs = false; int Arch = 0; int _PacketSize = 0]':
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:2094:15:   required from 'void Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(const DataMapper&, const LhsScalar*, const RhsScalar*, Index, Index, Index, Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::ResScalar, Index, Index, Index, Index) [with LhsScalar = double; RhsScalar = autodiff::forward::Dual<double, double>; Index = int; DataMapper = Eigen::internal::blas_data_mapper<autodiff::forward::Dual<double, double>, int, 0, 0, 1>; int mr = 1; int nr = 4; bool ConjugateLhs = false; bool ConjugateRhs = false; Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::ResScalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:198:15:   required from 'static void Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::run(Index, Index, Index, const LhsScalar*, Index, const RhsScalar*, Index, Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar*, Index, Index, Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar, Eigen::internal::level3_blocking<LhsScalar, RhsScalar>&, Eigen::internal::GemmParallelInfo<Index>*) [with Index = int; LhsScalar = double; int LhsStorageOrder = 0; bool ConjugateLhs = false; RhsScalar = autodiff::forward::Dual<double, double>; int RhsStorageOrder = 0; bool ConjugateRhs = false; int ResInnerStride = 1; Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:230:14:   required from 'void Eigen::internal::gemm_functor<Scalar, Index, Gemm, Lhs, Rhs, Dest, BlockingType>::operator()(Index, Index, Index, Index, Eigen::internal::GemmParallelInfo<Index>*) const [with Scalar = autodiff::forward::Dual<double, double>; Index = int; Gemm = Eigen::internal::general_matrix_matrix_product<int, double, 0, false, autodiff::forward::Dual<double, double>, 0, false, 0, 1>; Lhs = Eigen::Matrix<double, 2, -1>; Rhs = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Dest = Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>; BlockingType = Eigen::internal::gemm_blocking_space<0, double, autodiff::forward::Dual<double, double>, 2, 3, -1, 1, false>]'
PATH/src/eigen3/Eigen/src/Core/products/Parallelizer.h:114:7:   required from 'void Eigen::internal::parallelize_gemm(const Functor&, Index, Index, Index, bool) [with bool Condition = false; Functor = Eigen::internal::gemm_functor<autodiff::forward::Dual<double, double>, int, Eigen::internal::general_matrix_matrix_product<int, double, 0, false, autodiff::forward::Dual<double, double>, 0, false, 0, 1>, Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>, Eigen::internal::gemm_blocking_space<0, double, autodiff::forward::Dual<double, double>, 2, 3, -1, 1, false> >; Index = int]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:509:9:   required from 'static void Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo(Dest&, const Lhs&, const Rhs&, const Scalar&) [with Dest = Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>; Lhs = Eigen::Matrix<double, 2, -1>; Rhs = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::Scalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:445:20:   [ skipping 5 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]
PATH/src/eigen3/Eigen/src/Core/Matrix.h:332:31:   required from 'Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Matrix(const T&) [with T = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; _Scalar = autodiff::forward::Dual<double, double>; int _Rows = 2; int _Cols = 3; int _Options = 0; int _MaxRows = 2; int _MaxCols = 3]'
PATH/src/eigen3/Eigen/src/Core/AssignEvaluator.h:822:41:   required from 'void Eigen::internal::call_assignment(Dst&, const Src&, const Func&, typename Eigen::internal::enable_if<Eigen::internal::evaluator_assume_aliasing<Src>::value, void*>::type) [with Dst = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Src = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; Func = Eigen::internal::assign_op<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >; typename Eigen::internal::enable_if<Eigen::internal::evaluator_assume_aliasing<Src>::value, void*>::type = void*]'
PATH/src/eigen3/Eigen/src/Core/AssignEvaluator.h:808:18:   required from 'void Eigen::internal::call_assignment(Dst&, const Src&) [with Dst = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Src = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>]'
PATH/src/eigen3/Eigen/src/Core/PlainObjectBase.h:779:32:   required from 'Derived& Eigen::PlainObjectBase<Derived>::_set(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; Derived = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>]'
PATH/src/eigen3/Eigen/src/Core/Matrix.h:225:24:   required from 'Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::operator=(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; _Scalar = autodiff::forward::Dual<double, double>; int _Rows = -1; int _Cols = 3; int _Options = 0; int _MaxRows = -1; int _MaxCols = 3]'
C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:20:11:   required from here
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:18: error: cannot convert 'Eigen::internal::conj_helper<autodiff::forward::Dual<double, double>, double, false, false>::Scalar {aka autodiff::forward::Dual<double, double>}' to 'double' in assignment
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:44: error: no matching function for call to 'padd(autodiff::forward::Dual<double, double>&, double&)'
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
                                        ~~~~^~~~~~~
In file included from PATH/src/eigen3/Eigen/Core:160:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/GenericPacketMath.h:163:1: note: candidate: template<class Packet> Packet Eigen::internal::padd(const Packet&, const Packet&)
 padd(const Packet& a, const Packet& b) { return a+b; }
 ^~~~
PATH/src/eigen3/Eigen/src/Core/GenericPacketMath.h:163:1: note:   template argument deduction/substitution failed:
In file included from PATH/src/eigen3/Eigen/Core:313:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:44: note:   deduced conflicting types for parameter 'const Packet' ('autodiff::forward::Dual<double, double>' and 'double')
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
                                        ~~~~^~~~~~~
CMakeFiles\Test_AutoDiff.dir\build.make:62: recipe for target 'CMakeFiles/Test_AutoDiff.dir/main.cpp.obj' failed
CMakeFiles\Makefile2:74: recipe for target 'CMakeFiles/Test_AutoDiff.dir/all' failed

[allanleal]

Looks like something similar to this:
https://stackoverflow.com/questions/49919005/cppad-with-eigen

where Eigen and CppAD are used together. I'll try to get this fixed soon (hopefully today, if the solution is there already! :)). Thanks for reporting, Alexander.

[allanleal]

Looks like this is a common (and non-resolved) issue when using custom types in Eigen. Here is another stackoverflow question (without a definite solution):
https://stackoverflow.com/questions/49456181/using-eigen-with-custom-complex-number-designed-for-automatic-differentiation

Even adding the following piece of code in autodiff/forward/eigen.hpp does not solve the issue:

template<typename BinOp>
struct ScalarBinaryOpTraits<autodiff::dual, double, BinOp>
{
    typedef autodiff::dual ReturnType;
};

template<typename BinOp>
struct ScalarBinaryOpTraits<double, autodiff::dual, BinOp>
{
    typedef autodiff::dual ReturnType;
};

Would you please help me finding a solution for this, by possibly asking to Eigen developers? They most likely know how to resolve this.

[rath3t]

I also searched a bit and found
https://stackoverflow.com/questions/41494288/mixing-scalar-types-in-eigen
where they proposed to use
C = A.lazyproduct(B);
which indeed works but since it does not use optimization and is slow as mentioned in
https://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#title70 this is no option.

Interestingly it also works for fixed sizes, e.g. Matrix3dual,Matrix3d.

I'm reading right now through
https://eigen.tuxfamily.org/bz/show_bug.cgi?id=279
which also seems interesting.
and the PR:
https://bitbucket.org/eigen/eigen/pull-requests/194/relax-mixing-type-constraints-for-binary/diff

I'm confident to find a solution since a lot of people in the automatic diff community wasrequesting such a feature, e.g. Adol-C, CodiPack etc.

[rath3t]

I searched for several AD tools which also use eigen but i didn't fsomething useful in cppAD,adol-c, CodiPack, Stan math and https://gitlab.com/tesch1/cppduals/issues/3 has the same problem. Therefore i think casting is right know the only option.

Since the Eigen people are aware of that problem and they fixed it in eigen 3.0 for Arrays i don't know how to think about this.

Any thoughts about this?

[allanleal]

My opinion is that we keep this issue open and stay on top of the developments of Eigen in the direction of permitting the multiplication of matrices with different scalar types. Until then, I hope you can find an adequate temporary solution to your use case. ;)

Let me know if you eventually find out something new about this.

[ludkinm]

FYI, this also affects arbitrary operations when mixing Eigen types

#include "Eigen/Core"
#include "autodiff/reverse.hpp"
#include "autodiff/reverse/eigen.hpp"

int main(void){
  Eigen::VectorXvar Xv(10);
  Eigen::VectorXd Y(10);
  autodiff::var z = (Xv -Y).normSquared();
  return 0;
}

Would it be possible to make var a class template?
Then have a var with overloads for simple operations (i.e. those in the matrix cookbook)

[allanleal]

Hi @ludkinm - I'll check a fix for this. This should work for var - I have mixed VectorXdual and VectorXd in the past without issues.

[ludkinm]

Thanks. I am using Eigen3.3.7 from the latest stable tarball and the master branch of autodiff if that helps.

[ludkinm]

@allanleal, you nudge that VectorXdual works has helped me identify whats missing:
the Eigen helper struct ScalarBinaryOpTraits is not defined for var.
I'll make a PR to fix this.

[allanleal]

Excellent you've found this. Yes, a PR would be great - thanks.