[eigen] Using Eigen decompositions with ceres autodiff Jet data type.
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: [eigen] Using Eigen decompositions with ceres autodiff Jet data type.
• From: Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx>
• Date: Mon, 1 Jun 2020 14:45:13 +0300
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=mime-version:from:date:message-id:subject:to; bh=PPcC9M15BLQrjZXmWBbqqBYnZTZPzUeseP5QVudSBck=; b=MmOzwWBfSeU5GxH4pDCfp560rapL97yz1ifV61qEKpY549NoC4kWKYzWLHorXO19kO 4tJYXq9pM3AlGCa+k+5yjswhA4P3HAD/mpWxQREomc64/vAbdyRE2oBB1JWzhqV/g7yL w8ibwXkyi9+J+dU4QXKf5bxtU8hSBe3LN4/wyqKy8i/VVxgNjOfWMoQA3E4LQnw1MqQc xDqlqBbJjiYyjDfL44N1+jDsxWlNPvgDZgdhMlhNA0gH48gsIo68/vStVP5GHi6lxPbw CIlbgsiObUnuQi1f5/lbKRTIexE7bh99fLStrfLSuu6oFU5qET6bQAR8uY5IC/p9+yjK rDPQ==

Hi list, I am using Eigen 3.3.7 release with ceres solver 1.14.0 with autodiff Jet data type and I have some problems. I need to solve linear least square subproblem within variable projection algorithm, namely I need to solve LLS equation
A(p)*x = b
Where matrix A(p) depends on nonlinear parameters p:
x(p) = pseudo-inverse(A(p))*b;
x(p) will be optimized in nonlinear least squares fitting, so I need Jcobian. Rhs b is measured vector of doubles, e.g. VectorXd.. In order to use ceres's autodiff p must be of Jet type. Ceres provides corresponding traits for binary operations

#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
// Specifying the return type of binary operations between Jets and scalar types
// allows you to perform matrix/array operations with Eigen matrices and arrays
// such as addition, subtraction, multiplication, and division where one Eigen
// matrix/array is of type Jet and the other is a scalar type. This improves
// performance by using the optimized scalar-to-Jet binary operations but
// is only available on Eigen versions >= 3.3
template <typename BinaryOp, typename T, int N>
struct ScalarBinaryOpTraits<ceres::Jet<T, N>, T, BinaryOp> {
typedef ceres::Jet<T, N> ReturnType;
};
template <typename BinaryOp, typename T, int N>
struct ScalarBinaryOpTraits<T, ceres::Jet<T, N>, BinaryOp> {
typedef ceres::Jet<T, N> ReturnType;
};
#endif  // EIGEN_VERSION_AT_LEAST(3, 3, 0)

There two problems.
1. Small problem. In a function "RealScalar threshold() const" in SCDbase.h I have to replace "return m_usePrescribedThreshold ? m_prescribedThreshold
: diagSize* NumTraits<Scalar>::epsilon();" with "return m_usePrescribedThreshold ? m_prescribedThreshold
: Scalar(diagSize)* NumTraits<Scalar>::epsilon();"
This fix is similar Gael's fix of Bug 1403
2. It is less trivial. I expect that x(p) = pseudo-inverse(A(p))*b; is vector of Jet. And it is actually true for e.g SVD decompoazition
x(p) = VSU^T * b.
But if I use
JcobySVD<Matrix<Jet<double, 2>, Dynamic, Dynamic>> svd(A);
x(p) = svd.solve(b),
I got error message.
Here code for reproducing the error

// test_svd_jet.cpp
#include <ceres/jet.h>
using ceres::Jet;

int test_svd_jet()
{
typedef Matrix<Jet<double, 2>, Dynamic, Dynamic> Mat;
typedef Matrix<Jet<double, 2>, Dynamic, 1> Vec;
Mat A = MatrixXd::Random(3, 2).cast <Jet<double, 2>>();
VectorXd b = VectorXd::Random(3);
JacobiSVD<Mat> svd(A, ComputeThinU | ComputeThinV);
int l_rank = svd.rank();
Vec c = svd.matrixV().leftCols(l_rank)
* svd.matrixU().leftCols(l_rank).adjoint() * b; // *
Vec c1 = svd.solve(b.cast<Jet<double, 2>>()); // **
Vec c2 = svd.solve(b); // ***
return 0;
}
// End test_svd_jet.cpp

// * and // ** work fine an give the same results. // *** fails with VS 2019 error message
Eigen\src\Core\functors\AssignmentFunctors.h(24,1):
error C2679: binary '=': no operator found which takes
a right-hand operand of type 'const SrcScalar'
(or there is no acceptable conversion)
The error points to line //***. I thing that solution is of type VectorXd instead of Vec and there is problem with assignment of double to Jet. Derivatives are lost either. It should work similar to complex type. If A is complex matrix and b is real vector, x must be complex. There is something wrong with Type deduction in SVD or QR decomposition.

Do you have any idea of how to fix it.

Best regards,

Oleg Shirokobrod

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/