# Array issue with optimzing a function with quadgk

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?

using Distributions
using Random
using ForwardDiff
using Zygote
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

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

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


out = alpha * intdiffcdF1 + intdiffcdF2 = 0.03361344537815181
socialPlannerObjective([1.3, 0.3], [10, linearCost, uniform_dist]) = 0.03361344537815181
out = alpha * intdiffcdF1 + intdiffcdF2 = 0.40183246073298395
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::T) where T at essentials.jl:168

Stacktrace:
[4] percolate_down! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:28 [inlined]
[6] percolate_down! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:18 [inlined]
[8] heappop! at /Users/svass/.julia/packages/DataStructures/GvsTk/src/heaps/arrays_as_heaps.jl:59 [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
[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
[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
[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
[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



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).

In particular, if your interval is finite, you should use Gauss-Legendre with fixed nodes/weights, and transform according to θ₁, θ₂ etc. Then integration is simply a dot product, which is differentiable.

5 Likes

Do you have any tips on how to do this?

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 is an implementation (taken from MATLAB’s file exchange).

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

1 Like

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) 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, so that it integrates with Zygote etcetera, would be great.

4 Likes

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)
sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
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 )

3 Likes

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.

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,

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

1 Like