# ForwardDiff gradient iterated quadqk

**URL:** <https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735>\
**Category:** Specific Domains\
**Tags:** question\
**Created:** [June 17, 2024, 6:32am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735 "2024-06-17T06:32:39Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 6:32am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/1 "2024-06-17T06:32:39Z")

</div>

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

```julia
# 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)

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 17, 2024, 6:41am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/2 "2024-06-17T06:41:01Z")

</div>

> [@Domenico\_Lahaye](#):
>
> Not sure about limitation (if any) of usage of ForwardDiff and quadgk. See MWE below.

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

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 8:02am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/3 "2024-06-17T08:02:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 17, 2024, 8:08am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/4 "2024-06-17T08:08:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 8:54am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/5 "2024-06-17T08:54:08Z")

</div>

Modified MWE to force single input argument. Code now is

```julia
# 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

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 17, 2024, 9:02am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/6 "2024-06-17T09:02:28Z")

</div>

> [@Domenico\_Lahaye](#):
>
> ` bfluxvec = zeros(N)`

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](https://juliadiff.org/ForwardDiff.jl/stable/user/limitations/)). 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:

```julia
bfluxvec = zeros(eltype(m), N)

```

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 11:27am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/7 "2024-06-17T11:27:05Z")

</div>

Many thx!

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 11:42am UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/8 "2024-06-17T11:42:30Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RayleighLord](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rayleighlord/32/33592_2.png) [@RayleighLord](https://discourse.julialang.org/u/RayleighLord)\
**Post date:** [June 17, 2024, 12:41pm UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/9 "2024-06-17T12:41:21Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [June 17, 2024, 2:37pm UTC](https://discourse.julialang.org/t/forwarddiff-gradient-iterated-quadqk/115735/10 "2024-06-17T14:37:47Z")

</div>

Thx! Much appreciated.