# Advice on numerical integration of a very "sparse" function

**URL:** <https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176>\
**Category:** Numerics\
**Tags:** question\
**Created:** [November 30, 2018, 10:00am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176 "2018-11-30T10:00:55Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 10:00am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/1 "2018-11-30T10:00:55Z")

</div>

I have a function that when sampled evenly across its two domains looks something like this:

 ![screenshot_2018-11-30_09%3A32%3A24](https://global.discourse-cdn.com/julialang/original/3X/c/3/c378df91b0e346560cc381c36b2e0272693fa7d6.png)  
Almost everything is zero…  
I need to integrate it. I tried [HCubature.jl](https://github.com/stevengj/HCubature.jl) and [Cuba.jl](https://github.com/giordano/Cuba.jl) but they either returned zero or failed (see related [issue](https://github.com/stevengj/HCubature.jl/issues/13)). I’m not looking for some ultra accuracy, but it would be good to get a result that isn’t zero.  
Another approach would be to use [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) some how. While there should be a good way of doing this (see [this](https://github.com/JuliaMath/Interpolations.jl/issues/163) issue), there should also be a sub-optimal but simpler way of integrating this…

Does anyone have any suggestions on how…?

Thanks!

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 30, 2018, 10:10am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/2 "2018-11-30T10:10:55Z")

</div>

ApproxFun might work here:

```julia
@time f = Fun((x,y) -> exp(-1000*((x-0.4)^2+y^2)), (0..1)^2) # 0.05s
sum(f) # 0.0015707963267948962, in 0.0025s

```

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 10:16am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/3 "2018-11-30T10:16:15Z")

</div>

Can’t believe it, I got zero again…! Hmmmmmm

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 30, 2018, 10:19am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/4 "2018-11-30T10:19:06Z")

</div>

Is `norm(f.coefficients) == 0`? If so it’s somehow missed the hump. Dividing your rectangle in two so a corner is at the hump should fix it.

Otherwise, maybe the integral is zero?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2018, 10:22am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/5 "2018-11-30T10:22:24Z")

</div>

Do you know in advance where the bump is?

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [November 30, 2018, 10:32am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/6 "2018-11-30T10:32:47Z")

</div>

It seems a bit too simple, but `sum()` might work over the sampling?

A better estimate might be with a [1 1; 1 1] / 4 filter, estimating box heights from the corners

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 10:49am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/7 "2018-11-30T10:49:10Z")

</div>

> [@dlfivefifty](#):
>
> `norm(f.coefficients) == 0`

was zero, but

> [@dlfivefifty](#):
>
> Dividing your rectangle in two so a corner is at the hump should fix it.

resulted in this:

```julia
julia> @time f = Fun((x,y) -> getsignal([x,y]), (0.82..1)^2)
┌ Warning: Maximum number of coefficients 1048577 reached in constructing Fun.                         
└ @ ApproxFun ~/.julia/packages/ApproxFun/DkMmw/src/Fun/constructors.jl:141
 16.831609 seconds (209.32 M allocations: 5.272 GiB, 13.14% gc time)
Fun(Chebyshev(0.82..1.0)⊗Chebyshev(0.82..1.0),[0.0165691, 0.031849, -0.0267369, 0.0281701, -0.05141, 0.0114935, 0.0226335, -0.0455183, 0.0221598, 0.00361066 … 6.72985e-9, 1.8547e-8, -1.45913e-7, 2.26589e-7, -1.73169e-7, 5.40419e-8, 7.07458e-9, 2.9705e-8, -1.05139e-7, 3.52394e-8])

julia> sum(f)
0.00016545352622070335

```

That sum call was also pretty slow…

> [@rveltz](#):
>
> Do you know in advance where the bump is?

No I do not. But I could uniformly sample about 100 points and then I can get it…

> [@y4lu](#):
>
> but `sum()` might work over the sampling?

This might not work, since I need to compare the results on this integrand to others where the number of hits (as in a location where the function is not equal to zero) will be different and therefore be very susceptible to noise.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 30, 2018, 11:13am UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/8 "2018-11-30T11:13:10Z")

</div>

You’re lucky that your function is only 2d.

So your first question should be: How do you discover the needle in the haystack (regions where the function is far from zero)? An easy approach might be to use IntervalRootFinding.jl for `f-epsilon`. That abuses the fact that you have structural knowledge of `f`:

You don’t just have an oracle for `f` or derivatives; no, you have a machine that computes `f`! This of course only works for some functions (for most reasonable functions you can write down in julia); and it is only fast for certain functions. Once you discovered the needle in the haystack, you could zoom in for HCubature.

Alternatively, you can take a lot of samples. That may miss narrow spikes with large mass; depending on your `f`, one or the other way may be more efficient.

Note that the usefulness of IntervalRootFinding depends not only on your mathematical function (the relation between input and output), but depends crucially on the explicit algorithm used to compute `f`.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 30, 2018, 12:06pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/9 "2018-11-30T12:06:42Z")

</div>

> [@yakir12](#):
>
> I’m not looking for some ultra accuracy, but it would be good to get a result that isn’t zero.

Presumably you are looking for _some_ accuracy, otherwise

```julia
gimmeresult() = 42.0

```

would work (as it isn’t zero :troll:).

Can you provide an MWE that is similar to your `getsignal`?

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 12:12pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/10 "2018-11-30T12:12:56Z")

</div>

> [@Tamas\_Papp](#):
>
> Can you provide an MWE that is similar to your `getsignal` ?

hmm, `getsignal` is the result of ray tracing a scallop eye. It’s two dimensions are elevation and azimuth angles of the light leaving a point source. At larger viewing angles only a fraction of the light enters the eye and ends up as a “signal”. I could set it all up for you to explore…  
Right now I’m testing Interpolations and then some sort of Trapezoidal rule… Maybe Julia has already a `trapz` function?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 30, 2018, 12:18pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/11 "2018-11-30T12:18:38Z")

</div>

If the function is “smooth” (continuous and maybe reasonably differentiable, or a good approximation thereof), you should be much better off with adaptive cubature and/or ApproxFun.jl. A naive trapezoidal rule is kind of a last resort for badly behaved functions, at a large cost of accuracy. The fact that the aforementioned methods fail suggests that exploration using an MWE could be helpful.

Also, for debugging you could try to plot/integrate a 1D slice, eg one along the y axis that contains the “bump”.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 12:24pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/12 "2018-11-30T12:24:47Z")

</div>

It’s smooth where it’s \> 0, everywhere else it’s zero. I think the bumps might even have just one maxima… But I’m not sure. It just irritates me that while I can quickly see where all the bumps are by evaluating 100x100 points, neither of the Cuba-related packages pick those up. Cause then it should be fine. Finding the needle first, like @foobar_lv2 suggests, might be most robust, but to me, it seems like an overkill a bit. Might be wrong.  
I’ll try and set it up then…!

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 30, 2018, 12:27pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/13 "2018-11-30T12:27:03Z")

</div>

If you know the bounds, perhaps you can brute force an approximate integral using Sobol sequences,

[https://github.com/stevengj/Sobol.jl](https://github.com/stevengj/Sobol.jl)

which should also be a reasonable quick & dirty solution for finding the (proximity of the) bump.

---

<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:** [November 30, 2018, 1:04pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/14 "2018-11-30T13:04:59Z")

</div>

> [@yakir12](#):
>
> It just irritates me that while I can quickly see where all the bumps are by evaluating 100x100 points

We can add some keyword argument to force HCubature to subdivide the domain at the beginning to ensure that it is sampled more finely to start with.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 1:11pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/15 "2018-11-30T13:11:00Z")

</div>

I’m testing Cuba.jl’s xgiven to see if that helps (from [giordano’s reply](https://github.com/stevengj/HCubature.jl/issues/13#issuecomment-443156747) )…  
Also, trying Sobol.jl, cool!

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [November 30, 2018, 1:39pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/16 "2018-11-30T13:39:00Z")

</div>

Maybe try some more options with cuba? Here is an example where cuba finds even a much smaller 2d peak out of the box:

```julia
using Cuba

function f!(v, out)
    # tiny circle of radius r
    # expected integral is 1
    r = 1e-3
    x,y = v
    if (x-0.5)^2 + (y-0.5)^2 < r^2
        out[1] = 1 / (pi * r^2)
    else
        out[1] = 0
    end
end

@show vegas(f!,2)

```

```julia
vegas(f!, 2) = Component:
 1: 0.920180950098721 ± 0.0004369999395848204 (prob.: 1.0)
Integrand evaluations: 1007500
Fail: 1
Number of subregions: 0

```

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 1:45pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/17 "2018-11-30T13:45:22Z")

</div>

🤷‍♂️

```julia
julia> vegas(f!,2)
Component:                                                                                             
 1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Fail: 0
Number of subregions: 0

```

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 30, 2018, 1:56pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/18 "2018-11-30T13:56:26Z")

</div>

Your function was too easy: 0.5 is an obvious sample point. E.g.

```julia
julia> function f!(v, out)
           # tiny circle of radius r
           # expected integral is 1
           r = 1e-3
           xm = 0.7654
           ym = 0.4567
           x,y = v
           if (x-xm)^2 + (y-ym)^2 < r^2
               out[1] = 1 / (pi * r^2)
           else
               out[1] = 0
           end
       end
julia> @show vegas(f!, 2)
vegas(f!, 2) = Component:
 1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Fail: 0
Number of subregions: 0
Component:
 1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Fail: 0
Number of subregions: 0

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 30, 2018, 2:04pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/19 "2018-11-30T14:04:11Z")

</div>

This is a nice example: for an arbitrary numerical integration method, one can choose a `(xm, ym)` and `r` small enough so that this will fail. For a smoothly differentiable function, one could find the peak numerically (especially after a meaningful transformation) and work around that, or use MCMC, but not for this one.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [November 30, 2018, 2:08pm UTC](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176/20 "2018-11-30T14:08:30Z")

</div>

So I used `Sobol` to improve on the sampling, and with just 25 points I already hit all the live areas. I then feed those into `Cuba` via the `xgiven` optional argument of `divonne` and it works! Pretty promising…

[Next page](https://discourse.julialang.org/t/advice-on-numerical-integration-of-a-very-sparse-function/18176.md?page=2)