# Does integral function related optimization supported in Julia?

**URL:** <https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836>\
**Category:** Optimization (Mathematical)\
**Tags:** integral-function, nonlinear-optimizati\
**Created:** [February 18, 2023, 1:55pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836 "2023-02-18T13:55:51Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 18, 2023, 1:55pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/1 "2023-02-18T13:55:51Z")

</div>

Hi all,  
I have a nonlinear optimization problem that involves integral functions, for instance:

`F(x, a) = $$\int_0^x a*(x - 4y^3)dy + \int_0^x\int_0^3y (7s^2)dsdy$$`

How can I define this objective (involves single and double integral) function for Julia optimization packages like JuMP, Optim.jl etc.?

Thank you very much!

---

<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:** [February 18, 2023, 3:41pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/2 "2023-02-18T15:41:47Z")

</div>

> [@Jian\_ZUO](#):
>
> How can I define this objective (involves single and double integral) function for Julia optimization packages like JuMP, Optim.jl etc.?

Use Integrals.jl

> **[GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical...](https://github.com/SciML/Integrals.jl)**
>
> A common interface for quadrature and numerical integration for the SciML scientific machine learning organization - GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical int...

It adds what’s required for automatic differentation.

> **[Differentiating Integrals · Integrals.jl](https://docs.sciml.ai/Integrals/stable/tutorials/differentiating_integrals/)**
>
> Documentation for Integrals.jl.

Then use it with Optimization.jl

> **[Optimization.jl: A Unified Optimization Package · Optimization.jl](https://docs.sciml.ai/Optimization/stable/)**
>
> Documentation for Optimization.jl.

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 18, 2023, 3:56pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/3 "2023-02-18T15:56:42Z")

</div>

Thank you Chris for the suggestions.  
Sorry I am a new user of Integrals.jl. I have checked the differential part and tried to define the integral function. Note in my case I need to define the integral as a function of its upper (or lower) limits. For instance:

```julia
using Integrals
f(x, p) = @. x^3 + 2x*p
# p = ones(1)
function intf(p)
    lb = zeros(1)
    ub = p
    prob = IntegralProblem(f, lb, ub, p)
    solve(prob)[1]
end
intf(1)

```

This is not working. Can you help to clarify the possible way to define it in Integrals.jl?  
Thanks!

---

<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:** [February 18, 2023, 4:30pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/4 "2023-02-18T16:30:29Z")

</div>

> [@Jian\_ZUO](#):
>
> This is not working. Can you help to clarify the possible way to define it in Integrals.jl?

What error did you get?

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 18, 2023, 4:32pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/5 "2023-02-18T16:32:47Z")

</div>

```julia
MethodError: no method matching hcubature_(::Integrals.var"#34#36"{IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, ::Vector{Float64}, ::Int64, ::typeof(LinearAlgebra.norm), ::Float64, ::Float64, ::Int64, ::Int64)
Closest candidates are:
  hcubature_(::F, ::AbstractVector{T}, !Matched::AbstractVector{S}, ::Any, ::Any, ::Any, ::Any, ::Any) where {F, T<:Real, S<:Real} at ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:125
  hcubature_(::F, !Matched::StaticArraysCore.SVector{n, T}, !Matched::StaticArraysCore.SVector{n, T}, ::Any, ::Any, ::Any, ::Any, ::Any) where {F, n, T<:Real} at ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:49
  hcubature_(::Any, !Matched::Tuple{Vararg{Real, n}}, !Matched::Tuple{Vararg{Real, n}}, ::Any, ::Any, ::Any, ::Any, ::Any) where n at ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:131

Stacktrace:
  [1] hcubature(f::Integrals.var"#34#36"{IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, a::Vector{Float64}, b::Int64; norm::Function, rtol::Float64, atol::Float64, maxevals::Int64, initdiv::Int64)
    @ HCubature ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:179
  [2] __solvebp_call(prob::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::HCubatureJL{typeof(LinearAlgebra.norm)}, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Vector{Float64}, ub::Int64, p::Int64; reltol::Float64, abstol::Float64, maxiters::Int64)
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:157
  [3] __solvebp_call(prob::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::HCubatureJL{typeof(LinearAlgebra.norm)}, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Vector{Float64}, ub::Int64, p::Int64)
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:139
  [4] __solvebp(::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:122
  [5] __solvebp(::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::HCubatureJL{typeof(LinearAlgebra.norm)}, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::Vector{Float64}, ::Int64, ::Int64)
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:122
  [6] solve(prob::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::HCubatureJL{typeof(LinearAlgebra.norm)}; sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, do_inf_transformation::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:108
  [7] #solve#26
    @ ~/.julia/packages/Integrals/leYds/src/Integrals.jl:83 [inlined]
  [8] solve(prob::IntegralProblem{false, Int64, typeof(f), Vector{Float64}, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    @ Integrals ~/.julia/packages/Integrals/leYds/src/Integrals.jl:60
  [9] intf(p::Int64)

```

---

<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:** [February 19, 2023, 8:33pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/6 "2023-02-19T20:33:04Z")

</div>

```julia
using Integrals
f(x, p) = @. x^3 + 2x*p
# p = ones(1)
function intf(p)
    prob = IntegralProblem(f, 0, p, p)
    solve(prob, QuadGKJL())[1]
end
intf(1)

```

The error was saying that you needed to supply an algorithm:

```julia
ERROR: ArgumentError: No integration algorithm `alg` was supplied as the second positional argument.
Reccomended integration algorithms are:
For scalar functions: QuadGKJL()
For ≤ 8 dimensional vector functions: HCubatureJL()
For > 8 dimensional vector functions: MonteCarloIntegration.vegas(f, st, en, kwargs...)
See the docstrings of the different algorithms for more detail.

```

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 19, 2023, 9:09pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/7 "2023-02-19T21:09:32Z")

</div>

Yes, now it works for calculating the function value! Thank you!  
Sorry, I realized that there is also an error in calculating the derivative:

```julia
using Integrals
using ForwardDiff
f(x, p) = @. x^3 + 2x*p
# p = ones(1)
function intf(p)
    prob = IntegralProblem(f, 0, p, p)
    return solve(prob, QuadGKJL(), reltol = 1e-6, abstol = 1e-6)[1]
end
ForwardDiff.derivative(intf, 2.0)

```

I get:

```julia
Output exceeds the [size limit](command:workbench.action.openSettings?[). Open the full output data [in a text editor](command:workbench.action.openLargeOutput?f045d3a8-b688-4854-9a8b-e9a78ff3fb4d)

MethodError: no method matching kronrod(::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}}, ::Int64) Closest candidates are: kronrod(!Matched::Type{T}, ::Integer) where T<:AbstractFloat at ~/.julia/packages/QuadGK/BYxcx/src/gausskronrod.jl:150 Stacktrace: [1] macro expansion @ ~/.julia/packages/QuadGK/BYxcx/src/gausskronrod.jl:259 [inlined] [2] _cachedrule @ ~/.julia/packages/QuadGK/BYxcx/src/gausskronrod.jl:259 [inlined] [3] cachedrule @ ~/.julia/packages/QuadGK/BYxcx/src/gausskronrod.jl:264 [inlined] [4] do_quadgk(f::Integrals.var"#30#31"{IntegralProblem{false, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, typeof(f), Int64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}}, s::Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing) @ QuadGK ~/.julia/packages/QuadGK/BYxcx/src/adapt.jl:7 [5] #46 @ ~/.julia/packages/QuadGK/BYxcx/src/adapt.jl:219 [inlined] [6] handle_infinities(workfunc::QuadGK.var"#46#47"{Float64, Float64, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::Integrals.var"#30#31"{IntegralProblem{false, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, typeof(f), Int64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}}, s::Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(intf), Float64}, Float64, 1}}) @ QuadGK ~/.julia/packages/QuadGK/BYxcx/src/adapt.jl:118 [7] #quadgk#45 @ ~/.julia/packages/QuadGK/BYxcx/src/adapt.jl:218 [inlined] [8] #quadgk#44 @ ~/.julia/packages/QuadGK/BYxcx/src/adapt.jl:213 [inlined] [9] #__solvebp_call#29 @ ~/.julia/packages/Integrals/leYds/src/Integrals.jl:134 [inlined] [10] #__solvebp#28 @ ~/.julia/packages/Integrals/leYds/src/Integrals.jl:122 [inlined]

...

@ Main ~/Julia_projects/Maintenance_EMS/RUL_dist.ipynb:7 [13] derivative(f::typeof(intf), x::Float64) @ ForwardDiff ~/.julia/packages/ForwardDiff/QdStj/src/derivative.jl:14 [14] top-level scope

```

---

<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:** [February 19, 2023, 10:01pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/8 "2023-02-19T22:01:54Z")

</div>

Interesting. Your code in theory should be fine, or the following:

```julia
using Integrals, ForwardDiff
f(x, p) = x[1]^3 + 2x[1]*p[1]
# p = ones(1)
function intf(p)
    prob = IntegralProblem(f, [0.0], [p], [p])
    solve(prob, HCubatureJL())[1]
end
intf(1)
ForwardDiff.derivative(intf, 1.0)

```

@ArnoStrouwen can you look into why the derivative rule didn’t catch these cases?

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 19, 2023, 11:51pm UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/9 "2023-02-19T23:51:24Z")

</div>

Also, I am not sure about the rules in Integrals.jl for defining (the initial formula in my question):  
 ![intg](https://global.discourse-cdn.com/julialang/original/3X/9/7/9797c5bb23866a94390b69443f408e6cce56b747.png)  
I try to code it as follows:

```julia
f1(x, y, p) = p*(y-4*x^3)
f2(x, y, p) = (y*7*x^2)
function F(x, p)
    prob1 = IntegralProblem(f1, [0], [x], p)
    prob2 = IntegralProblem(f1, f2, [0, 0], [3, x])
    val1 = solve(prob, QuadGKJL(), reltol = 1e-6, abstol = 1e-6)[1]
    val2 = solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1]
    return val1 + val2
end

```

What is the correct way to define such a function? Thanks!

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [February 21, 2023, 10:16am UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/10 "2023-02-21T10:16:16Z")

</div>

Wouldn’t it be simpler to compute the analytical value of F(x, p)?

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 21, 2023, 10:25am UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/11 "2023-02-21T10:25:09Z")

</div>

Hi, thanks for the comment.  
The formula I put here is as an example, in my case, the objective is a complicated integral function that is hard to tackle analytically.

---

<div class="post-metadata">

**Author:** ![ArnoStrouwen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnostrouwen/32/8699_2.png) [@ArnoStrouwen](https://discourse.julialang.org/u/ArnoStrouwen)\
**Post date:** [February 22, 2023, 1:12am UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/12 "2023-02-22T01:12:40Z")

</div>

Currently, there are no automatic differentiation rules dealing with sensitivities of the integration domain bounds. I’ll try to get that fixed before next week.  
In the meantime, you can use this stopgap solution:  
I think from calculus, the following holds:  
\frac{d}{dp}\int\_{h(p)}^{g(p)}f(x,p)dx = \int\_{h(p)}^{g(p)}\frac{d}{dp}f(x,p)dx + f(g(p),p)\frac{d}{dp}g(p) - f(h(p),p)\frac{d}{dp}h(p)  
Thus, for the example @ChrisRackauckas posted, you would get:

```julia
using Integrals
using FiniteDiff
using ForwardDiff

using Integrals, ForwardDiff
f(x, p) = x[1]^3 + 2x[1]*p[1]
# p = ones(1)
function intf(p)
    prob = IntegralProblem(f, [0.0], [p], [p])
    solve(prob, HCubatureJL())[1]
end
p = 1.0
intf(p)
FiniteDiff.finite_difference_derivative(intf, p) #4

function intf_stopgap(p)
    prob = IntegralProblem(f, [0.0], [ForwardDiff.value(p)], [p])
    solve(prob, HCubatureJL())[1]
end
lb(p) = 0.0
ub(p) = p

ForwardDiff.derivative(intf_stopgap,p) +
    f([ub(p)],[p])*ForwardDiff.derivative(ub,p) -
    f([lb(p)],[p])*ForwardDiff.derivative(lb,p)# 4

```

The trick here is to use `ForwardDiff.value` to remove the derivative parts of the dual number of the bounds, and to manually add what the bounds contribute by the chain rule.

Another issue, even when you are dealing with scalar values, it is best to pass both bounds and parameters as one element vectors, like in the above example, as the AD for scalars could also use some work.

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 22, 2023, 7:22am UTC](https://discourse.julialang.org/t/does-integral-function-related-optimization-supported-in-julia/94836/13 "2023-02-22T07:22:11Z")

</div>

Thank you very much for pointing out these!  
I am going to try the stopgap solution and keep following the progress  
of Integrals.jl!  
Thanks again for the help!