Try to compile your code in debug mode with the type assertions on. On Wed, Jun 3, 2020 at 1:14 PM Rasmus Munk Larsen <[email protected]> wrote:
> Are you saying that you compute the decomposition in one type and solve > with a RHS of a different type? Why do you say that VS^-1U^T * b should be > Matrix<T>? That makes an assumption about type coercion rules. In fact, you > cannot generally mix types in Eigen expressions without explicit casting, > and U.adjoint() * b should fail if the types are different. > > On Wed, Jun 3, 2020 at 11:33 AM Oleg Shirokobrod < > [email protected]> wrote: > >> Rasmuss, I do not quite understand this issue. Decomposition solve should >> propagate scalar type of a matrix but not scalar type of its argument. >> Example: >> template <typename T> Matrix<T> A; >> VectorXd b; >> A.jcobiSVD().solve(b) should be of type Matrix<T> but it is not. Type of >> result is Matrix<double>. If we make SVD decomposition of A = USV^T and >> express result as VS^-1U^T * b, than result will be of type Matrix<T>. >> Which is correct and differs from result of solve which uses the same >> algorithm but more complex result’s type deduction. This is the problem. >> >> On Wed 3. Jun 2020 at 20.19, Rasmus Munk Larsen <[email protected]> >> wrote: >> >>> OK, I opened https://gitlab.com/libeigen/eigen/-/issues/1909 for this. >>> >>> On Tue, Jun 2, 2020 at 11:06 PM Oleg Shirokobrod < >>> [email protected]> wrote: >>> >>>> Yes. At the time of computing only 1d observation (VectorXd) is known. >>>> >>>> On Tue, Jun 2, 2020 at 9:42 PM Rasmus Munk Larsen <[email protected]> >>>> wrote: >>>> >>>>> OK, so b is declared as VectorXf or some other type with >>>>> ColsAtCompileTime=1? >>>>> >>>>> On Tue, Jun 2, 2020 at 11:27 AM Oleg Shirokobrod < >>>>> [email protected]> wrote: >>>>> >>>>>> >>>>>> Yes, b is measured spectrum that is 1d array. I have to get x for 1d >>>>>> array at a time. I fit sum of peak models to 1d rhs. 1d array of peak >>>>>> model >>>>>> values is one column of matrix A. >>>>>> >>>>>> On Tue 2. Jun 2020 at 20.06, Rasmus Munk Larsen <[email protected]> >>>>>> wrote: >>>>>> >>>>>>> Why do you say that? You could be solving for multiple >>>>>>> right-hand sides. Is b know to have 1 column at compile time? >>>>>>> >>>>>>> On Tue, Jun 2, 2020 at 1:31 AM Oleg Shirokobrod < >>>>>>> [email protected]> wrote: >>>>>>> >>>>>>>> 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 >>>>>>>>>>>>>> >>>>>>>>>>>>>>
