Can't "multiply" Eigen::Tensors

General
1

Hey folks,

I have the following code:

template <typename T>
T my_dot (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
    const Eigen::array<const Eigen::IndexPair<int>, 1> product_dims = {Eigen::IndexPair<int>(0,0)};
    const Eigen::Tensor<double, 1> a_tensor = Eigen::TensorMap<const Eigen::Tensor<const double, 1>> (a.data(), a.size());;
    const Eigen::Tensor<T, 1> b_tensor = Eigen::TensorMap<const Eigen::Tensor<const T, 1>> (b.data(), b.size());
    const Eigen::Tensor<T, 0> ret = a_tensor.contract(b_tensor, product_dims);
    return ret();
}

template <typename T>
T my_dot2 (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
    return (a*b).sum();
}

my_dot is templated so as to allow automatic differentiation via stan. So the substitution I’m most interested in is T=stan::math::var.

My issue is that my_dot2 works perfectly with autodiff types, but my_dot won’t compile.

/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:1555:15: error: no matching function for call to ‘Eigen::internal::gebp_traits<stan::math::var, double, false, false>::madd(Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SLhsPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SAccPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket&)’
               straits.madd(A0,B_0,C0,B_0);
               ^~~~~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:435:28: note: candidate: template<class LhsPacketType, class RhsPacketType, class AccPacketType> void Eigen::internal::gebp_traits<_LhsScalar, _RhsScalar, _ConjLhs, _ConjRhs>::madd(const LhsPacketType&, const RhsPacketType&, AccPacketType&, AccPacketType&) const [with LhsPacketType = LhsPacketType; RhsPacketType = RhsPacketType; AccPacketType = AccPacketType; _LhsScalar = stan::math::var; _RhsScalar = double; bool _ConjLhs = false; bool _ConjRhs = false]
   EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, AccPacketType& tmp) const
                            ^~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:435:28: note:   template argument deduction/substitution failed:
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:1555:15: note:   deduced conflicting types for parameter ‘AccPacketType’ (‘stan::math::var’ and ‘Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket {aka double}’)
               straits.madd(A0,B_0,C0,B_0);
               ^~~~~~~

Or if I change my_dot by swapping the order of arguments:

template <typename T>
T my_dot (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
    const Eigen::array<const Eigen::IndexPair<int>, 1> product_dims = {Eigen::IndexPair<int>(0,0)};
    const Eigen::Tensor<double, 1> a_tensor = Eigen::TensorMap<const Eigen::Tensor<const double, 1>> (a.data(), a.size());;
    const Eigen::Tensor<T, 1> b_tensor = Eigen::TensorMap<const Eigen::Tensor<const T, 1>> (b.data(), b.size());
    const Eigen::Tensor<T, 0> ret = b_tensor.contract(a_tensor, product_dims);
    return ret();
}

then I get

/usr/include/eigen3/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h:418:117: error: no matching function for call to ‘Eigen::internal::general_matrix_vector_product<long int, stan::math::var, Eigen::internal::TensorContractionInputMapper<stan::math::var, long int, 1, Eigen::TensorEvaluator<const Eigen::Tensor<stan::math::var, 1, 0, long int>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 1, true, false, 16>, 0, false, double, Eigen::internal::TensorContractionInputMapper<double, long int, 0, Eigen::TensorEvaluator<const Eigen::Tensor<double, 1>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 2, true, true, 16>, false, 0>::run(const Index&, const Index&, LhsMapper&, RhsMapper&, Eigen::TensorContractionEvaluatorBase<Eigen::TensorEvaluator<const Eigen::TensorContractionOp<const std::array<const Eigen::IndexPair<int>, 1>, const Eigen::Tensor<stan::math::var, 1, 0, long int>, const Eigen::Tensor<double, 1> >, Eigen::DefaultDevice> >::Scalar*&, const Index&, const Scalar&)’
     internal::general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,false,RhsScalar,RhsMapper,false>::run(
     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
         rows, cols, lhs, rhs,
         ~~~~~~~~~~~~~~~~~~~~~                                                                                        
         buffer, resIncr, alpha);
         ~~~~~~~~~~~~~~~~~~~~~~~                                                                                      
In file included from /usr/include/eigen3/Eigen/Core:473:0,
                 from /usr/include/eigen3/Eigen/Dense:1,
                 from [redacted]
                 from [redacted]
/usr/include/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h:88:24: note: candidate: static void Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::run(Index, Index, const LhsMapper&, const RhsMapper&, Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::ResScalar*, Index, RhsScalar) [with Index = long int; LhsScalar = stan::math::var; LhsMapper = Eigen::internal::TensorContractionInputMapper<stan::math::var, long int, 1, Eigen::TensorEvaluator<const Eigen::Tensor<stan::math::var, 1, 0, long int>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 1, true, false, 16>; bool ConjugateLhs = false; RhsScalar = double; RhsMapper = Eigen::internal::TensorContractionInputMapper<double, long int, 0, Eigen::TensorEvaluator<const Eigen::Tensor<double, 1>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 2, true, true, 16>; bool ConjugateRhs = false; int Version = 0; Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::ResScalar = stan::math::var]
 EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run(
                        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h:88:24: note:   no known conversion for argument 7 from ‘const Scalar {aka const stan::math::var}’ to ‘double’

I feel like this is probably an issue with Eigen as opposed to stan, but I haven’t heard back from them: https://forum.kde.org/viewtopic.php?f=74&t=141239&sid=f2808fe4aa7300c1e35456338fd801d0

Do you guys happen to know how to work around this issue? For now, I’m doing for loops…

2

const Eigen::Tensor<T, 0> ret = a_tensor.contract(b_tensor, product_dims);

a_tensor has scalar type double and b_tensor has scalar type stan::math::var.

I don’t think it knows how to contract tensors of these things together, even though we know all the scalar ops are in place to do this.

Can you try casting a_tensor to an Eigen::Tensor<T> manually and see what happens? Maybe there’s an easy way to do it with Issue casting C++ Eigen::Matrix types via templates - Stack Overflow .

Hope that helps

3

Exactly. To make these things efficient, we’d need to provide specializations for our autodiff types. Otherwise, you have to promote.