MATLAB: I want a numerical or fast multidimensional integral

numarical integration

I have the following code:

 if true
clc
clear all
close all
tic
syms x1 x2 x3 x4 x5
warning off;
X=[x1-12.568; x2+9.568; x3-11.486; x4-4; x5-0.8];%this is [x-E[x]] matrix.
Cx=[39.4573    0.2350    1.5190   25.6028    1.6754; 0.2350    0.0089    0.0285    0.5958    0.0344;...
    1.5190    0.0285    0.2337    1.8325    0.1193; 25.6028    0.5958    1.8325  667.1131    5.5158;...
    1.6754    0.0344    0.1193    5.5158    1.3541];%this the covariance matrix.
DeterminantCx=det(Cx)
CxInverse=inv(Cx);%this is the inverse matrix of the covariance matrix (Cx).
DeterminantCxInverse=det(CxInverse)
Xtranspose=transpose(X);
Exponent=Xtranspose*CxInverse*X;
fJointDensity=(1/(sqrt(DeterminantCx)*(2*pi)^(5/2)))*exp(-(Exponent)/2);
I1=int(fJointDensity,x1,-inf,inf);
I2=int(I1,x2,-inf,inf)
I3=int(I2,x3,-inf,inf)
I4=int(I3,x4,-inf,inf)
I5=int(I4,x5,-inf,inf)
t=['Time Elapsed=',num2str(toc), ' sec'];
disp(t)
end

The expected output for I5 is 1, but the problem is that performing the integrations using "int" takes a lot of time. This program is an example of my main problem, but if I can run this program very fast then I can solve my problem.

So I am searching for a numerical integration or a fast integration method.

Thanks a lot.

Hadi.

Best Answer

(I've got a bit more time now, and since Mike has not chimed in, I'll offer an answer.)

The point is, there are things you MIGHT do, but how can we possibly help you? What do you know or not know about numerical integration? What characteristics does your function have? Are the limits really infinite, or are they finite? The problem is, one could write a book about integration. Oh, thats right, there are already lots of books on the subject. But without knowing something about what you have to solve, the best answer is to buy a faster computer.

There are no n-fold integrations in MATLAB, but I recall that Mike Hosea has posted a higher degree integration tool on the FEX, here:

http://www.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

It handles 4,5,6 fold integrations, I think even allowing infinite limits.

That is the good news. The bad news is to expect it to be seriously slow. Iterated integration tools like this tend to take a lot of function evaluations. The reason is simple. To understand, lets make a simple test...

function [f,nx] = testfun(x)
persistent nevals
if isempty(nevals)
  nevals = 0;
end
nevals = nevals + numel(x);
f = 1/sqrt(2*pi)*exp(-x.^2/2);
if nargout > 1
  nx = nevals;
end

That little function counts the number of function evaluations integral makes for a simple Gaussian.

integral(@testfun,-inf,inf)
ans =
            1
[f,nx] = testfun(0)
f =
      0.39894
nx =
   211

So to compute the integral of ONE gaussian in one dimension, it took 211 function evaluations.

But you have a 5-fold integral. Each of these integrations is of similar complexity, so you might expect to need

211^5
ans =
   4.1823e+11

so about 4e11 function evaluations. That is 400 billion evals. If your function has 600 lines of code, how long will it take to execute? If it is a bit slow, and you get only 1 eval per second, then it would take

211^5/30e6
ans =
        13491

YEARS to finish. Even if one call to your 600 line function takes 0.001 seconds to execute, that is still over 13 YEARS to terminate.

Even if some of those integrations are simpler, and take only half as many calls, go get a cup of coffee and a good book, as you will be waiting for a bit. "War and Peace" might be a good choice of books for this one, or better yet, the complete works of Shakespeare, and expect to re-read whatever book you have quite a few times over. The good news is when it gets done, you will be a true Shakespeare scholar.

My point about not telling us anything is still valid. Are there symmetries you can take advantage of? Or, if your function has characteristics like that of a Gaussian as you show, do you know anything about Gauss-Hermite quadrature? This could easily be a boon for you, or a bust. I cannot know, since you have told us nothing.

Related Solutions

MATLAB: Multidimensional numerical integration! is there any function to do that?!!

You do not say if the domain of integration is a simple one, perhaps a hypercube in R^6. It is not really relevant, but if the domain is more complex, that just makes it worse. And nothing that Meysam has posed will help you.

Numerical integration in a high number of dimensions is a nasty thing. And, sadly, 6 dimensions is what I would call moderately high. Why? Lets see what happens when you perform a single dimensional numerical integration. Commonly, tools like integral, or one of the quad tools, will require between 100 and 1000 function evaluations to compute an integral estimate in one dimension. That number of evals will depend on many things, but mainly I would think of the complexity of the integrand and the integration tolerance.

If your function is at all complicated, with locally high derivatives or even singularities, spikes, sharp bends, etc., then the integration will need to resolve where those problems lie, and use many points to pin down the best estimate of the integral. And naturally, the integration tolerance is important, since if you insist on a better estimate, then it will take more time and effort to get it for you. remember that any numerical integration, if it is to be accurate, must put many points in any region of importance. But adaptive numerical integrations do not know what your function looks like. All they do is evaluate the function, and then try to infer where the complexities lie from those function evals. Any such tool treats your objective asa simple black box. You pass numbers in, get numbers out. But inside the kernel, it is a black box that you cannot see into.

You might think that well, it makes sense then to just use a tool like Simpson's rule or trapezoidal rule. Just evaluate the function at a string of points. But then you must recognize that not much is happening at alost all of those points, so most of them are totally wasted function evals. So you need to accept the idea that 100-1000 points for a 1-d problem is about as good as you can do in a classical sense, IF the integration knows nothing about the integrand.

What happens in multiple dimensions? Now it gets exponentially nasty. An adaptive tool, in two dimensions will typically require at least on the order of 100^2 - 1000^2 function evaluations, so often at least 10000 function evaluations. In n dimensions (here n==6) we are now talking about potentially 100^6 = 10^12 function evals. And that may be on the low side. If your function is complicated, as I said, you may need to worry about something on the order of 1000^6 = 10^18 evals. This is the curse of dimensionality for nuemrical integration. 10^12 function evaluations is a TRILLION function evals. That will take some serious time to do.

So, what can you do? You do have several options. There is one numerical integration that does not suffer as badly from such a curse - Monte Carlo integration. In fact, the beauty of Monte Carlo is it can survive a 6 dimensional problem in a finite amount of time. And Monte Carlo is easy to perform too. The problem is that Monte Carlo integration will not give you a highly accurate estimate. In fact, you can estimate how the integration error for Monte Carlo performs with the number of evals. It tends to decrease as 1/sqrt(n), where n is the number of sample points, and that behavioral decrease will not depend on the number of dimensions!

A very important idea is often if you can do the integration analytically, in at least one of the necessary dimensions. If you can decrease the dimensions, of course then things get exponentially better.

Finally, you can consider tools like a tensor product Gaussian integration. But that typcially requires you know a very specific behavior of the integrand, and we have been told nothing in this respect. Sorry, but if you really want help, then I would need to know more about your problem.

So, without such specific information, I would strongly suggest you consider a Monte Carlo integration, which is actually trivially easy to compute. As a simple example, consider the function

Z = x1 + x2 + x3 + x^4 + x5 + x6

over the region from 0 to 1 in each of the 6 dimensions present. I've chosen it as an example, because it is a simple one where we know the answer to compare it to. Since we know

syms x
int(x,[0,1])

So 1/2 in one dimension. In 6 dimensions, the integral will be 3. A Monte carlo just has us sum the function values, divide by the area of the integration domain (here that ares is diff(limits)^6), and then divide by the number of sampels.

fun = @(X) sum(X,2);
intfun = @(N,limits) sum(fun(rand(N,6)*diff(limits) + limits(1)))/diff(limits)^6/N;
format long g
intfun(100,[0,1])
ans =
          3.02889547585076
intfun(10000,[0,1])
ans =
          3.00735201210059
intfun(10000000,[0,1])
ans =
          2.99998533404415

So as we increased the sample size from 100 to 10 million, the error of integration goes down slowly. But you can see how trivially easy it is to write a Monte Carlo numerical integration.

I could have done far better if I knew the integrand itself of course.

MATLAB: Is it possible to calculate a 10- Dimensional integral

Is it possible to compute the 10-d integral shown here, or something like that, but not shown, and far more complicated? Almost always when someone asks a question as you have posed, they really have a far different and far more complicated function.

So rather trivially, the 10-fold integral can be done symbolically, yielding 155/6. (Do I really need to write it out?) TRY IT!

syms a b c d e f g h i j
K = expand((a+b+c+d+e+f+g+h+i+j).^2);
K = int(K,j,[0 1]);
K = int(K,i,[0 1]);
K = int(K,h,[0 1]);
K = int(K,g,[0 1]);
K = int(K,f,[0 1]);
K = int(K,e,[0 1]);
K = int(K,d,[0 1]);
K = int(K,c,[0 1]);
K = int(K,b,[0 1]);
K = int(K,a,[0 1])
K =
155/6
double(K)
ans =
        25.8333333333333

Can you solve the problem numerically as a nested integration on some far more general nonlinear function? Well, in theory yes, but not really in practice.

Typically numerical integration tools are adaptive routines. Call them in a nested form though, and they take time to do. So a 1-d integration, even on a quite simple function still takes on the order of 100 function evaluations (but it may require many more evals too.) Expect a 2-d numerical integration to require something like 10000 function evals, etc. So I would predict that a 10-fold nested adaptive integration will use something on the rough order of 100^10=1e20 function evaluations.

That is a LOT of time to do. If some of those dimensions end up being a bit more nasty to integrate on, or if your kernel is a bit more complex in reality, then it might take more time. Even if your computer can evaluate that kernel a billion times per second, we are still talking thousands of years to finish the task. Of course, in a few thousand years from now, if you are still around, you might have a REALLY fast computer.

So no, 10-fold nested adaptive integrations are a waste of time to even try.

Typically when you start to look for very high order nested integrations, Monte Carlo integration is a far better bet. So in way less than a second of computing time, I get:

double(K)
ans =
          25.8333333333333
X= rand(10000000,10);
sum(sum(X,2).^2)/10000000
ans =
           25.826309982603

Not bad for maybe a second of CPU time.

A nice thing about the Monte Carlo is you can even compute uncertainty estimates around the result. Thats basic statistics, and I have the flu so my mind is a bit wasted. You can learn to do that.

The reason why Monte Carlo is so good for high dimensional problems is while the nested integration requires exponential time in the number of levels of nesting, the Monte Carlo has an error estimate that is proportional to sqrt(N), where N is the number of function evals taken.

Could I have done this in a different way? Yes, of course. I could, for example, use a 10-fold nested Gauss-Legendre numerical integration. The nested Gauss-Legendre is also exponential in the number of nested levels but the base will be lower than 100. Here, with 5 points in each dimension,

5^10
ans =
     9765625

it would require evaluation of the kernel at roughly 10 million points. But if my function is adequately well represented by a low order polynomial, this should be accurate enough. The trivially simple kernel you have posed would be exactly integrable using a 2^10 set of points as a nested Gauss-Legendre quadrature. I'll get you started with this for the nodes:

nodes = (sqrt(3)/3*[-1 1]+1)/2
nodes =
       0.211324865405187         0.788675134594813
XGL = nodes(dec2bin(0:1023) - '0' + 1);
sum(sum(XGL,2).^2)/1024
ans =
        25.8333333333332

So pretty much full precision of a double. But that requires knowing the kernel was truly a very low order polynomial.

Best Answer

Related Solutions

MATLAB: Multidimensional numerical integration! is there any function to do that?!!

MATLAB: Is it possible to calculate a 10- Dimensional integral

Related Question