Accuracy is a bit brittle

[giordano]

The README of this package convinced me to use FiniteDifferences instead of Calculus in Measurements for cases where I need a numerical derivative (for functions without known derivative or automatic differentiation isn't possible, e.g. calling into C libraries). However I found that accuracy of FiniteDifferences is quite brittle (sometimes strongly depends on the number of points with big even/odd jumps, see example below) and I need to use a largish number of grid points to achieve the same accuracy I have with Calculus, which in comparison is also much faster (and makes me question the change of the backend). For example, the derivative of cosc in 0 is -(pi^2)/3, with Julia v1.5.3 I get:

julia> using Calculus, FiniteDifferences, BenchmarkTools

julia> -(pi ^ 2) / 3
-3.289868133696453

julia> (p -> FiniteDifferences.central_fdm(p, 1)(cosc, 0)).(2:10)
9-element Array{Float64,1}:
 -3.2575126788312536 # 2
  0.0                # 3
 -3.2894144933275746 # 4
 -3.9730287078618036 # 5 
 -3.289860720421136  # 6
 -3.2898376423854927 # 7
 -3.2898681297554884 # 8
 -3.289868469475128  # 9
 -3.2898681348743732 # 10

julia> Calculus.derivative(cosc, 0) 
-3.2898678494407765

To achieve an accuracy of the same order as Calculus I need 8 grid points, and this is the performance comparison:

julia> @btime FiniteDifferences.central_fdm(8, 1)($cosc, 0)
  2.340 μs (25 allocations: 992 bytes)
-3.2898681297554884

julia> @btime Calculus.derivative($cosc, 0)
  60.591 ns (0 allocations: 0 bytes)
-3.2898678494407765

Is this behaviour really expected?

[nickrobinson251]

Doesn't sound expected to me. Would you mind trying with FiniteDifferences v0.10, to see if this changed between versions?

cc @wesselb

[wesselb]

Ouch, this looks terrible. A quick investigation shows that the step size estimation behaves badly:

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(2, 1), cosc, 0.0)
(0.1, 1.005e-39)

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(3, 1), cosc, 0.0)
(2.502562156920058e-14, 5.993857119001894e-27)

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(4, 1), cosc, 0.0)
(0.1, 1.5000166666666666e-39)

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(5, 1), cosc, 0.0)
(5.502059929042815e-9, 3.4078145715984384e-32)

I'm suspecting that the issue is that both the argument is 0.0 and cosc(0.0) = 0.0.

Fixing the former appears to resolve the issue:

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(4, 1), cosc, 1.0)
(7.595498177230801e-5, 5.846742366173149e-12)

julia> FiniteDifferences.central_fdm(4, 1)(cosc, 1) - 2
4.800604358479177e-13

julia> Calculus.derivative(cosc, 1) - 2
-1.0873346667494843e-10

Moreover, fixing the latter also appears to resolve the issue:

julia> FiniteDifferences.estimate_step(FiniteDifferences.central_fdm(4, 1), x -> cosc(x) + 1, 0.0)
(0.0009866749526150414, 4.500866355967255e-13)

julia> FiniteDifferences.central_fdm(4, 1)(x -> cosc(x) + 1, 0) + (pi ^ 2) / 3
4.785150053976395e-11

julia> Calculus.derivative(cosc, 0) + (pi ^ 2) / 3
2.8425567633050264e-7

I'll further look into this.

Performance of the package is definitely something that at some point needs to be addressed.

[wesselb]

I should also say that the package since recently supports Richardson extrapolation with Richardson.jl, which is an alternative that does not break in this case:

julia> @btime FiniteDifferences.extrapolate_fdm(FiniteDifferences.central_fdm(2, 1, adapt=0), f, 0.0)[1] + (pi ^ 2) / 3
  1.898 μs (29 allocations: 1.44 KiB)
-3.455082886461014e-9

[giordano]

I confirm #125 addresses my complaints, thanks @wesselb!!