Recommendations for evaluating a 6-D Integral

I've have to integrate a 6-D function that doesn't seem to have an analytical solution.

Can anyone recommend a decent numerical integration/sampling technique or script that might be able to handle this. I've implemented a script that runs a 3-D quadrature on a function that calls another 3-D quadrature to obtain its value but it's using a for loop for the first quadrature and so it's very slow. It takes about 8 hours to evaluate the integral - the quadratures meshgrid is 70x70x70 and evaluates through each point individually.

I have access to a machine with quite a bit of memory (24GB) that should be able to handle fairly large matrices.

Any help at all would be greatly appreciated.

Answer 1

[the cyclist]

I believe that the most efficient algorithms for the evaluation of high-dimensional integrals use Monte Carlo techniques. I recall that there is a good discussion of these techniques (and why they are superior in higher dimensions) in the book Numerical Recipes. I don't have my copy handy, so I can't be more specific.

Answer 2

[Teja Muppirala]

Here's a simple idea. If your integrand is more or less smooth. Evaluate it on a number of different coarse grids, and then extrapolate to estimate the answer.

If your grid spacing is h, then the error should converge with order h^2.

For example:

Integrate F = exp(x+y+z+u+v+w) from 0 to 1 on all 6 variables.

The analytic solution is (exp(1) - 1)^6. We will estimate this using 4 different grids, and then extrapolate to get the final result.

F = @(x,y,z,u,v,w) exp(x+y+z+u+v+w);
I = []; %Estimate ov each grid
h = []; %Spacing of each grid
for n = 5:8; 
    pts = linspace(0,1,n+1);
    pts = conv(pts,[1 1]/2,'valid');
    [X,Y,Z,U,V,W] = ndgrid(pts);
    V = F(X,Y,Z,U,V,W);
    h(end+1) = 1/n;
    I(end+1) = sum(V(:))/(n^6);
end;

P = polyfit(h,I,numel(h));
EstValue = P(end)  %The constant term is what we want
TrueValue = (exp(1)-1)^6
RelativeError = (EstValue - TrueValue)/TrueValue

[Simon]

That's an interesting approach.

Unfortunately, I seems like my function may not be smooth enough for it (or maybe I'm doing something wrong).

The variable names are r, theta, phi, s, gamma and delta. Their limits are zeros to infinity, pi. 2pi, infinity, pi and 2pi, respectively.

(I've approximated infinity as 70 in this case!).

Is dividing by n^6 still the right thing to do considering that I changed the intervals?

a = [1;0;0]

F = @(r,theta,phi,s,psi,delta) exp(-(a(1)-r.*sin(theta).*cos(phi)+s.*sin(psi).*cos(delta)).^2).*...

exp(-(a(2)-r.*sin(theta).*sin(phi)+s.*sin(psi).*sin(delta)).^2).*...

r.*sin(theta).*s.*sin(psi).*exp(-(a(3)-r.*cos(theta)+s.*cos(psi)).^2);

I = []; %Estimate ov each grid

h = []; %Spacing of each grid

for n = 5:20;

%Set Grid Limits and Generate Grid

rlimit = 70;

rpts = linspace(0,rlimit,n+1);

rpts = conv(rpts,[1 1]/2,'valid');

thetapts = linspace(0,pi,n+1);

thetapts = conv(thetapts,[1 1]/2,'valid');

phipts = linspace(0,2*pi,n+1);

phipts = conv(phipts,[1 1]/2,'valid');

slimit = 70;

spts = linspace(0,slimit,n+1);

spts = conv(spts,[1 1]/2,'valid');

gammapts = linspace(0,pi,n+1);

gammapts = conv(gammapts,[1 1]/2,'valid');

deltapts = linspace(0,2*pi,n+1);

deltapts = conv(deltapts,[1 1]/2,'valid');

[X,Y,Z,U,V,W] = ndgrid(rpts,thetapts,phipts,slimit,gammapts,deltapts);

%Evaluate Function

V = F(X,Y,Z,U,V,W);

h(end+1) = 1/n;

I(end+1) = sum(V(:))/(n^6);

end;

P = polyfit(h,I,numel(h));

EstValue = P(end) %The constant term is what we want

Answer 3

[Teja Muppirala]

Assuming that "gamma" and "psi" are the same thing, are you sure your integrand is correct? Look at when

theta = psi = pi/2

phi = delta = pi

r = s

Then your integrand turns into (after a bit of symbolic math)

s^2*exp(1)

Just try it: F(70,pi/2,pi,70,pi/2,pi)

you get exp(-1)*70^2

F(7000,pi/2,pi,7000,pi/2,pi)

you get exp(-1)*7000^2

So as r and s get big, it blows up to infinity.

[Simon]

Hmmm... I'll have a look at it closely tomorrow and get back to you.

And you are correct, gamma and psi are the same thing.