ForwardDiff errors when evaluating type: Array{Float64}

Specific Domains Numerics
1

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
2

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
3

I wrote a mini-package

for type calculations like these.