The Gauss-Newton algorithm optimizes functions

$$ E(x) = \sum f(x)^2 $$

by approximating f as (locally) linear, in which case the Hessian of $E$ is approximated as

$$ H = 2 \sum {J_f}^T J_f $$

Now if I introduce a robust cost function in place of the squared cost above, I can similarly approximate the Hessian of $E$ using a linear approximation to $f$. For example, using the Cauchy robustifier

$$ E(x) = \sum \log\bigl(1+\frac{f^2}{\sigma^2}\bigr) $$

The hessian for $E$ assuming $f$ is linear is:

$$ H = \sum \frac{2 J^T J}{\sigma^2 + f^2} - \frac{4 (J^T f) (J^T f)^T}{(\sigma^2 + f^2)^2} $$

So is it a good idea to use this Hessian to solve the normal equations during gradient descent? Can I still use the Levenberg-Marquardt damping trick? Are there better options than this?

Some extra details of my problem:

- $x$ is small (6 dimensions parametrizing Fundamental matrices)
- $f$ is the Sampson error (a geometric error measure for two-image correspondences)
- I expect to have 100-400 residual terms in the summation
- I have a very low computation budget (a few milliseconds on a mobile device)