johnmyleswhite/Calculus.jl — pull request #4 by timholy (2013-01-24 17:32:01+00:00)
Title: Improved finite differencing

I made some improvements in the finite differencing. Here are some brief explanations:

- Note that the new formula, `2 * sqrt(1e-12) * (1 + norm(x))`, for `abs(x) < 1` is approximately equal to 2e-6. The formula I used previously, `eps(max(1.0, abs(x)))^(1/3)`, for `abs(x) < 1` is equal to 6e-6, not very different. So contrary to your in-code comment, the small `x` behavior actually can't be very different for these two.
- `1e-12` implicitly assumes `Float64`, whereas using the type-dispatched version of `eps` does not. I definitely run optimization with `Float32`s sometimes (when using huge datasets where all the computations are `Float32`).
- The proper power scaling depends on the order of the derivative, and whether centered or forward differencing is used. You'll see the power I used in the revised version is 1/2, 1/3, or 1/4 depending on the circumstance. This is particularly complex for the Hessian.

I also removed the dtype parameter from `finite_difference_hessian(f, x)` because you basically have to use centered differencing for direct Hessian evaluation.

For your cases in `test/finite_difference.jl`, I think you'll mostly see improved accuracy. For central differencing the change is quite modest (usually less than one order of magnitude), but for forward differencing and the Hessian it can be dramatic. For example, with the old version
```
julia> norm(Calculus.finite_difference_hessian(fx, gx, [0.0, 0.0], :central) - [-sin(0.0) 0.0; 0.0 -cos(0.0)])
6.938866148331613e-6
```
and with the new
```
julia> norm(Calculus.finite_difference_hessian(fx, [0.0, 0.0]) - [-sin(0.0) 0.0; 0.0 -cos(0.0)])
7.450580596923828e-9
```
Presumably it should also be more robust in real-world problems, although that remains to be seen under usage.