Re: [eigen] Using Eigen decompositions with ceres autodiff Jet data type

On Mon, Jun 1, 2020 at 8:44 AM Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx> wrote:

Hi Sameer,

Thank you for quick reply and advices.
1. I implemented using both decompositions SVD and QR. QR decomposition has the same problem as SVD for Jet but explicit _expression_ for solution is much more complex.
2. I implemented not only CostFunctor with operator()() for autodiff but also CostFunction with Evaluate method. In this method I calculate residuals and Jacobian using method similar to your proposal. This is main idea of variable projection method. Using pseudo-inverse solution from SVD or QR decomposition eliminate linear parameters (project linear variable out) and then calculate Jcobian using closed _expression_ for residuals. But calculating Jacobian is not trivial 1),2) and time consuming. I am going to profile both approaches. Physical meaning of my simplified model is linear combination of gaussian. Linear parameters are their amplitudes, nonlinear are peak positions and widths. When we use unconstraint linear problem, there is closed _expression_ for residuals which we can differentiate. If we put constraints on linear parameters, e.g. nonnegativity of them, we do not have closed expressions for residuals, we need differentiate some algorithm, for example Quadratic Programming with constrains. There are articles which offer calculating Jcobian of close task near optimum without constraints. I think it is beneficial in this case calculate Jcobian of algorithm solving exact task.
3. BW. I compiled ceres solver, glog and gflags with VS2019. I do not remember if I tweaked a bit CMake VS project. I can provide solutions for these projects. But I added derivative of erf and erfc functions, which I use and which are part of C++11 standard. Besides I have to enclose in to pragmas
#ifdef _MSC_VER
#pragma optimize( "", off )
#endif // _MSC_VER

#ifdef _MSC_VER
#pragma optimize( "", on )
#endif // _MSC_VER

class GOOGLE_GLOG_DLL_DECL CheckOpMessageBuilder and class GOOGLE_GLOG_DLL_DECL LogMessageVoidify in glog::logging.h without this VS compiler in release mode optimizes out some member functions. After that linker gives unresolved external error message.

References
1) G. H. Golub and V. Pereyra, The differentiation of pseudo-inverses and nonlinear least squares whose variables separable. SIAM J. Numer. Anal. Vol. 10. Num. 2.1973.
2). O Leary. Rust. Variable Projection for Nonlinear Least Squares Problems.

Again thank you very much. I will think of your advice.

Oleg

On Mon, Jun 1, 2020 at 4:57 PM Sameer Agarwal <sameeragarwal@xxxxxxxxxx> wrote:
Oleg,
Two ideas:

1. You may have an easier time using QR factorization instead of SVD to solve your least squares problem.
2. But you can do better, instead of trying to solve linear least squares problem involving a matrix of Jets, you are better off, solving the linear least squares problem on the scalars, and then using the implicit function theorem to compute the derivative w.r.t the parameters and then applying the chain rule.

i.e., start with min |A x = b|

the solution satisfies the equation

A'A x - A'b = 0.

solve this equation to get the optimal value of x, and then compute the jacobian of this equation w.r.t A, b and x. and apply the implicit theorem.

Sameer

On Mon, Jun 1, 2020 at 4:46 AM Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx> wrote:
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.singularValues().head(l_rank).asDiagonal().inverse()
         * 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