Hi Oleg and Sameer, A faster option than SVD, but more robust than QR (since it also handles the under-determined case) is the complete orthogonal decomposition that I implemented in Eigen a few years ago.
https://eigen.tuxfamily.org/dox/classEigen_1_1CompleteOrthogonalDecomposition.html (Looks like the docstring is broken - oops!) It appears to also be available in the 3.3 branch: https://gitlab.com/libeigen/eigen/-/blob/3.3/Eigen/src/QR/CompleteOrthogonalDecomposition.h Rasmus On Mon, Jun 1, 2020 at 6:57 AM Sameer Agarwal <[email protected]> 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 <https://en.wikipedia.org/wiki/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 < > [email protected]> 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 >> <http://eigen.tuxfamily.org/bz/show_bug.cgi?id=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 >> >>
