# "vectorized" QuadGK?

**URL:** <https://discourse.julialang.org/t/vectorized-quadgk/39853>\
**Category:** Performance\
**Tags:** question\
**Created:** [May 20, 2020, 7:51pm UTC](https://discourse.julialang.org/t/vectorized-quadgk/39853 "2020-05-20T19:51:20Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [May 20, 2020, 7:51pm UTC](https://discourse.julialang.org/t/vectorized-quadgk/39853/1 "2020-05-20T19:51:20Z")

</div>

I’ll start off by saying that [QuadGK](https://github.com/JuliaMath/QuadGK.jl) is really fast. But I want to know how to make it faster for what I think is a relatively common usecase, which is its vectorized form. Broadcasting the `quadgk` call itself performs as we’d expect (slightly slower than two times a single call), whereas if broadcasting the function that we are integrating (so that the integrand is vector valued) results in a slowdown.

In this specific case, there should be speedups available because although the integrand is vector valued it is basically copies of similar integrals which differ only by the parameter `a` (so that it seems like `quadgk` could reuse some computations).

```julia
julia> using QuadGK, BenchmarkTools

julia> f(x,a) = exp(-a * x^2)
f (generic function with 1 method)

julia> quad1(a) = quadgk(x ->f(x,a), 0, 1, rtol=1e-12)
quad1 (generic function with 1 method)

julia> quadvec(a) = quadgk(x ->f.(x,a), 0, 1, rtol=1e-12)
quadvec (generic function with 1 method)

julia> avec = [0.5; 1.0 ;1.5]
3-element Array{Float64,1}:
 0.5
 1.0
 1.5

julia> @btime quad1($1.0)
  976.333 ns (39 allocations: 992 bytes)
(0.746824132812427, 0.0)

julia> # broadcasting quadgk
       @btime quad1.($avec)
  2.329 μs (87 allocations: 2.25 KiB)
3-element Array{Tuple{Float64,Float64},1}:
 (0.8556243918921488, 2.55351295663786e-15)
 (0.746824132812427, 0.0)
 (0.6633509458403348, 1.1102230246251565e-16)

julia> # broadcasting f
       @btime quadvec($avec)
  5.054 μs (157 allocations: 16.36 KiB)
([0.8556243918921488, 0.746824132812427, 0.6633509458403348], 1.3401577416544657e-16)

```

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [May 20, 2020, 8:08pm UTC](https://discourse.julialang.org/t/vectorized-quadgk/39853/2 "2020-05-20T20:08:42Z")

</div>

Hmm ok I guess that it does just fine! It just depends on how many integrals it needs to compute.

```julia
julia> avecB = rand(1000)
1000-element Array{Float64,1}:
 0.9109829394262448
 0.7401979941321699
 0.004438422199036562
 0.7622184357652708
 0.3752119597124437
 0.0870799867237182
 0.7612865115744936
 0.8395578809491078
 ⋮
 0.48755227092235875
 0.06770987500160586
 0.029118998270370566
 0.8758783271579778
 0.1163149703976627
 0.9930090806821483
 0.09210918209385932

julia> # broadcasting quadgk
       @btime quad1.($avecB)
  217.334 μs (7171 allocations: 191.77 KiB)

julia> # broadcasting quadgk
       @btime quadvec($avecB)
  106.152 μs (50 allocations: 373.16 KiB)

```

Would still appreciate suggestions for improvements here, but less necessary.

---

<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:** [May 20, 2020, 8:49pm UTC](https://discourse.julialang.org/t/vectorized-quadgk/39853/3 "2020-05-20T20:49:15Z")

</div>

> [@tbeason](#):
>
> `quadvec(a) = quadgk(x ->f.(x,a), 0, 1, rtol=1e-12)`

If you want this to be fast for small vectors `a` you should make `a` into an `SVector` (from the StaticArrays package), so that all the loops are unrolled and you aren’t allocating zillions of tiny arrays on the heap.

(If `a` is a large vector, the overhead of heap-allocating temporary arrays etcetera become negligible.)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [May 20, 2020, 9:17pm UTC](https://discourse.julialang.org/t/vectorized-quadgk/39853/4 "2020-05-20T21:17:46Z")

</div>

Great! That completely switches the conclusion for small `n` and even shaves off some time from the `n=1000` case. I don’t know yet what my problem size will be, but will have to keep this in mind.