When Zygote.jl meets Quadgk.jl: Mutating arrays is not supported -- called setindex!(::Vector{QuadGK.Segment})

General Usage
,
1

I have a problem calculating the derivative of a function that involves an integral operation. Examples are as follows:

using QuadGK
using Zygote

function test()
    f(t)=quadgk(x->sin(x)*exp(t*cos(x))  , 0, pi, atol=1e-16)[1]
    @show f(1)
    g=gradient(x -> f(x), 1.0)[1]
    @show g
end
test()

The results are as follows:

f(1) = -4.428138526555854
ERROR: Mutating arrays is not supported -- called setindex!(::Vector{QuadGK.Segment{Float64, Float64, Float64}}, _...)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] (::Zygote.var"#441#442"{Vector{QuadGK.Segment{Float64, Float64, Float64}}})(#unused#::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/lib/array.jl:71
  [3] (::Zygote.var"#2337#back#443"{Zygote.var"#441#442"{Vector{QuadGK.Segment{Float64, Float64, Float64}}}})(Δ::Nothing)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
  [4] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:45 [inlined]
  [5] (::typeof(∂(percolate_up!)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
  [6] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:48 [inlined]
  [7] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:73 [inlined]
  [8] (::typeof(∂(heappush!)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
  [9] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:63 [inlined]
 [10] (::typeof(∂(adapt)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [11] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:35 [inlined]
 [12] (::typeof(∂(do_quadgk)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [13] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:179 [inlined]
 [14] (::typeof(∂(λ)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [15] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:113 [inlined]
 [16] (::typeof(∂(handle_infinities)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [17] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
 [18] (::typeof(∂(#quadgk#27)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [19] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [20] (::Zygote.var"#1750#back#214"{Zygote.var"#212#213"{Tuple{NTuple{9, Nothing}, Tuple{}}, typeof(∂(#quadgk#27))}})(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [21] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
 [22] (::typeof(∂(quadgk##kw)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [23] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [24] (::Zygote.var"#1750#back#214"{Zygote.var"#212#213"{Tuple{Tuple{Nothing, Nothing, Nothing}, Tuple{Nothing, Nothing}}, typeof(∂(quadgk##kw))}})(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [25] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:173 [inlined]
 [26] (::typeof(∂(#quadgk#26)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [27] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [28] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [29] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:173 [inlined]
 [30] (::typeof(∂(quadgk##kw)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [31] Pullback
    @ ./REPL[3]:2 [inlined]
 [32] (::typeof(∂(f)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [33] Pullback
    @ ./REPL[3]:4 [inlined]
 [34] (::Zygote.var"#56#57"{typeof(∂(#2))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface.jl:41
 [35] gradient(f::Function, args::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface.jl:76
 [36] test()
    @ Main ./REPL[3]:4
 [37] top-level scope
    @ REPL[4]:1

Zygote.buffer() doesn’t seem to solve this problem. I want to consult what good idea does everybody have to be able to solve the problem.

2

In your example f(x) is not actually a function of x, but rather a constant as your range of integration is fixed:

\int_0^\pi e^x \cdot \sin x \cdot \cos x \; dx \approx -4.428

Thus, you don’t need to take any derivative.
In case, your integration range depends on x I would suggest adding an adjoint for quadgk which makes use of the fundamental theorem of calculus. This should be way more efficient than trying to autodiff through a numeric integral.

3

QuadGK is not differentiable.

1 Like
4

I think this is what GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical integration for the SciML scientific machine learning organization etc. is intended for, but not sure the exact status of available algorithms. Looks like QuadGK is available, though? Integral Solver Algorithms · Integrals.jl

Regadless, you don’t want to differentiate the quadrature directly and need an adjoint rule implementing the calculus at a higher level.

5

yeah it defines the differentiation overloads. See:

It doesn’t correct differentiation of the integral limits with forward mode right now, but it does do them with reverse. And the derivative w.r.t. parameters is handled in both cases.

1 Like
6

I’m terribly sorry that it is my mistake. I think my function is going to look something like this,

using QuadGK
using Zygote

function test()
    f(t)=quadgk(x->sin(x)*exp(t*cos(x))  , 0, pi, atol=1e-16)[1]
    @show f(1)
    g=gradient(x -> f(x), 1.0)[1]
    @show g
end
test()
  • What I really need is something like this
7

It’s great that Integrals.jl solved my problem. Thank you very much!

1 Like