Speeding up integral3 command is taking too long to execute

Question

DM 2015-2-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/180517-speeding-up-integral3-command-is-taking-too-long-to-execute

编辑： Mike Hosea 2015-2-28

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mike Hosea 2015-2-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/180517-speeding-up-integral3-command-is-taking-too-long-to-execute#answer_169399

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Mike Hosea 2015-2-26

编辑：Mike Hosea 2015-2-26

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 2015-2-28

在 MATLAB Online 中打开

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)));

Mike Hosea 2015-2-28

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 2015-2-28

编辑：Mike Hosea 2015-2-28

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Speeding up integral3 command is taking too long to execute

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Speeding up integral3 command is taking too long to execute

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论