ForwardDiff gradient iterated quadqk

Specific Domains
1

Not sure about limitation (if any) of usage of ForwardDiff and quadgk. See MWE below.

# define the mesh with N elements e_k = [x_k, x_{k+1}] for 1 <= k <= N 
# the mesh contains N+1 nodes of which N-1 internal and 2 (left and right) boundary nodes  
N = 5; h = 1/N;
# define mesh points 
xmesh = Vector(0:h:1)
# define evaluation points 
xeval = Vector(0:h:1) 
# define midpoints 
xmid = Vector(h/2:h:1-h/2)

# define the magnetization coefficients  
m = ones(N)  

# define kernel  
function kernel(x,xp) 
    return x+xp 
end 

# compute integral by quadrature over second input (xp) over all elements e_k 
# on e_k the integral is x-dependent and stored in a array of functions 
function foo(x,m,xmesh)
   N = length(m) 
   fooarray = Array{Function}(undef, N)
   for k=1:N 
      fooarray[k] = (x,m,xmesh) -> quadgk(xp -> m[k]*kernel(x,xp), xmesh[k], xmesh[k+1])[1]
   end
   return sum([fooarray[k](x,m,xmesh) for k=1:N])
end 

# average of gradient of foo over e_k  
function avgfoo(m,xmesh)
    N = length(m) 
    bfluxvec = zeros(N)
    for k=1:N
        bfluxvec[k] = quadgk(x -> foo(x,m,xmesh), xmesh[k], xmesh[k+1];atol=1e-2, rtol=1e-2)[1] 
    end 
    return bfluxvec  
end

# evaluate foo in the nodes 
foosamp = foo.(xeval,Ref(m),Ref(xmesh)) 

# evaluate avgddxfoo on the elements 
avgfoosamp = avgfoo(m,xmesh)

p1 = plot(xeval, foosamp)
p2 = plot(xmid, avgfoosamp)
plot(p1,p2,layout=(1,2))

# compute the Jacobian wrt parameters m 
Cp = ForwardDiff.jacobian(avgfoo, m)
2

What’s the problem here, does it error? If so, can you provide the stack trace?

1 Like
3

MethodError: no method matching foo(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5}})
Closest candidates are:
foo(::Any, ::Any) at In[20]:4
foo(::Any, ::Any, ::Any) at In[85]:21

Stacktrace:
[1] vector_mode_dual_eval!(f::typeof(foo), cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5}}}, x::Vector{Float64})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24
[2] vector_mode_jacobian(f::typeof(foo), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
[3] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
[4] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(foo), Float64}, Float64, 5}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
[5] top-level scope
@ In[85]:51

4

Here the issue is that the Jacobian can only be computed for a function with one argument. You’re trying to do it with avgfoo(m, xmesh)m, which has two.

1 Like
5

Modified MWE to force single input argument. Code now is

# define the mesh with N elements e_k = [x_k, x_{k+1}] for 1 <= k <= N 
# the mesh contains N+1 nodes of which N-1 internal and 2 (left and right) boundary nodes  
N = 5; h = 1/N;
# define mesh points 
xmesh = Vector(0:h:1)
# define evaluation points 
xeval = Vector(0:h:1) 
# define midpoints 
xmid = Vector(h/2:h:1-h/2)

# define the magnetization coefficients  
m = ones(N)  

# define kernel  
function kernel(x,xp) 
    return x+xp 
end 

# compute integral by quadrature over second input (xp) over all elements e_k 
# on e_k the integral is x-dependent and stored in a array of functions 
function foo(x,m,xmesh)
   N = length(m) 
   fooarray = Array{Function}(undef, N)
   for k=1:N 
      fooarray[k] = (x,m,xmesh) -> quadgk(xp -> m[k]*kernel(x,xp), xmesh[k], xmesh[k+1])[1]
   end
   return sum([fooarray[k](x,m,xmesh) for k=1:N])
end 

# average of gradient of foo over e_k  
function avgfoo(m)
    N = length(m) 
    bfluxvec = zeros(N)
    for k=1:N
        bfluxvec[k] = quadgk(x -> foo(x,m,xmesh), xmesh[k], xmesh[k+1];atol=1e-2, rtol=1e-2)[1] 
    end 
    return bfluxvec  
end

# evaluate foo in the nodes 
foosamp = foo.(xeval,Ref(m),Ref(xmesh)) 

# evaluate avgddxfoo on the elements 
avgfoosamp = avgfoo(m) 

p1 = plot(xeval, foosamp)
p2 = plot(xmid, avgfoosamp)
plot(p1,p2,layout=(1,2))

# compute the Jacobian wrt parameters m 
Cp = ForwardDiff.jacobian(avgfoo, m)

Resulting stack trace is

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
…

Stacktrace:
[1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5})
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5}, i1::Int64)
@ Base ./array.jl:966
[3] avgfoo(m::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5}})
@ Main ./In[90]:35
[4] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24 [inlined]
[5] vector_mode_jacobian(f::typeof(avgfoo), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:125
[6] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:0
[7] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(avgfoo), Float64}, Float64, 5}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/jacobian.jl:19
[8] top-level scope
@ In[90]:51

6

Now the problem is here. ForwardDiff propagates a special type of number called Dual through your code to compute derivatives (see Limitations of ForwardDiff · ForwardDiff). It requires your code to be open to such generality, but here you create a vector bfluxvec of Float64. As a result, ForwardDiff cannot fill it with Dual numbers.
The fix here is to make the eltype of bfluxvec depend on the input:

bfluxvec = zeros(eltype(m), N)
1 Like
7

Many thx!

8

What are options to mitigate the single argument restriction that ForwardDiff imposes?

9

If you have a function f(x, p) that you want to differentiate with respect to x for some fixed parameters p, you can always do

der = ForwardDiff.jacobian(x -> f(x, p), x0)

which computes the jacobian at x0. Note that this is only creating an auxiliary function of one argument that has the function of several arguments encapsulated inside, but is usually convenient to use this shorthand notation.

1 Like
10

Thx! Much appreciated.