Levenberg-Marquadt near the minima for non-zero-residual problems - MathOverflow

Levenberg-Marquadt near the minima for non-zero-residual problems

2012-09-03T21:22:03Z

I'm using the LM algorithm to do gradient descent in a model fitting context. I'm minimizing:

$$ c(x) = \sum ( f_i(x) - y_i )^2 $$

I'm noticing that after a few steps when I'm close to the minima, I often wind up with a zero somewhere on the diagonal of the Hessian -- which I approximate in the usual Gauss-Newton way:

$$ H \approx J^T J $$

where $J$ is the jacobian of the residuals, $f_i$. Now, mathematically speaking, the Hessian will only have a zero on the diagonal if I've reached the true minimum. But with limited precision computers, one or more diagonal elements can easily go to zero, even when I'm quite far from the true optimum, since the curvature of different dimensions can easily differ by many orders of magnitude. This means that I wind up halting the gradient descent even when I'm not very close to the optimum.

As a side note, I am of course using the LM damping trick in which the diagonal is multiplied by $(1+\lambda)$, but this makes no difference when a diagonal element is numerical zero.

What is the correct thing to do here? Some options:

Switch to first-order gradient descent near the optimum
If one element on the diagonal is zero while the gradient is far from zero, then add epsilon to this element.
If one element on the diagonal is zero while the gradient is far from zero, then eliminate the corresponding parameter from optimization (i.e. delete the i-th row and col from $H$), and solve the normal equations for the remaining parameters.

Answer by Thomas Klimpel for Levenberg-Marquadt near the minima for non-zero-residual problems

2012-09-03T22:32:52Z

It's not true that you've reached the true minimum if the Hessian has a zero on the diagonal. It just means that the problem doesn't depend on the corresponding variable.

One popular LM implementation uses a QR decomposition (probably of $J$) with pivoting. If a diagonal element of $R$ becomes very small, it will be set to zero. This corresponds more or less to your last suggestion, with the difference that it is more general, actually works in practice, and is quite robust. If you are more interested in theoretical considerations, you could instead solve the linear least squares problem by a singular value decomposition, and set the singular values smaller than certain threshold to zero. This should clarify that this proceeding is the exact opposite of adding an epsilon to zero diagonal elements.

Answer by Brian Borchers for Levenberg-Marquadt near the minima for non-zero-residual problems

2012-09-04T00:33:19Z

The technique of multiplying the diagonal by $(1+\lambda)$ fails when you've got a 0 on the diagonal of $J^{T}J$. Having a zero on the diagonal of $J^{T}J$ can happen when you're far away from the optimal solution (in which case it tells you nothing about what might be true at optimality), or it can happen at an optimal solution (where the singularity of $J^{T}J$ can be a suggestion that the optimal solution is not unique.)

As Thomas stated in his answer, a common approach is to use the QR factorization to solve the least squares problem at each iteration and zero out small diagonal elements of $R$.

If you're stuck using the Cholesky factorization (e.g. if your problem is large and sparse and the Cholesky factorization can reasonably be computed but the QR factorization can't be computed), then you can avoid this problem with zeros on the diagonal by adding a small multiple of the identity to $J^{T}J$ rather than by multiplying the diagonal elements by $(1+\lambda)$.

You should also look carefully at the scaling of your problem- bad scaling frequently leads to poor convergence in practice.