ForwardDiff errors when evaluating type: Array{Float64}

My questions are asked at the bottom of this post, but their background is developed through this post - which is perhaps a bit pedantic.

I’m writing a p-method FEM code. Instead of explicitly coding high-order shape functions and Jacobian matrices, I’m writing basis function routines and utilizing auto-differentiation via ForwardDiff. Here’s my Lagrange basis function routine.

``````function lagrangeBasis(k,variate)
node = range(-1,1,length=k+1)
L = Array{Float64,1}(undef,k+1)    ## Question regards this line
for j = 1:k+1
L[j] = 1.0
for m = 1:k+1
if j != m
L[j] *= ((variate-node[m])/(node[j]-node[m]))
end
end
end
return L
end
``````

Calling this function returns an array of `Float64` values for each basis function, in this case the 9 basis functions associated with the degree 8 Lagrange basis:

``````julia> lagrangeBasis(8,0.)
9-element Array{Float64,1}:
0.0
-0.0
0.0
-0.0
1.0
0.0
-0.0
0.0
-0.0
``````

And benchmarking this code is giving me acceptable performance:

``````julia> BenchmarkTools.@btime(lagrangeBasis(8,\$rand()));
417.225 ns (1 allocation: 160 bytes)
``````

The problem I have is when I try to evaluate the derivatives (or Jacobians in multi-variate cases) using the ForwardDiff package:

``````julia> ForwardDiff.derivative(x->lagrangeBasis(8,x),rand())
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##20#21")),Float64},Float64,1})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:741
Float64(::Int8) at float.jl:60
...
Stacktrace:
[1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##20#21")),Float64},Float64,1}) at .\number.jl:7
[2] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##20#21")),Float64},Float64,1}, ::Int64) at .\array.jl:767
[3] lagrangeBasis(::Int64, ::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##20#21")),Float64},Float64,1}) at .\REPL[23]:8
[4] #20 at .\REPL[29]:1 [inlined]
[5] derivative(::getfield(Main, Symbol("##20#21")), ::Float64) at C:\Users\gregj\.julia\packages\ForwardDiff\N0wMF\src\derivative.jl:14
[6] top-level scope at none:0
``````

I’ve found a workaround, which is to replace the array initialization in the basis function routine:

``````function lagrangeBasis(k,variate)
node = range(-1,1,length=k+1)
L = Array{Any,1}(undef,k+1)  ## MODIFIED
for j = 1:k+1
L[j] = 1.0
for m = 1:k+1
if j != m
L[j] *= ((variate-node[m])/(node[j]-node[m]))
end
end
end
return L
end
``````

Now evaluating this routine returns an array of `Any` values

``````julia> lagrangeBasis(8,0.)
9-element Array{Any,1}:
0.0
-0.0
0.0
-0.0
1.0
0.0
-0.0
0.0
-0.0
``````

And it takes considerably longer (~6x) to evaluate:

``````julia> BenchmarkTools.@btime(lagrangeBasis(8,\$rand()));
2.517 μs (145 allocations: 2.41 KiB)
``````

But at least now I can compute the derivatives / Jacobians:

``````julia> ForwardDiff.derivative(x->lagrangeBasis(8,x), 0.)
9-element Array{Float64,1}:
0.014285714285714285
-0.15238095238095237
0.7999999999999999
-3.2
-1.1102230246251565e-15
3.2
-0.7999999999999998
0.15238095238095237
-0.014285714285714284
``````

From these examples, I have two questions

1. Is there a simple & performant way to use the ForwardDiff package when the evaluated function returns an array of `Float64`, such as in the first `lagrangeBasis` function I provided?
2. Maybe a quick solution to the first question is to use multiple dispatch on the `lagrangeBasis` function: write one version that gets called when simply evaluating the function, and another version that gets called if ForwardDiff is evaluating the function:
``````function lagrangeBasis(k,variate::Float64)
node = range(-1,1,length=k+1)
L = Array{Float64,1}(undef,k+1)
...
end

function lagrangeBasis(k,variate::ForwardDiffType)  # <- ??????
node = range(-1,1,length=k+1)
L = Array{Any,1}(undef,k+1)
...
end
``````

What type would I need to pass into `lagrangeBasis` to make this second case work?

1 Like

By doing this you not only don’t allow `L` to accept dual numbers, but you also don’t adapt to other different number types — e.g. what if `variate` is a different precision and you want to keep the values to that precision, or what if `variate` is a dimensionful quantity?

In general, the strength of Julia is that you can write code that is simultaneously fast and type-generic, and it’s well worth making an effort to achieve this when you are writing reusable library code.

One option would be to do something like:

``````node = range(-oneunit(variate), oneunit(variate), length=k+1)
L = Vector{typeof(zero(variate)/oneunit(variate))}(undef, k+1)
``````

so that you let Julia compute the type of the entries when allocating `L`, and so that it uses the type of `variate` to compute the type of `node`.

Note, by the way, that you might look into the barycentric formula: https://people.maths.ox.ac.uk/trefethen/barycentric.pdf (and note especially that equally spaced nodes are usually a poor choice).

3 Likes

I wrote a mini-package

for type calculations like these.