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.



Levenberg-Marquadt near the minima for non-zero-residual problems
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.


