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