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. 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 >>> >>>
