Hi Rasmus, I have just tested COD decomposition in Eigen library. It arises the same problem. This is defect of Eigen decomposition module type reduction of result of solve method. If template <typename T> Matrix<T, Dynamic, Dynamic> A; and ArraXd b;, than x = A.solve(b) should be of type <typename T> Matrix<T, Dynamic, 1.>.
I like the idea to use COD as an alternative to QR or SVD and I added this option to my code. On Tue, Jun 2, 2020 at 10:36 AM Oleg Shirokobrod <[email protected]> wrote: > 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 >>>>>> >>>>>>
