# Array issue with optimzing a function with quadgk

**URL:** <https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837>\
**Category:** Optimization (Mathematical)\
**Created:** [June 22, 2020, 4:45am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837 "2020-06-22T04:45:04Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [June 22, 2020, 4:45am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/1 "2020-06-22T04:45:04Z")

</div>

I’m trying to minimize a function that involves an integral, the integrand of which is one of the optimizing parameters. I can compute the function and get the gradient with Zygote, but I’m getting an odd conversion error when I try to run an optimization routine. Might anyone know what is going on?

```julia
using Distributions
using Random
using ForwardDiff
using Zygote
using QuadGK
using Optim

theta_lower = 1.0
theta_upper = 5.0
uniform_dist = Uniform(theta_lower, theta_upper)

function linearCost(θ)
    
    out = 2.5 - 0.3*θ
    
    return(out)
    
end

function getAlpha(θ₁, θ₂, B, costFn, Fdist)
    
    intcdF1 = quadgk(θ -> (linearCost(θ)*pdf(Fdist,θ)), θ₁, θ₂, rtol=1e-8)[1]
    intcdF2 = quadgk(θ -> (linearCost(θ)*pdf(Fdist,θ)), θ₂, theta_upper, rtol=1e-8)[1]
    
    alpha_num = B + θ₂*(1.0-cdf(Fdist,θ₂)) - intcdF2
    alpha_denom = intcdF1 + θ₂*(1.0-cdf(Fdist,θ₂)) - θ₁*(1-cdf(Fdist,θ₁))
    
    α = alpha_num / alpha_denom
    
    return(α)
    
end

function socialPlannerObjective(x, params)
    
    θ₁, Δθ = x
    B, costFn, Fdist = params
    
    θ₂ = θ₁ + Δθ
    
    intdiffcdF1= quadgk(θ -> ((θ-linearCost(θ))*pdf(Fdist,θ)), θ₁, θ₂, rtol=1e-8)[1]
    intdiffcdF2 = quadgk(θ -> ((θ-linearCost(θ))*pdf(Fdist,θ)), θ₂, theta_upper, rtol=1e-8)[1]
    
    alpha = getAlpha(θ₁, θ₂, B, costFn, Fdist)
    
    @show out = alpha*intdiffcdF1 + intdiffcdF2 
    
    return(out)
    
end
    
gradsocialPlannerObjective(y) = Zygote.gradient(x -> socialPlannerObjective(x,[10, linearCost, uniform_dist]), y)[1]

## test
@show socialPlannerObjective([1.3, 0.3], [10, linearCost, uniform_dist])
@show gradsocialPlannerObjective([1.3, 0.6])

optimize(x -> socialPlannerObjective(x,[10, linearCost, uniform_dist]), x->gradsocialPlannerObjective(x), [1.3, 0.6], LBFGS(); inplace = false)

```

Output being received:

```julia
out = alpha * intdiffcdF1 + intdiffcdF2 = 0.03361344537815181
socialPlannerObjective([1.3, 0.3], [10, linearCost, uniform_dist]) = 0.03361344537815181
out = alpha * intdiffcdF1 + intdiffcdF2 = 0.40183246073298395
gradsocialPlannerObjective([1.3, 0.6]) = [2.7550300156245733, 1.3382479921054806]
out = alpha * intdiffcdF1 + intdiffcdF2 = 0.40183246073298395
out = alpha * intdiffcdF1 + intdiffcdF2 = 0.40183246073298395

MethodError: Cannot `convert` an object of type Array{QuadGK.Segment{Float64,Float64,Float64},1} to an object of type QuadGK.Segment{Float64,Float64,Float64}
Closest candidates are:
  convert(::Type{T}, !Matched::QuadGK.Segment) where T<:QuadGK.Segment at /Users/svass/.julia/packages/QuadGK/jmDk8/src/evalrule.jl:10
  convert(::Type{T}, !Matched::T) where T at essentials.jl:168
  QuadGK.Segment{Float64,Float64,Float64}(::Any, !Matched::Any, !Matched::Any, !Matched::Any) where {TX, TI, TE} at /Users/svass/.julia/packages/QuadGK/jmDk8/src/evalrule.jl:3

Stacktrace:
 [1] setindex!(::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Int64) at ./array.jl:782
 [2] adjoint at /Users/svass/.julia/packages/Zygote/YeCEW/src/lib/array.jl:60 [inlined]
 [3] _pullback at /Users/svass/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:47 [inlined]
 [4] percolate_down! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:28 [inlined]
 [5] _pullback(::Zygote.Context, ::typeof(DataStructures.percolate_down!), ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Int64, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}, ::Int64) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [6] percolate_down! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:18 [inlined]
 [7] _pullback(::Zygote.Context, ::typeof(DataStructures.percolate_down!), ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Int64, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [8] heappop! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:59 [inlined]
 [9] adapt at /Users/svass/.julia/packages/QuadGK/jmDk8/src/adapt.jl:36 [inlined]
 [10] _pullback(::Zygote.Context, ::typeof(QuadGK.adapt), ::var"#64#66"{Uniform{Float64}}, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Float64, ::Float64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [11] do_quadgk at /Users/svass/.julia/packages/QuadGK/jmDk8/src/adapt.jl:28 [inlined]
 [12] _pullback(::Zygote.Context, ::typeof(QuadGK.do_quadgk), ::var"#64#66"{Uniform{Float64}}, ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [13] #28 at /Users/svass/.julia/packages/QuadGK/jmDk8/src/adapt.jl:159 [inlined]
 [14] handle_infinities at /Users/svass/.julia/packages/QuadGK/jmDk8/src/adapt.jl:93 [inlined]
 [15] _pullback(::Zygote.Context, ::typeof(QuadGK.handle_infinities), ::QuadGK.var"#28#29"{Nothing,Float64,Int64,Int64,typeof(LinearAlgebra.norm)}, ::var"#64#66"{Uniform{Float64}}, ::Tuple{Float64,Float64}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [16] #quadgk#27 at /Users/svass/.julia/packages/QuadGK/jmDk8/src/adapt.jl:157 [inlined]
 [17] adjoint at /Users/svass/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:168 [inlined]
 [18] _pullback at /Users/svass/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:47 [inlined]
 [19] #quadgk at ./none:0 [inlined]
 [20] _pullback(::Zygote.Context, ::QuadGK.var"#kw##quadgk", ::NamedTuple{(:rtol,),Tuple{Float64}}, ::typeof(quadgk), ::var"#64#66"{Uniform{Float64}}, ::Float64, ::Float64) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [21] socialPlannerObjective at ./In[8]:42 [inlined]
 [22] _pullback(::Zygote.Context, ::typeof(socialPlannerObjective), ::Array{Float64,1}, ::Array{Any,1}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [23] #67 at ./In[8]:52 [inlined]
 [24] _pullback(::Zygote.Context, ::var"#67#68", ::Array{Float64,1}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [25] _pullback(::Function, ::Array{Float64,1}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface.jl:39
 [26] pullback(::Function, ::Array{Float64,1}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface.jl:45
 [27] gradient(::Function, ::Array{Float64,1}) at /Users/svass/.julia/packages/Zygote/YeCEW/src/compiler/interface.jl:54
 [28] gradsocialPlannerObjective(::Array{Float64,1}) at ./In[8]:52
 [29] #70 at ./In[8]:58 [inlined]
 [30] (::NLSolversBase.var"#gg!#2"{var"#70#72"})(::Array{Float64,1}, ::Array{Float64,1}) at /Users/svass/.julia/packages/NLSolversBase/mGaJg/src/objective_types/inplace_factory.jl:21
 [31] (::NLSolversBase.var"#fg!#8"{var"#69#71",NLSolversBase.var"#gg!#2"{var"#70#72"}})(::Array{Float64,1}, ::Array{Float64,1}) at /Users/svass/.julia/packages/NLSolversBase/mGaJg/src/objective_types/abstract.jl:13
 [32] value_gradient!!(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /Users/svass/.julia/packages/NLSolversBase/mGaJg/src/interface.jl:82
 [33] value_gradient!(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at /Users/svass/.julia/packages/NLSolversBase/mGaJg/src/interface.jl:69
 [34] value_gradient!(::Optim.ManifoldObjective{OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}}, ::Array{Float64,1}) at /Users/svass/.julia/packages/Optim/L5T76/src/Manifolds.jl:50
 [35] (::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}},Array{Float64,1},Array{Float64,1},Array{Float64,1}})(::Float64) at /Users/svass/.julia/packages/LineSearches/WrsMD/src/LineSearches.jl:84
 [36] (::HagerZhang{Float64,Base.RefValue{Bool}})(::Function, ::LineSearches.var"#ϕdϕ#6"{Optim.ManifoldObjective{OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}},Array{Float64,1},Array{Float64,1},Array{Float64,1}}, ::Float64, ::Float64, ::Float64) at /Users/svass/.julia/packages/LineSearches/WrsMD/src/hagerzhang.jl:140
 [37] HagerZhang at /Users/svass/.julia/packages/LineSearches/WrsMD/src/hagerzhang.jl:101 [inlined]
 [38] perform_linesearch!(::Optim.LBFGSState{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1},Float64,Array{Float64,1}}, ::LBFGS{Nothing,InitialStatic{Float64},HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.ManifoldObjective{OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}}) at /Users/svass/.julia/packages/Optim/L5T76/src/utilities/perform_linesearch.jl:56
 [39] update_state!(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Optim.LBFGSState{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1},Float64,Array{Float64,1}}, ::LBFGS{Nothing,InitialStatic{Float64},HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}) at /Users/svass/.julia/packages/Optim/L5T76/src/multivariate/solvers/first_order/l_bfgs.jl:198
 [40] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::LBFGS{Nothing,InitialStatic{Float64},HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}, ::Optim.LBFGSState{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1},Float64,Array{Float64,1}}) at /Users/svass/.julia/packages/Optim/L5T76/src/multivariate/optimize/optimize.jl:57
 [41] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::LBFGS{Nothing,InitialStatic{Float64},HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at /Users/svass/.julia/packages/Optim/L5T76/src/multivariate/optimize/optimize.jl:33
 [42] #optimize#95 at /Users/svass/.julia/packages/Optim/L5T76/src/multivariate/optimize/interface.jl:129 [inlined]
 [43] (::Optim.var"#kw##optimize")(::NamedTuple{(:inplace,),Tuple{Bool}}, ::typeof(optimize), ::Function, ::Function, ::Array{Float64,1}, ::LBFGS{Nothing,InitialStatic{Float64},HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at ./none:0 (repeats 2 times)
 [44] top-level scope at In[8]:57

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 22, 2020, 5:35am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/2 "2020-06-22T05:35:51Z")

</div>

Instead of treating this as a technical problem, I would recommend that you consider the underlying conceptual issue: QuadGK uses _adaptive_ quadrature, which is tricky to differentiate (it may be differentiable locally, but not continuous).

I would recommend using simple Gaussian quadrature instead, eg

> **[GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for...](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)**
>
> Julia package for Gaussian quadrature. Contribute to JuliaApproximation/FastGaussQuadrature.jl development by creating an account on GitHub.

In particular, if your interval is finite, you should use [Gauss-Legendre](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature) with fixed nodes/weights, and transform according to `θ₁`, `θ₂` etc. Then integration is simply a dot product, which is differentiable.

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [August 17, 2020, 8:41pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/3 "2020-08-17T20:41:54Z")

</div>

> [@Tamas\_Papp](#):
>
> transform according to `θ₁` , `θ₂` etc.

Do you have any tips on how to do this?

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [August 17, 2020, 10:15pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/4 "2020-08-17T22:15:15Z")

</div>

Nevermind, I rolled my own. Ideally, FastGaussQuadrature would allow for arbitrary intervals [a,b] (and I may submit a PR along these lines), but in absence of that, [here](https://github.com/OU-PhD-Econometrics/fall-2020/blob/master/ProblemSets/PS4-mixture/lgwt.jl) is an implementation (taken from MATLAB’s [file exchange](https://www.mathworks.com/matlabcentral/fileexchange/4540-legendre-gauss-quadrature-weights-and-nodes)).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [August 18, 2020, 11:13am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/5 "2020-08-18T11:13:47Z")

</div>

> [@tyleransom](#):
>
> Ideally, FastGaussQuadrature would allow for arbitrary intervals [a,b] (and I may submit a PR along these lines)

I think the intention of that package is to support “canonical” intervals for each quadrature, which the user can transform easily.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 18, 2020, 12:25pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/6 "2020-08-18T12:25:34Z")

</div>

> [@Tamas\_Papp](#):
>
> Instead of treating this as a technical problem, I would recommend that you consider the underlying conceptual issue: QuadGK uses _adaptive_ quadrature, which is tricky to differentiate (it may be differentiable locally, but not continuous).

If you assume that the tolerance is set low enough that QuadGK is computing a good approximation of the exact integral, then you could express the derivative of the integral (via the [Leibniz rule](https://en.wikipedia.org/wiki/Leibniz_integral_rule)) as another integral and then compute that integral via QuadGK. (This is sometimes called a “differentiate then discretize” approach.)

A PR to provide this interface to QuadGK via [ChainRules.jl](https://github.com/JuliaDiff/ChainRules.jl), so that it integrates with Zygote etcetera, would be great.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 18, 2020, 12:35pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/7 "2020-08-18T12:35:14Z")

</div>

Quadrature.jl wraps the quadrature libraries and adds a derivative rule: [https://github.com/SciML/Quadrature.jl](https://github.com/SciML/Quadrature.jl) .

```julia
using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba
f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
dp1 = Zygote.gradient(testf,p)
dp2 = FiniteDiff.finite_difference_gradient(testf,p)
dp3 = ForwardDiff.gradient(testf,p)
dp1[1] ≈ dp2 ≈ dp3

```

We currently don’t handle derivatives w.r.t. the bounds of the integral though, but that’s something we plan to add soon (also, there’s a slightly better way to compute what we’re doing since we aren’t caching the primal evaluation, but 🤷‍♂️)

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [August 19, 2020, 2:53pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/8 "2020-08-19T14:53:39Z")

</div>

> [@Tamas\_Papp](#):
>
> I think the intention of that package is to support “canonical” intervals for each quadrature, which the user can transform easily.

Thanks, Tamas. While it is immediately clear how to transform the quadrature abscissae, it is not clear to me how to transform the quadrature weights without re-evaluating the procedure used to produce the [-1,1] weights.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [August 19, 2020, 3:12pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/9 "2020-08-19T15:12:19Z")

</div>

Unless there is some complication I am not seeing, if quadgk can evaluate

\int\_a^b f(x) w(x) dx

and I need

\int\_c^d g(x) dx

I usually

1. transform from g to the relevant domain using [change of variables](https://en.wikipedia.org/wiki/Integration_by_substitution),

2. take out the w in some nice way from the result.

---

<div class="post-metadata">

**Author:** ![tyleransom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tyleransom/32/530_2.png) [@tyleransom](https://discourse.julialang.org/u/tyleransom)\
**Post date:** [August 19, 2020, 3:20pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/10 "2020-08-19T15:20:12Z")

</div>

I hadn’t thought about it in this way; thanks for sharing.

---

<div class="post-metadata">

**Author:** ![chisquared](https://avatars.discourse-cdn.com/v4/letter/c/87869e/32.png) [@chisquared](https://discourse.julialang.org/u/chisquared)\
**Post date:** [September 20, 2020, 10:10pm UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/11 "2020-09-20T22:10:02Z")

</div>

In case anyone’s reading this and has the same problem, [Expectations.jl](https://github.com/QuantEcon/Expectations.jl) offers a really convenient wrapper for [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl).

**Edit:** I should add that, as the name suggests, Expectations.jl is convenient only for integrating functions with respect to some univariate probability distribution (as in the OP). You can make it work for arbitrary integrals, but that’s probably a lot more complicated than just using FastGaussQuadrature directly.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 14, 2023, 9:48am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/12 "2023-06-14T09:48:25Z")

</div>

> [@ChrisRackauckas](#):
>
> We currently don’t handle derivatives w.r.t. the bounds of the integral though, but that’s something we plan to add soon (also, there’s a slightly better way to compute what we’re doing since we aren’t caching the primal evaluation, but 🤷‍♂️)

Hello, has this been implemented? I can’t still get it working…

```julia
using Integrals,ForwardDiff
k=100; a = 10 # Parameters of the demand function
demand_price(r;k=k,a=a) = k*exp(-a*r)
welfare_manual(r;k=k,a=a) = - k*(exp(-a*r) - 1) / a
function welfare_auto(r::T;k=k,a=a) where {T}
    h(q,unused;k=k,a=a) = demand_price(q,k=k,a=a)
    prob = IntegralProblem(h, zero(T), r)
    sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
    return sol.u
end

welfare_manual(0.05) ≈ welfare_auto(0.05) # ok 
ForwardDiff.derivative(welfare_manual,0.5)
ForwardDiff.derivative(welfare_auto,0.5) # ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(welfare_auto), Float64}, Float64, 1})

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2023, 3:07am UTC](https://discourse.julialang.org/t/array-issue-with-optimzing-a-function-with-quadgk/41837/13 "2023-06-15T03:07:44Z")

</div>

Many derivatives were added but not all. It’s still a bit work in progress on some aspects of differentiating boundaries of higher dimensional integrals.