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.