Rasmus, I wiil have a look at COD. Brad, I did not try CppAD.I am working in given framework: ceres nonlinear least squares solver + ceres autodiff + Eigen decomposition modules SVD or QR. The problem is not just on autodiff side. The problem is that Eigen decomposition modul does not work properly with autodiff type variable.
Thank you everybody for advice. On Mon, Jun 1, 2020 at 8:41 PM Rasmus Munk Larsen <[email protected]> wrote: > > > On Mon, Jun 1, 2020 at 10:33 AM Patrik Huber <[email protected]> > wrote: > >> Hi Rasmus, >> >> This is slightly off-topic to this thread here, but it would be great if >> you added your COD to the list/table of decompositions in Eigen: >> https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html >> >> First, it would make it easier for people to find, and second, it would >> also help a lot to see on that page how the algorithm compares to the >> others, to be able to choose it appropriately. >> > > Good point. Will do. > > >> >> >> Unrelated: @All/Maintainers: It seems like lots (all) of the images on >> the documentation website are broken? At least for me. E.g.: >> >> [image: image.png] >> >> >> Best wishes, >> Patrik >> >> On Mon, 1 Jun 2020 at 17:59, Rasmus Munk Larsen <[email protected]> >> wrote: >> >>> 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 >>>>> >>>>>
