1
$\begingroup$

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.
Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ "It just means that the problem doesn't depend on the corresponding variable" isn't a correct statement. Consider the example of minimizing $(x_{1}^{2}-0)^{2}+(x_{2}^{2}-0)^{2}$. Clearly, $x=0$ is the unique optimal solution, but $J^{T}J$ is 0 at $x=0$. $\endgroup$
    – Brian Borchers
    Sep 4, 2012 at 0:40
  • $\begingroup$ @BrianBorchers Considering that AlexFlint writes that the zero is only occurring after a few steps, your explanation might be closer to what is actually happening in his case than my suggestion. However, if this should be the case, then there is something strange going on. I believe this can happen if the model fitting problem degenerates into an optimization problem, i.e. if the model can't really fit the measurement data, and has to try to maximize (or minimize) the model component which can't be fitted. But note that LM is not very efficient for this type of problem. $\endgroup$
    – Thomas Klimpel
    Sep 4, 2012 at 8:37
  • $\begingroup$ @ThomasKlimpel Thanks, this is helpful. Is there a text you can recommend on this topic? I've consulted a couple of convex optimization books but they (understandably) do not go into this level of detail about this specific setting. Also, can you elaborate on "LM is not very efficient for this type of problem" - my problem is certainly overdetermined and I do not expect a perfect fit. I also expect there to be (infrequent) outliers, which I expect to have large residuals at the optimum (I am using a robust cost function for this reason). $\endgroup$
    – Alex Flint
    Sep 5, 2012 at 14:07
  • $\begingroup$ @Alex; The $J^{T}J$ approximation to the hessian of $c(x)$ works well when $f_{i}(x)-y_{i}$ is small, but if the residuals are large then the approximation can degrade. This is because the second order terms that go into the Hessian are dropped when we approximate it with $J^{T}J$, and those terms have $f_{i}(x)-y_{i}$ factors. $\endgroup$
    – Brian Borchers
    Sep 5, 2012 at 18:01
2
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ By the way, it's easy to compute the gradient of $c(x)$ given $J$ and the values of $f_{i}(x)-y_{i}$. You should always check that the gradient is reasonably small as part of your termination criteria. $\endgroup$
    – Brian Borchers
    Sep 4, 2012 at 4:19
  • $\begingroup$ @BrianBorchers. Thanks, yes that is definitely something I should be doing. I'm aware of scaling issues but unsure precisely how to scale my problem given that I'm using a robust error function. What is the relevant quantity to normalize in general? $\endgroup$
    – Alex Flint
    Sep 5, 2012 at 16:31
  • $\begingroup$ This is confusing- are you minimizing a sum of squares or are you using some other objective (such as the Huber measure)? If you're not minimizing a sum of squares than LM isn't appropriate in the first place. Assuming that you are doing nonlinear least squares, and assuming that you have a reasonable assessment of the uncertainty in the measurements (e.g. $y_{i}$ is known with uncertainty $\sigma_{i}$, then you should normalize by $c(x)=\sum ((f_{i}(x)-y_{i})/\sigma_{i})^{2}$. You should also scale the parameters $x$ so that they're all of similar magnitude. $\endgroup$
    – Brian Borchers
    Sep 5, 2012 at 17:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.