0

I have a function f(x) = 1/(x + a+ b*I*sign(x)) and I want to calculate the integral of

dx dy dz f(x) f(y) f(z) f(x+y+z) f(x-y - z)

over the entire R^3 (b>0 and a,- b are of order unity). This is just a representative example -- in practice I have n<7 variables and 2n-1 instances of f(), n of them involving the n integration variables and n-1 of them involving some linear combintation of the integration variables. At this stage I'm only interested in a rough estimate with relative error of 1e-3 or so.

I have tried the following libraries :

  • Steven Johnson's cubature code: the hcubature algorithm works but is abysmally slow, taking hundreds of millions of integrand evaluations for even n=2.
  • HintLib: I tried adaptive integration with a Genz-Malik rule, the cubature routines, VEGAS and MISER with the Mersenne twister RNG. For n=3 only the first seems to be somewhat viable option but it again takes hundreds of millions of integrand evaluations for n=3 and relerr = 1e-2, which is not encouraging.

For the region of integration I have tried both approaches: Integrating over [-200, 200]^n (i.e. a region so large that it essentially captures most of the integral) and the substitution x = sinh(t) which seems to be a standard trick.

I do not have much experience with numerical analysis but presumably the difficulty lies in the discontinuities from the sign() term. For n=2 and f(x)f(y)f(x-y) there are discontinuities along x=0, y=0, x=y. These create a very sharp peak around the origin (with a different sign in the various quadrants) and sort of 'ridges' at x=0,y=0,x=y along which the integrand is large in absolute value and changes sign as you cross them. So at least I know which regions are important. I was thinking that maybe I could do Monte Carlo but somehow "tell" the algorithm in advance where to focus. But I'm not quite sure how to do that.

I would be very grateful if you had any advice on how to evaluate the integral with a reasonable amount of computing power or how to make my Monte Carlo "idea" work. I've been stuck on this for a while so any input would be welcome. Thanks in advance.

gaugeinvariance
  • 1
  • 3

1 Answers1

0

One thing you can do is to use a guiding function for your Monte Carlo integration: given an integral (am writing it in 1D for simplicity) of ∫ f(x) dx, write it as ∫ f(x)/g(x) g(x) dx, and use g(x) as a distribution from which you sample x.

Since g(x) is arbitrary, construct it such that (1) it has peaks where you expect them to be in f(x), and (2) such that you can sample x from g(x) (e.g., a gaussian, or 1/(1+x^2)).

Alternatively, you can use a Metropolis-type Markov chain MC. It will find the relevant regions of the integrand (almost) by itself. Here are a couple of trivial examples.

Peter O.
  • 32,158
  • 14
  • 82
  • 96
ev-br
  • 24,968
  • 9
  • 65
  • 78
  • Thank you very much for your link to the Markov chain MC. With respect to your first point, wouldn't that be just like VEGAS, only worse? VEGAS is already sampling according to the value of the entire integrand rather the crude approximation to it obtained by ignoring the f's with the cross-terms? – gaugeinvariance May 03 '13 at 14:57
  • Having never used VEGAS, I've only very vague idea of its strengths and/or weaknesses. A simple guided sampling may or may not be better: I just don't know how it behaves in presence of singularities. Devil's in the details, and here a simplest guided sampling might or might not be better. – ev-br May 03 '13 at 18:23