Created attachment 383 [details]
There is a flaw in the sparse version of lmqrsolv in unsupported's LevenbergMarquardt module that makes the optimization assert with
Eigen/src/SparseCore/TriangularSolver.h:78: static void Eigen::internal::sparse_solve_triangular_selector<Lhs, Rhs, Mode, 2, 1>::run(const Lhs&, Rhs&) [with Lhs = Eigen::Block<Eigen::SparseMatrix<double, 1>, -1, -1, false>; Rhs = Eigen::Matrix<double, -1, 1>; int Mode = 2]: Assertion `it && it.index()==i' failed.
in certain cases. lmqrsolv iterates and updates only the non-zero entries in R, which does not suffice. Attached patch provides a fix, though in a trivial and less than optimal way.
Unfortunately, I do not have the time to go through the theory to see, which zero entries must not be ignored. Please let me know if I can assist in testing.
Just some background information:
The original documentation in MINPACK gives some more insight into lmqrsolv, for example at http://www.netlib.no/netlib/sminpack/qrsolv.f. The task is to solve a system of linear equations A*x=b with A=(B*P|D*P)^T composed of
* a matrix B, for which we computed a QR factorization before
* and a (permuted) diagonal matrix D.
The current implementation updates the QR factorization of B by Givens rotation to a QR factorization of entire A.
Does this update by Givens rotations make sense for a sparse matrix? Or would it be better to compute a fresh QR decomposition for A with SparseQR?
-- GitLab Migration Automatic Message --
This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.
You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/657.