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