Rasmus,
I expected that ceres-solver traits do the job
#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)
And they do, but only for direct multiplication of Jet matrix and double
vector. It seems that it is not enough for the method solve.
On Thu, Jun 4, 2020 at 8:50 PM Rasmus Munk Larsen <[email protected]>
wrote:
> Oleg,
>
> Basically my response is that A.jacobiSVD().solve( b.cast<T>()) is what
> you should be doing - we made a conscious decision to disallow code to
> rely on implicit casting. Just speculation: Perhaps add your own
> specialization of BinOp(T, double) that does the right thing and is faster?
>
> On Wed, Jun 3, 2020 at 10:35 PM Oleg Shirokobrod <
> [email protected]> wrote:
>
>> 1. I would like to have autodiff ability, so I cannot use double for both
>> A and b. If I cast b: A.jacobiSVD().solve( b.cast<T>()) everything works
>> fine, but BinOp(T,T) is more expensive than BinOp(T, double). I would like
>> to keep b as a vector of doubles.
>> 2. T=Jet is ceres solver autodiff implementation type. There is a trait
>> definition for Jet binary operations for type deduction such that
>> type(Jet*double) = Jet and so on. It works when I do direct multiplication
>> VS^-1U^T
>> * b. It works similar to complex scalar matrices and double rhs and there
>> is the same problem for complex scalar cases.
>> 3. I think that the mixed type deduction rule should give the same type
>> for VS^-1U^T * b and for A.jcobianSVD().solve(b); where A = USV^T
>> because both use the same algorithm.
>> 4. Unless there are serious reasons, deduction rules should be similar to
>> scalar type equations. complex<double> A; double b; x = A^-1 * b; type(x) =
>> complex<double>.
>>
>> On Wed, Jun 3, 2020 at 11:16 PM Rasmus Munk Larsen <[email protected]>
>> wrote:
>>
>>> 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
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
