gradient only accurate to 10 digits ?

[GravityAssisted]

Here is a function I like to calculate a gradient of:

function getW33(kv :: Vector)
    k = kv[1]
    sqrt2l= sqrt(2.0)
          t2 = k*k
          t3 = acos(t2 - 1.0)
          t4 = 2.0 - t2
          t6 = 2*π - t3
          t5 = 1.0/t4
          W =   (t6 * sqrt(t5) - k)*t5
    return W
end

Its actual derivative at kv = [-0.5] is -2.8185256628482382
but using the function gradient I get an answer which is only 10 digits accurate:

g = gradient(getW33)
g([-0.5])[1]
> -2.8185256629946482

Using complex step derivative I get full 16 digits of accuracy.

kv = [complex(-0.5,1E-30)]
imag(getW33(kv)) / 1E-30
> -2.8185256628482387

I was under the impression that the ForwardDiff was accurate to machine precision, is that wrong ?

I am on Julia 0.4, MacOSX, LLVM 3.3

thanks,
Nitin

[papamarkou]

Hi @GravityAssisted, forward mode AD avoids the truncation errors of numerical differentiation, and is in general accurate. It may get affected by the actual implementation of the evaluated function, as this may cause propagation of machine precision errors, so that's where I would start from. That being said, forward mode has been proven to be backward stable in the sense of Wilkinson, which means that even small perturbations of the original function due to machine precision eps should still yield accurate derivatives. A reference for this theoretical point is perhaps Griewank's recent paper (2014), see here.

[GravityAssisted]

@scidom I just did a @code_warntype on the function. Looks like its calling a finite differencing function ??

 @code_warntype g([-0.5])
Variables:
  x::Array{Float64,1}
  ##g#9007::Array{Float64,1}

Body:
  begin  # /Users/arora/.julia/v0.4/Calculus/src/derivative.jl, line 5:
      GenSym(0) = (Base.arraylen)(x::Array{Float64,1})::Int64
      ##g#9007 = (top(ccall))(:jl_alloc_array_1d,(top(apply_type))(Base.Array,Float64,1)::Type{Array{Float64,1}},(top(svec))(Base.Any,Base.Int)::SimpleVector,Array{Float64,1},0,GenSym(0),0)::Array{Float64,1}
      (Calculus.finite_difference!)(f::F,x::Array{Float64,1},##g#9007::Array{Float64,1},dtype::Symbol)::Void
      return ##g#9007::Array{Float64,1}
  end::Array{Float64,1}

That would explain the 10 digits of precision. Roundoff errors cant be that bad on this function.

Seems fishy to me.

[GravityAssisted]

Just found the error... The code is not even calling the ForwardDiff library !! Its calling the inbuilt Julia gradient function. Very suspicious that It didn't warn me of that when I did "using ForwardDiff".

Now, if I do the following I match up-to 16 digits.

g = ForwardDiff.gradient(getW33)
g([-0.5])[1]
> -2.8185256628482382

The module should be warning me of this or this is something that should be in documentation as I am sure lot of people might be using the gradient function directly and accumulating errors needlessly.

[KristofferC]

erm: http://www.juliadiff.org/ForwardDiff.jl/perf_diff.html#gradients

When calling ForwardDiff.jl’s gradient method, you must use the fully qualified method name 
ForwardDiff.gradient(...) in order to avoid conflict with Base.gradient.

[papamarkou]

Great, good to see you spotted the error @GravityAssisted; since the warning is already documented as @KristofferC mentioned I suppose we can close this issue.

[GravityAssisted]

@KristofferC , @scidom thanks, my bad. I should have read the documentation more carefully.

[papamarkou]

Glad to see that you got it sorted @GravityAssisted; if not mistaken, @jrevels has updated the examples to use import instead of using in an attempt to eliminate namespace conflict errors.

[jrevels]

@jrevels has updated the examples to use import instead of using in an attempt to eliminate namespace conflict errors

Yup. It's a very easy mistake to make (I still make it myself sometimes when doing quick hacks in the REPL). I've taken to always using import instead of using in library code in order to mitigate accidental scope issues.