# Performance problems in one-dimensional complex integration

**URL:** <https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912>\
**Category:** Numerics\
**Tags:** question, performance, integral\
**Created:** [June 27, 2023, 7:01pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912 "2023-06-27T19:01:44Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![javier\_mj](https://avatars.discourse-cdn.com/v4/letter/j/e56c9b/32.png) [@javier\_mj](https://discourse.julialang.org/u/javier_mj)\
**Post date:** [June 27, 2023, 7:01pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/1 "2023-06-27T19:01:45Z")

</div>

Hello to all  
I am currently moving some code to Julia with the intention of improving performance. However, within a function I am running into a bottleneck when calculating the integral in one dimension of a complex function in the range [0, 1000].

This function is as follows

```julia
using HCubature

S₀ = 4000.0
K = 100.0
V₀ = 0.83500
r = 0.03
κ = 2.39271
θ = 0.01887
σ = 1.48821
λ = -0.38643
ρ = -0.97271
τ = 1.0

function price_algorithm(S₀::Float64, V₀::Float64, κ::Float64, θ::Float64, 
    σ::Float64, ρ::Float64, λ::Float64, K::Float64, τ::Float64, r::Float64)

    a = κ * θ
    b = κ + λ

    function price_function(ϕ::T) where {T<:Number}

        c = ρ * σ * ϕ * 1im

        d = sqrt((c - b)^2 + (ϕ * 1im + ϕ^2) * σ^2)

        g = (b - c + d) / (b - c - d)

        part1 = exp(r * ϕ * 1im * τ)
        part2 = S₀^(ϕ * 1im) * ((1 - g * exp(d * τ)) / (1 - g))^(-2 * a / σ^2)
        part3 = exp(a * τ * (b - c + d) / σ^2 + V₀ * (b - c + d) * ((1 - exp(d * τ)) / (1 - g * exp(d * τ))) / σ^2)

        return part1 * part2 * part3

    end

    function price_integral(ϕ::Float64)

        num_integral = exp(r * τ) * price_function(ϕ - 1im) - K * price_function(ϕ)
    
        dem_integral = 1im * ϕ * K^(1im * ϕ)
    
        return num_integral / dem_integral
    
    end

    integral, _ = hquadrature(price_integral, 0, 1000)

    return (S₀ - K * exp(-r * τ)) / 2 + real(integral) / π

end

```

The performance measured through `@benchmark` is.

```julia
julia> using BenchmarkTools

julia> @benchmark price_algorithm($S₀, $V₀, $κ, $θ, $σ, $ρ, $λ, $K, $τ, $r)

BenchmarkTools.Trial: 110 samples with 1 evaluation.
 Range (min … max): 39.698 ms … 90.912 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 42.571 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 45.532 ms ± 7.154 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▄▁█ ▃         
  ▇███▄█▃▅▅▃▃▃▄▆▄▃▄▄▆▁▄▄▄▅▁▄▄▁▄▁▃▁▄▃▁▃▁▁▁▁▄▁▃▁▁▁▁▁▁▁▃▁▁▁▁▁▃▁▄ ▃
  39.7 ms Histogram: frequency by time 64.5 ms <

 Memory estimate: 129.31 KiB, allocs estimate: 6.

```

Which seems to me to be extremely slow.

I would appreciate any recommendations on how to integrate this function in a more efficient way.

**Note** : Although I know the `QuadGK` package and its function to integrate in one dimension `quadgk`, this function gives a domain error for very small values close to zero (for example `1.424047472694446089e-303`).

---

<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:** [June 27, 2023, 8:02pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/2 "2023-06-27T20:02:18Z")

</div>

> [@javier\_mj](#):
>
> Although I know the `QuadGK` package and its function to integrate in one dimension `quadgk`, this function gives a domain error for very small values close to zero (for example

~~Probably that’s because you aren’t passing an absolute tolerance.~~ Or any tolerance at all, for that matter. The default is a _relative_ tolerance of `√ε` (about 8 digits), which may be more accurate than you need, and is definitely more accurate than you need if the integral mostly cancels.

_Update:_ Actually it’s because the imaginary part of your integral is blowing up for your test values. If I add a `@show` before the `integral, _ = ...` in your posted code, I get:

```julia
(integral, _) = hquadrature(price_integral, 0, 1000, atol = 0.001) = (6453.582138997572 - Inf*im, Inf)

```

So, you are potentially timing garbage, and `quadgk` is complaining about the divergence.

---

<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:** [June 27, 2023, 8:25pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/3 "2023-06-27T20:25:02Z")

</div>

> [@stevengj](#):
>
> Actually it’s because your integral is blowing up for your test values.

One useful thing I would always recommend when you are investigating this sort of thing is to plot your integrand. Here, the basic problem is that the imaginary part of your integrand blows up:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/0/b08d99178ca9474d387ee68f64afef7f9223cdc6.png)  
Since you only need the real part of the integral at the end, I would suggest taking the real part _inside_ the integral instead of afterwards — then you don’t waste effort trying to integrate the divergent term that you discard anyway.

Furthermore, rather than integrating from 0 to 1000, which I’m guessing is an approximation for integrating from 0 to \infty, I would just integrate directly [from `0` to `Inf` in QuadGK](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Improper-integrals:-Infinite-limits), which should typically be much more efficient.

If you are still getting NaNs, you’ll want to carefully audit your integrand function to figure out what’s going on, because you may be having floating-point overflow or similar issues that will affect correctness (but can generally be eliminated by careful rewriting).

---

<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:** [June 27, 2023, 8:37pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/4 "2023-06-27T20:37:54Z")

</div>

> [@stevengj](#):
>
> Furthermore, rather than integrating from 0 to 1000, which I’m guessing is an approximation for integrating from 0 to \infty∞\infty, I would just integrate directly [from `0` to `Inf` in QuadGK](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Improper-integrals:-Infinite-limits), which should be much more efficient.
> 
> If you are still getting NaNs, you’ll want to carefully audit your integrand function to figure out what’s going on, because you may be having floating-point overflow or similar issues that will affect correctness (but can generally be eliminated by careful rewriting).

In particular, I noticed that your exponentials are overflowing when `real(d * τ)` gets large, producing `NaN`, but in this regime the integrand should be exponentially small (magnitude \< 10^{-300}) and hence you can just return zero. So I added a check:

```julia
real(d * τ) > 700 && return zero(c)

```

after the `d = ...` line in `price_function`. Then I changed your integration line to:

```julia
integral, _ = quadgk(real ∘ price_integral, 0, Inf)

```

The computation time (as reported by `@btime`) sped up by 7x on my machine.

This is still using the default tolerance of `rtol = sqrt(eps())`. If I change it to pass `rtol = 1e-3` (3 significant digits), then I get another 8x speedup (down to ≈ 0.3 ms).

---

<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:** [June 27, 2023, 8:56pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/5 "2023-06-27T20:56:23Z")

</div>

> [@javier\_mj](#):
>
> ```julia
> price_algorithm(S₀::Float64, V₀::Float64, κ::Float64, θ::Float64, 
> σ::Float64, ρ::Float64, λ::Float64, K::Float64, τ::Float64, r::Float64)
> 
> ```

By the way, all of these type declarations accomplish nothing for performance. They just make your code less readable and less generic. Read the [manual section on “Argument-Type Declarations”](https://docs.julialang.org/en/v1/manual/functions/#Argument-type-declarations).

---

<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:** [June 27, 2023, 9:00pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/6 "2023-06-27T21:00:14Z")

</div>

There are various additional optimizations that can help to speed up your integrand. e.g. a simple one is to use `cis(x)` for e^{ix}, e.g. I got a 20% speedup by using the following in `price_function` instead of generic exponentiation:

```julia
part1 = cis(r * ϕ * τ)
part2 = cis(logS₀ * ϕ) * ((1 - g * exp(d * τ)) / (1 - g))^(-2 * a / σ^2)

```

where `logS₀ = log(S₀)`.

To do significantly better you might have to do more analysis of your integrand. e.g. if you can analytically work out the asymptotic decay rate, and can rescale it to e^{-x}, then you can probably use a [Gauss–Laguerre quadrature scheme](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/#Gauss-Laguerre-quadrature) to good effect.

---

<div class="post-metadata">

**Author:** ![javier\_mj](https://avatars.discourse-cdn.com/v4/letter/j/e56c9b/32.png) [@javier\_mj](https://discourse.julialang.org/u/javier_mj)\
**Post date:** [June 27, 2023, 9:23pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/7 "2023-06-27T21:23:32Z")

</div>

Thank you very much for your advice, I have just noticed an 100x performance improvement over the original code.  
I will try to do more analysis of my integrated to further improve the results.

```julia
@btime price_algorithm($S₀, $V₀, $κ, $θ, $σ, $ρ, $λ, $K, $τ, $r)
379.400 μs (3 allocations: 1.69 KiB)

```

---

<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:** [June 27, 2023, 9:26pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/8 "2023-06-27T21:26:17Z")

</div>

> [@javier\_mj](#):
>
> I have just noticed an 8x performance improvement over the original code.

Isn’t 379µs a 100x improvement over your original code (40ms), not 8x?

---

<div class="post-metadata">

**Author:** ![javier\_mj](https://avatars.discourse-cdn.com/v4/letter/j/e56c9b/32.png) [@javier\_mj](https://discourse.julialang.org/u/javier_mj)\
**Post date:** [June 27, 2023, 9:26pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/9 "2023-06-27T21:26:57Z")

</div>

Sorry, I misspelled the number

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 27, 2023, 10:33pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/10 "2023-06-27T22:33:20Z")

</div>

Some more issues:

1. `exp(d * τ)` repeats multiple times in the same function. I don’t know whether Julia is able to fix this via common subexpression elimination, however in any case it would make sense to deduplicate your code.

2. Your functions `price_function` and `price_integral` capture many variables from `price_algorithm`, and you do it in such a way that (still AFAIK) causes huge performance problems in Julia. [This](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured) is the relevant section in the manual.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 27, 2023, 10:35pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/11 "2023-06-27T22:35:43Z")

</div>

Julia isn’t yet smart enough to remove the redundant `exp(d * τ)` calls. We have the analysis (the conditions are exactly the same as dead code elimination), but no one has written the pass that will remove the redundancy. [Rebase of julia effects to LLVM by gbaraldi · Pull Request #50188 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/50188) probably would let LLVM remove the redundant call though.

---

<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:** [June 27, 2023, 10:43pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/12 "2023-06-27T22:43:41Z")

</div>

> [@Oscar\_Smith](#):
>
> probably would let LLVM remove the redundant call though.

Probably not for complex arguments?

---

<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:** [June 27, 2023, 10:44pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/13 "2023-06-27T22:44:28Z")

</div>

> [@nsajko](#):
>
> Your functions `price_function` and `price_integral` capture many variables from `price_algorithm`, and you do it in such a way that (still AFAIK) causes huge performance problems in Julia. [This](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured) is the relevant section in the manual.

I don’t think that’s an issue here? That would show up as lots of additional allocations in `@btime`, which aren’t happening here.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 27, 2023, 10:45pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/14 "2023-06-27T22:45:31Z")

</div>

On the Julia side the effects are good

```julia
julia> Base.infer_effects(exp, (Complex{Float64},))
(+c,+e,!n,+t,+s,+m,+i)

```

So if we taught LLVM about the Julia effects which that PR does, LLVM should be able to remove the call.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 27, 2023, 10:46pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/15 "2023-06-27T22:46:10Z")

</div>

I guess that the compiler sees that the captures are read-only so it doesn’t box them? Nice!

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 27, 2023, 10:52pm UTC](https://discourse.julialang.org/t/performance-problems-in-one-dimensional-complex-integration/100912/16 "2023-06-27T22:52:23Z")

</div>

> [@javier\_mj](#):
>
> ```julia
> function price_integral(ϕ::Float64)
> 
> num_integral = exp(r * τ) * 
> 
> ```

As far as I see you could move the `exp(r * τ)` to outside `price_integral`.