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