Speeding up integral3 command is taking too long to execute

I am using MATLAB to execute a triple integral using integral3 and it is running very slow. I was wondering if there ways to speed it. Below is my program. Note that some of the variables in my code are nonnegative constants: A1, a1, B1, b2, lambda1, lambda2, delta1, delta2, theta..

    clear all
    syms gamma1
    syms gamma2
    syms z
    syms v 

    Nt=16; sigmanoise=10^(-7.9); c3=0.129; c1=(1-c3)/2;a2=0;b2=0; 
    a1=0.0030; b1= 0.0030; A1=  1.5625e-04,A2=0;  B1= 7.8125e-05;B2=0;
    theta= 3.1623;lambda1=  4.9736e-05;lambda2=0;p1=1;p2=0; alpha1=2; alpha2=4;delta1=2/alpha1; delta2=2/alpha2;beta1=0.025; beta2=0.025;

    a= gamma1^-1+gamma2^-1+2*gamma1^(-0.5)*gamma2^(-0.5);
    laplacesgi=(exp(+2*pi*j.*z*a)-1)./(2*pi*j*z);
    laplacesgi=matlabFunction(laplacesgi);

    laplacenoi=exp(-2*pi*j.*z*theta*sigmanoise/Nt);
    laplacenoi=matlabFunction(laplacenoi);

    interfere= @(gamma1,gamma2,v,z)( (1 -2*c1-c3./(1+2*pi*j*z*theta*v.^(-1))).*(A1.*(v).^(delta1-1).*exp(-a1.*(v).^ (delta1./2))+B1.*(v).^(delta2-1) .*(1-exp(-b1.*(v).^ (delta2./2))))); 
    gscalar =@(gamma1,gamma2,z)integral(@(v)(interfere(gamma1,gamma2,v,z)),gamma2,inf);
    g = @(gamma1,gamma2,z)arrayfun(gscalar,gamma1,gamma2,z);

    lp= A1*(gamma1)^(delta1-1)*exp(-a1*(gamma1)^ (delta1/2))+B1*(gamma1)^(delta2-1)*(1-exp(-b1*(gamma1)^ (delta2/2)))+A2*gamma1^(delta1-1)*exp(-a2*gamma1^(delta1/2))+ B2*gamma1^(delta2-1)*(1-exp(-b2*gamma1^ (delta2/2)));%;
    dk1=((2*pi*lambda1))/(beta1^2)*(1-exp(-a1*(gamma2)^(delta1/2))*(1+(gamma2)^(delta1/2)*a1))+ pi*lambda1*gamma2^(delta2)*p1^delta2-((2*pi*lambda1)/(beta1^2))*(1-exp(-b1*(gamma2)^(delta2/2))*(1+(gamma2)^(delta2/2)*b1));
    dk2=((2*pi*lambda2))/(beta2^2)*(1-exp(-a2*(gamma2)^(delta1/2))*(1+(gamma2)^(delta1/2)*a2))+ pi*lambda2*gamma2^(delta2)*p2^delta2-((2*pi*lambda2)/(beta2^2))*(1-exp(-b2*(gamma2)^(delta2/2))*(1+(gamma2)^(delta2/2)*b2));
    dk=dk1+dk2;
    lcp= A1*(gamma2)^(delta1-1)*exp(-a1*(gamma2)^ (delta1/2))+B1*(gamma2)^(delta2-1)*(1-exp(-b1*(gamma2)^ (delta2/2)))+A2*gamma2^(delta1-1)*exp(-a2*gamma2^ (delta1/2))+ B2*gamma2^(delta2-1)*(1-exp(-b2*gamma2^(delta2/2)));%;

    pdflast=lp*lcp*exp(-dk);
    pdflast=matlabFunction(pdflast);
    pdflast= @(gamma1,gamma2)arrayfun(pdflast,gamma1,gamma2);
    gamma2min=@(gamma1)gamma1;
    warning('off','MATLAB:integral:MinStepSize');
    T = integral3(@(gamma1,gamma2,z)(laplacenoi(z).*laplacesgi(gamma1,gamma2,z).*pdflast(gamma1,gamma2).*exp(-g(gamma1,gamma2,z))),0,inf,@(gamma2)gamma2,inf,0.05,1000,'abstol',1e-3)

I appreciate any ideas or suggestions.

Answer 1

[Mike Hosea]

These two articles might help.

http://blogs.mathworks.com/loren/category/integration/

Try different variable orderings. Make sure your tolerances are sane (do this first). Also, consider whether it is possible to speed up the integrand function.

[DM]

Thanks, is 1e-30 a sane number for abstol?

[Mike Hosea]

First of all, there aren't many instances in which it makes any sense to supply just AbsTol. The actual tolerance applied is the LEAST restrictive of the absolute and relative tolerances. Normally it is the RelTol that controls the integration, but a reasonable AbsTol is needed for the cases where the true answer appears to be close to zero. It is a question of whether the size of the numbers are small because of the scaling (e.g. doing a terrestrial distance calculation in light years instead of nanometers, for example) or because of cancellation. In the former case, the leading zeros in a result like, say, 0.000000000012345 are probably "free" in a computational sense. When it is cancellation, however, they are not free--at least some of those leading zeros have to be worked for. In double precision it is only possible to compute at most around 15 of them, more likely fewer, depending on the specific problem.

An AbsTol of 1e-30 is basically asserting, in the context of double precision arithmetic, that there are more than 14 leading zeros to the right of the decimal point at all times during the calculation. If not, then it is not a "sane" absolute error tolerance.

[Mike Hosea]

Now that I see more of your problem, I would probably want to spend some time seeing if I could get rid of the arrayfun usage. I'll play with it a little bit later...bit too busy just now.

[DM]

Mike, thanks I just saw your comment (the second one), if you copy the code as is it should work with no syntax errors. I am looking forward for your help !

[Mike Hosea]

I wanted to see whether I could use MATLAB Coder to speed this up. MATLAB Coder doesn't support anonymous functions, INTEGRAL, or INTEGRAL3, but it does have QUADGK, and there's a trick using PERSISTENT that gets around the lack of anonymous functions. I could very well have made a mistake in the translation, and maybe you don't even have a MATLAB Coder license, but perhaps the approach, eschewing SYMS at least and coding the integrand directly, trying to re-use as much work as possible, might be of some help. If you don't have MATLAB Coder, you can still integrate @integrand directly without compiling it. It might be faster (or slower).

Anyway, here it is:

>> tic; T = integral3(@integrand_mex,0,inf,@(gamma1)gamma1,inf,0.05,1000,'AbStol',1e-5,'RelTol',1e-3), toc
T =
      5.329614020381315e+00 - 2.196183908954074e+02i
Elapsed time is 6463.696856 seconds.

The function integrand.m was compiled like so:

cfg = coder.config('mex');
cfg.IntegrityChecks = false;
codegen integrand -args {t,t,t}

function w = integrand(gamma1,gamma2,z)
persistent Nt r1 r2 r3 r4 a1 a2 b1 b2 delta2d2 delta1d2 twopij p1 p2 c1 c3 A1 B1 A2 B2 delta1 delta2 p1d2 p2d2 lambda1 lambda2
if isempty(Nt)
    Nt = 16;
    sigmanoise = 10^(-7.9);
    c3 = 0.129;
    c1 = (1 - c3)/2;
    a2 = 0;
    b2 = 0;
    a1 = 0.0030;
    b1 = 0.0030;
    A1 = 1.5625e-04;
    A2 = 0;
    B1 = 7.8125e-05;
    B2 = 0;
    theta = 3.1623;
    lambda1 = 4.9736e-05;
    lambda2 = 0;
    p1 = 1;
    p2 = 0;
    alpha1 = 2;
    alpha2 = 4;
    delta1d2 = 1/alpha1;
    delta2d2 = 1/alpha2;
    delta1 = 2*delta1d2;
    delta2 = 2*delta2d2;
    beta1 = 0.025;
    beta2 = 0.025;
    beta1sq = beta1*beta1;
    beta2sq = beta2*beta2;
    r1 = 2*pi*lambda1/beta1sq;
    r2 = 2*pi*lambda2/beta2sq;
    twopij = 2*pi*1j;
    r3 = -2*pi*1j*theta*sigmanoise/Nt;
    r4 = twopij*theta;
    p1d2 = p1^delta2;
    p2d2 = p2^delta2;
end
a = gamma1.^-1 + gamma2.^-1 + 2*gamma1.^(-0.5).*gamma2.^(-0.5);
twopijz = twopij*z;
laplacesgi = exp(twopijz.*a - 1)./twopijz;
laplacenoi = exp(r3*z);
g = coder.nullcopy(complex(gamma1));
for k = 1:numel(gamma1)
    interfere('set',z(k),a1,b1,c1,c3,r4,A1,B1,delta1,delta2,delta1d2,delta2d2);
    g(k) = quadgk(@interfere,gamma2(k),inf);
end
g1td1 = gamma1.^delta1d2;
g1td2 = gamma1.^delta2d2;
g2td1 = gamma2.^delta1d2;
g2td2 = gamma2.^delta2d2;
a1g2td1 = a1*g2td1;
a2g2td1 = a2*g2td1;
b1g2td2 = b1*g2td2;
b2g2td2 = b2*g2td2;
ea1g2td1 = exp(-a1g2td1);
ea2g2td1 = exp(-a2g2td1);
eb1g2td2 = exp(-b1g2td2);
eb2g2td2 = exp(-b2g2td2);
g1d1 = gamma1.^(delta1 - 1);
g1d2 = gamma1.^(delta2 - 1);
g2d1 = gamma2.^(delta1 - 1);
g2d2 = gamma2.^(delta2 - 1);
lp = ...
    A1*g1d1.*exp(-a1*g1td1) + ...
    A2*g1d1.*exp(-a2*g1td1) + ...
    B1*g1d2.*(1 - exp(-b1*g1td2)) + ...
    B2*g1d2.*(1 - exp(-b2*g1td2));
dk1 = ...
    r1*(1 - ea1g2td1.*(1 + a1g2td1)) + ...
    pi*lambda1*gamma2.^delta2*p1d2 - ...
    r1*(1 - eb1g2td2.*(1 + b1g2td2));
dk2 = ...
    r2*(1 - ea2g2td1.*(1 + a2g2td1)) + ...
    pi*lambda2*gamma2.^delta2*p2d2 - ...
    r2*(1 - eb2g2td2.*(1 + b2g2td2));
dk = dk1 + dk2;
lcp = ...
    A1*g2d1.*ea1g2td1 + ...
    A2*g2d1.*ea2g2td1 + ...
    B1*g2d2.*(1 - eb1g2td2) + ...
    B2*g2d2.*(1 - eb2g2td2);
pdflast = lp.*lcp.*exp(-dk);
w = ...
    laplacenoi.* ...
    laplacesgi.* ...
    pdflast.* ...
    exp(-g);

function y = interfere(v,zi,a1i,b1i,c1i,c3i,r4i,A1i,B1i,delta1i,delta2i,delta1d2i,delta2d2i)
persistent z a1 b1 c1 c3 r4 A1 B1 delta1 delta2 delta1d2 delta2d2
if ischar(v)
    z = zi;
    a1 = a1i;
    b1 = b1i;
    c1 = c1i;
    c3 = c3i;
    r4 = r4i;
    A1 = A1i;
    B1 = B1i;
    delta1 = delta1i;
    delta2 = delta2i;
    delta1d2 = delta1d2i;
    delta2d2 = delta2d2i;
    return
elseif isempty(z)
    z = 0;
    a1 = 0;
    b1 = 0;
    c1 = 0;
    c3 = 0;
    r4 = complex(0);
    A1 = 0;
    B1 = 0;
    delta1 = 0;
    delta2 = 0;
    delta1d2 = 0;
    delta2d2 = 0;
end
y = (1 - 2*c1 - c3./(1 + r4*z./v)).* ...
    (A1*v.^(delta1 - 1).*exp(-a1*v.^delta1d2) + B1*v.^(delta2 - 1).*(1 - exp(-b1*v.^delta2d2)));

[DM]

Mike, thank you very much, can I use this program even though I don't have MATLAB coder?

[Mike Hosea]

Yeah, it should work, but trying things out it looked like compiling it with MATLAB Coder speeded it up about 10x. I'm running a slightly modified version now that uses PARFOR for the loop computing quadgk.

[Mike Hosea]

BTW, just FYI (since you don't have MATLAB Coder), the loop doing quadgk can be a "parfor" in MATLAB Coder. You don't need the Parallel Computing Toolbox for that since it's only multithreading (i.e., not distributed). I do not think it would be a good idea to make it a parfor in MATLAB, but I don't have that much experience with the Parallel Computing Toolbox, so I don't have a good feeling for how much overhead there will be.

But with Coder's multithreading, I get (I've got several spare cores, so YMMV)

>> tic; T = integral3(@integrand_mex,0,inf,@(gamma1)gamma1,inf,0.05,1000,'AbStol',1e-5,'RelTol',1e-3), toc
T =
      5.329614020381315e+00 - 2.196183908954074e+02i
Elapsed time is 1434.025262 seconds.