Levenberg-Marquardt

[stum...@gmail.com]

Have you considered implementing a Levenberg-Marquardt algorithm (like that implemented in MATLAB's lsqnonlin) for solving a nonlinear least squares problem in Manopt? Levenberg-Marquardt is more efficient than a trust-region algorithm because it only requires gradients.

[BM]

Hello, could we do that by modifying the Hessian appropriately, e.g., with grad*grad', and then calling the trust region solver?

[Nicolas Boumal]

Hello,

That would be a nice thing to implement indeed.

We do not have an implementation in Manopt, mostly because the primary purpose of the toolbox is to provide generic optimization algorithms (and a collection of geometries and generically useful tools), whereas Levenberg-Marquardt is specific to nonlinear least squares problems.

If you are interested in working on this, I think you should be able to use the geometries and tools available in Manopt to get a good head start on such a project. If you want help along the way, we will be happy to continue the conversation here. Actually, I have quite a few nonlinear least squares problems I'm currently working on, so I'd welcome some numerical insight here.

It looks like such a tool might take as input a function handle for the residue computation, as well as two function handles for operators representing the Jacobian of the residue function, both as is and its adjoint (Jv and J^Tw).

Best wishes,

Nicolas

[stum...@gmail.com]

Rather than implementing Levenberg-Marquardt, which is an older algorithm, it seems like you could modify your existing trust region algorithm to solve nonlinear least squares problems, i.e. min_x .5*r'(x)*r(x), where r(x) is a column vector of m residual functions evaluated at x \in R^n. Instead of defining a cost function, the problem would define the m residual functions and the gradient of each residual function; optionally, the gradients could be approximated via finite-differences if the gradients are not provided. Instead of computing the Hessian vector product H(x)*u in your trust region algorithm, you would instead compute [J(x)]'*(J(x)*u), where J(x) is the Jacobian of r(x) so that the rows of J(x) are the gradients of the residual functions. This technique works because [J(x)]'*J(x) is a very good approximation of the Hessian matrix H(x) for the special case when the cost function is a nonlinear least squares problem: .5*r'(x)*r(x).

Here is some information about MATLAB's nonlinear least squares trust region algorithm: https://www.mathworks.com/help/optim/ug/least-squares-model-fitting-algorithms.html#brrzgus

[stum...@gmail.com]

MATLAB's nonlinear least squares trust region algorithm doesn't work if the number m of residual functions is less than n, the dimension of x. In this case, Levenberg-Marquardt should be used instead.

https://www.mathworks.com/help/optim/ug/lsqnonlin.html#buumikd-6

[stum...@gmail.com]

If m>=n, it may be possible to realize a nonlinear least squares trust region algorithm immediately using Manopt's current trust region solver by defining the cost function to be .5*[r(x)]'*r(x), the gradient to be [J'(x)]*r(x), and the (approximate NLS) Hessian to be [J(x)]'*(J(x)*u), where r is the column vector of m residual functions and J is the m x n Jacobian of r. Obviously, lots of computational efficiency could be realized by using the Manopt store structure, saving r(x) and J(x).

[BM]

Modifying the Hess call to the trust-region solver was exactly what I had in my mind. We would be glad to make it work for your case. Let us know if you have further queries.

Regards,

On Monday, April 30, 2018 at 5:00:39 AM UTC+5:30, stum...@gmail.com wrote:

[Nicolas Boumal]

Hello,

That sounds nice -- you could have a tool that takes as input a manifold M (obtained from a factory) and function handles for r(x), J(x)u and J(x)'v (the matrix-vector product); it could either return a problem structure which can then be passed to trustregions, or it could call trustregions itself and return the answer -- if the latter, then the function should also take an options structure as input, to be passed to trustregions, and an initial guess. That's why it may be easier to do the former.

In any case, the problem should get problem.M = M and there should be a good implementation of cost / grad / approxhess that uses the store. Perhaps it is the given function handles that should use the store, for maximum flexibility.

(About problem.approxhess: if no Hessian is given, trustregions will look for an approximate Hessian first and will find this one; if it is not given, then finite differences of the gradient are used by default. -- it is an approximation of the Riemannian Hessian, so there should potentially be a call to ehess2rhess, which requires access to the Euclidean gradient too. Much of this could indeed be cached.)

Best,

Nicolas