20

Let's say I'm computing some model parameters my minimizing the sum squared residuals, and I'm assuming my errors are Gaussian. My model produces analytical derivatives, so the optimizer does not need to use finite differences. Once the fit is complete, I want to calculate standard errors of the fitted parameters.

Generally, in this situation, the Hessian of the error function is taken to be related to the covariance matrix by: $$ \sigma^2 H^{-1} = C $$ where $\sigma^2$ is the variance of the residuals.

When no analytical derivatives of the error are available, it is typically impractical to compute the Hessian, so $J^TJ$ is taken as a good approximation.

However, in my case, I've got an analytical J, so it's relatively cheap for me to compute H by finite differencing J.

So, my question is this: Would it be more accurate to approximate H using my exact J and applying the above approximation, or to approximate H by finite differencing J?

Colin K
  • 477
  • 3
  • 9

1 Answers1

12

GOOD question. First, recall where this approximation $H \approx J^T J$ comes from. Let $(x_i, y_i)$ be your data points, $f(\cdot)$ be your model and $\beta$ be the parameters of your model. Then the objective function of the non-linear least squares problem is $\frac{1}{2} r^T r$ where $r$ is the vector of the residuals, $r_i = y_i - f(x_i, \beta)$. The exact Hessian of the objective function is $H = J^T J + \sum r_i \nabla^2 r_i$. So the error in this approximation is $H - J^T J = \sum r_i \nabla^2 r_i$. It's a good approximation when the residuals, themselves, are small; or when the 2nd derivative of the residuals is small. Linear least squares can be considered a special case where the 2nd derivative of the residuals is zero.

As for finite difference approximation, it is relatively cheap. To compute a central difference, you'll need to evaluate the Jacobian an additional $2n$ times (a forward difference will cost you $n$ additional evaluations, so I wouldn't bother). The error of the central difference approximation is proportional to $\nabla^4 r$ and $h^2$, where $h$ is the step size. The optimal step size is $h \sim \epsilon^\frac{1}{3}$, where $\epsilon$ is machine precision. So unless the derivatives of the residuals are blowing up, it's pretty clear that the finite difference approximation should be a LOT better. I should point out that, while the computation is minimal, the bookkeeping is nontrivial. Each finite difference on the Jacobian will give you one row of the Hessian for each residual. You'll then have to reassemble the Hessian using the formula above.

There is, however, a 3rd option. If your solver uses a Quasi-Newton method (DFP, BFGS, Bryoden, etc.), it is already approximating the Hessian at each iteration. The approximation can be quite good, as it uses the objective function and gradient values from every iteration. Most solvers will give you access to the final Hessian estimate (or its inverse). If that's an option for you, I would use that as the Hessian estimate. It's already computed and it's probably going to be a pretty good estimate.

Bill Woessner
  • 1,168
  • 11
  • 15
  • Excellent response, thank you. Justifying it with a comparison of the estimation error in each case is very enlightening. Can I ask how you know that $\epsilon^{1/3}$ is the optimal step for finite differences? I've never seen that before. – Colin K Oct 21 '15 at 19:51
  • 5
    That's an old trick to balance truncation error vs. round-off error. Obviously, to minimize truncation error, you want to make $h$ as small as possible. But once $h$ gets too small, you start to incur significant round-off error. The derivation is relatively straightforward. Assuming a central difference, the truncation error is proportional to $h^2 f'''(x)$. The round-off error is always proportional to $\frac{\epsilon f(x)}{h}$. Add the two and minimize over $h$. You get $h \sim \epsilon^\frac{1}{3}$. – Bill Woessner Oct 22 '15 at 14:02
  • 3
    This only holds for central differences. For forward differences, the optimal step size is $h \sim \epsilon^\frac{1}{2}$. There are other tricks, as well. For example, make sure you actually know what $h$ is. I know this sounds silly, but weird things can happen in floating point arithmetic. Here's a simple way to make sure you have the correct value of $h$: `h_actual = (x + h_desired) - x`. Mathematically, of course, $h_{actual} = h_{desired}$. But if you use values that cannot be exactly represented in floating point (like $h = 0.0001$), you'll see that is not the case. – Bill Woessner Oct 22 '15 at 14:33
  • Perhaps this content could be added to your answer, rather than the comments. That way, future users don't have to wade through an extended comment section to find material that directly bears on claims made in the answer. – Sycorax Oct 22 '15 at 15:35
  • 3
    Oh my goodness. A Quasi-Newton approximation of the Hessian can be a terrible estimate of the Hessian, and therefore result in a very poor estimate of the covariance matrix. It may serve well to facilitate progression of the algorithm to the optimum, but can be quite poor as an estimate of the Hessian. – Mark L. Stone Oct 29 '15 at 22:58