hcubature does not converge on discontinuous functions

[antoine-levitt]

As discussed in giordano/Cuba.jl#16, hcubature seems to fail for a discontinuous integrand (see MWE there). The integrand is the indicator function of a polygon, so a Riemann sum with NxN points should be accurate to 1/N, and the integral should be about 0.88408, with a small error on the last digit. When running Cuba vs Cubature for decreasing values of abstol (reltol and maxiters are fixed at 0), I get:

0.8840847072116771 ± 9.995513758054815e-5 (Cuba.jl)
0.8848895327256797 ± 9.998142861595564e-5 (Cubature.jl)

0.8840898028407925 ± 9.998969682594401e-6 (Cuba.jl)
0.8848941514595511 ± 9.999535960807015e-6 (Cubature.jl)

0.8840867311005485 ± 9.99974855348049e-7 (Cuba.jl)
0.8848940796576487 ± 9.999981972972212e-7 (Cubature.jl)

0.8840865507480177 ± 9.999859751764881e-8 (Cuba.jl)
0.8848940808860499 ± 9.999994401879008e-8 (Cubature.jl)

So Cuba looks about right, but Cubature appears to fail in that case.

[stevengj]

Related to #25 and JuliaMath/HCubature.jl#4, hopefully with the same fix?

[stevengj]

It would be good to have a simpler test case for this problem that does not involve external packages.

(That would be needed to include the test case in runtests.jl once it is fixed, and also makes it easier for me to run since e.g. the Interpolations package doesn't work yet in Julia 0.7.)

[antoine-levitt]

Unfortunately I cannot reproduce easily without Interpolations. I tried with HCubature, but HCubature master is at 0.7. Let's revisit once Interpolations is 0.7 compatible.

[gideonsimpson]

I'm struggling with something that might be related to this. It might be that this method won't work that well for discontinuous functions, but consider the problem of integrating

f = x-> Real.(x[1]+x[2]<1)

over $[0,1]^2$ (which should give an answer of 1/2). I get utterly abysmal performance (i.e., taking minutes to compute)

hcubature(f, [0.,0.], [1.,1.])

[giordano]

For the record

julia> using Cuba

julia> @time cuhre( (x,f) -> f[] = Real.(x[1]+x[2]<1))
  0.177037 seconds (992.96 k allocations: 47.165 MiB, 36.74% gc time)
Component:
 1: 0.49997978547152266 ± 4.996317974744375e-5 (prob.: 0.0)
Integrand evaluations: 170885
Number of subregions:  1315
Note: The desired accuracy was reached

[stevengj]

As I said above, this is fixed in the HCubature package:

julia> using HCubature

julia> count = 0
0

julia> hcubature([0,0],[1,1], rtol=1e-4, atol=1e-12) do x
          global count += 1
          x[1]+x[2]<1
       end
(0.49999990377497977, 4.9994842542085703e-5)

julia> count
214659

which is a number of evaluations within 25% of Cuba.cuhre (but > 2 orders of magnitude more accurate).

(Cuba.cuhre defaults to rtol=1e-4, atol=1e-12, which is much larger than the default tolerance in Cubature and HCubature, so you have to set the tolerance for a fair comparison.)

[stevengj]

Actually, that test problem is fine in the Cubature package too:

julia> using Cubature

julia> count = 0

julia> hcubature([0,0],[1,1], reltol=1e-4, abstol=1e-12) do x
                 global count += 1
                 x[1]+x[2]<1
              end
(0.5000000962250204, 4.9994842542085703e-5)

julia> count
214659

The only reason @gideonsimpson is getting "abysmal" performance is that you are using the default error tolerance of 1e-8.

Moral: When you have a non-smooth integrand, you have to set a higher error tolerance. It's really hard to get extremely good accuracy when integrating non-smooth functions. (And the error estimate is less accurate too.)

[antoine-levitt]

I just ran into a similar problem again. I get strange noise with Cubature sometimes that I don't get with Cuba. It doesn't seem to disappear when decreasing the tolerance or increasing the maximum number of iterations. MWE here (a bit ugly, but at least self-contained): https://gist.github.com/antoine-levitt/bfbb83be1bd4ecdf08b4b95b05aa1a4c

[antoine-levitt]

I got a problem where it was more pronounced but can't reproduce it anymore, and this should be fine as an MWE. To be clear the problem is around .5, where cubature returns a slightly wrong answer, even at very small tolerances.

[stevengj]

@antoine-levitt, are you running into something like #25? Can you try the HCubature package instead?

[antoine-levitt]

I don't know, possibly ? Hcubature does fix the slight bump, I'll use that next time I have to integrate something (right now I've switched to Cuba and it works well enough for what I want to do). Integrating discontinuous functions remains a huge pain. There are in the DFT literature methods specialized in integrating f(x) on the set g(x) >= 0 by interpolating g, but they are a pain to implement...

[stevengj]

Definitely discontinuous functions are hard to do well in general.

Yes, in principle if you have an good description of the boundary, e.g. as a level set of a smooth function, you can do much better, but I agree that they seem like a pain to implement. Would make a great Julia project though — this is the right language for quadrature algorithms.

[stevengj]

See e.g. High-Order Quadrature Methods for Implicitly Defined Surfaces and Volumes in Hyperrectangles

[antoine-levitt]

Thanks! I looked into papers from the math side on that topic a while ago and couldn't find any, but that one and the references within are very relevant!