Compilation error when multiplying matrix of Jets by matrix of doubles

[Miłosz Makowski]

Hi,

I am using Eigen for matrix multiplication in Ceres solver.
Suppose I have two matrices, one consists of Jet<double, N>, whereas the other consists of double. When I try to calculate the product of the two matrices, I'm getting a compilation error. In short, the error boils down to an attempt of Jet-to-double conversion.

Here's an example of code that fails to compile:

-------------------------------------
typedef ceres::Jet<double, 2> J;

Eigen::MatrixX<J> A;
A.resize(2, 2);
A << (J)1., (J)1., (J)1., (J)1.;

Eigen::MatrixXd B;
B.resize(2, 2);
B << 1., 1., 1., 1.;

Eigen::MatrixX<J> C;
C.noalias() = A * B;
-------------------------------------

Note, that in this simple case, I could use Matrix2<J> and Matrix2d instead, and it would work fine. However, I need to use matrices of much higher size, so it's not possible to use compilation-time defined size.

What's more, if both of matrices consist of Jets, the example compiles fine. However, for performance reason, I don't want to convert the matrix of doubles to Jets just to make it compile.

It's also worth noting, that the multiplication works, if matrix of Jets is a vector. For now, the workaround I'm using is to just loop over the rows like that:

-------------------------------------
C.resize(A.rows(), B.cols());
for (int i = 0; i < A.rows(); i++) C.row(i).noalias() = A.row(i) * B;
-------------------------------------

It's far from perfect though.


Does anyone know how to solve that issue?
Thanks in advance


P.S.
A fragment of compilation output of the first example looks like this:

*************************************
include/Eigen/src/Core/util/BlasUtil.h:43:79: error: cannot convert 'const ceres::Jet<double, 2>' to 'double' without a conversion operator
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE To run(const From& x) { return To(x); }

include/Eigen/src/Core/GeneralProduct.h:246:66: note: in instantiation of member function 'Eigen::internal::get_factor<ceres::Jet<double, 2>, double>::run' requested here
RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
[...]
include/Eigen/src/Core/products/GeneralMatrixMatrix.h:445:7: note: in instantiation of function template specialization 'Eigen::internal::generic_product_impl<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::Matrix<double, -1, -1, 0>, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>>' requested here
scaleAndAddTo(dst, lhs, rhs, Scalar(1));
[...]
main.cpp:24:15: note: in instantiation of function template specialization 'Eigen::NoAlias<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::MatrixBase>::operator=<Eigen::Product<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::Matrix<double, -1, -1, 0>, 0>>' requested here
[build]   C.noalias() = A * B;
*************************************

[Sameer Agarwal]

+Sergiu who has been looking and fixing a number of these issues.

What version of Ceres Solver are you using?

Sameer

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/5cc62058-d9e3-4e62-b31b-1d094d0da44an%40googlegroups.com.

[Miłosz Makowski]

Ceres Solver: 2.0.0

Eigen: 3.4.0

- Milosz

[Sergiu Deitsch]

It looks like dynamic size matrices use a variant of GEMM which requires partially specializing Eigen::internal::gebp_traits and the corresponding packet math for mixed type operations. Eigen provides such a specialization currently only for std::complex which evidently does not generalize to Jet.

Given the necessary interfaces live in an internal namespace which is not intended for general use and can change anytime, I don't believe it is currently feasible to provide a workaround in Ceres.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/f7a50673-b7c3-4d92-945e-d201e6e000den%40googlegroups.com.

[Sameer Agarwal]

Milosz, this may be worth raising with the Eigen folks.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CALjVmNco%2BiH_%2BTQhH7fdXRdsoxTA0hSMSkAOsjs_kWmxf1-K-Q%40mail.gmail.com.

[Miłosz]

Sure,

Thank you for your help

- Milosz

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CABqdRUAy0v2LnWqO-APCf7hXn6a8-eswjognER5pD8DGruQ-ig%40mail.gmail.com.

[Alex Stewart]

For reference, this is the original CL from Chris for this kind of case that followed the Eigen documentation here, which does address some use-cases, but not all.  We also provide (AFAIK a correct) specialisation of NumTraits for Jet types as Eigen discuss in their docs for custom scalar types here.  As Sergiu says, the traits it appears we might have to specialise seem to be in the internal namespace, implying that we shouldn't do so.  However, if you get a response from the Eigen devs about what we should be doing to make it work we would be very happy to update Ceres.  From some initial investigation though, it looks like we're not the first people to encounter this, and I haven't found any documented fixes besides explicit casting.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/8afa10b2-56e8-6d3b-fd46-49a44fb50e7d%40mxlabs.ai.

[Miłosz Makowski]

Thank you for your help and suggestions. I posted the question on Eigen's issue list: https://gitlab.com/libeigen/eigen/-/issues/2434

- Miłosz