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!
Use Integrals.jl
It adds what’s required for automatic differentation.
Then use it with Optimization.jl
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:
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!
What error did you get?
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)
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:
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.
Yes, now it works for calculating the function value! Thank you!
Sorry, I realized that there is also an error in calculating the derivative:
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:
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
Interesting. Your code in theory should be fine, or the following:
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?
Also, I am not sure about the rules in Integrals.jl for defining (the initial formula in my question):
I try to code it as follows:
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!
Wouldn’t it be simpler to compute the analytical value of F(x, p)?
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.
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:
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.
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!