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