How to reduce computation time

Hi,

The function "funct" takes up almost 0.06 seconds and it is the first responsible for a slow execution of an integral (more than 3 hours) when it is used in it with regards to y.

Any help please to reduce the computation time of the following function

function [out] = funct(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x)
nx=length(x);
out=zeros(1,nx);
A=(mu./pi).*(wmax+lmax);
B=(mu.*wmax.*lmax)./4;
C1=(hb-(2.*hmax))./(2*hmax);
C2=(hu+hb-(2.*hmax))./(2*hmax);
for i=1:nx
   S1=@(u,theta) u.*(1-exp(C1.*((A.*sqrt(u.^2+x(i).^2-(2.*x(i).*u.*cos(theta))))+B)));
   S2=@(theta) (1-exp(C1.*((A.*sqrt(y.^2+x(i).^2-(2.*x(i).*y.*cos(theta))))+B)));
   out(i)=lambdab.*y.*integral(S2,0,2*pi,'ArrayValued',true).*exp(-lambdab.*quad2d(S1,0,y,0,2*pi))./(1-exp(-lambdab.*quad2d(S1,0,1000,0,2*pi)));
end
end

[per isakson]

Help us to help by providing a typical set of input data for funct so we can run the profiler.

[Yassine Hmamouche]

Thank you so much for the feedback.

All input data are positive real numbers such that C1 and C2 are negative numbers. Hereafter some typical values:

mu=1; lmax=4;wmax=6;hmax=30;lambdab=1; hu=40; hb=10;y=4;x=2;

Best,

[Walter Roberson]

On my system, I time as closer to 0.012 .

[Yassine Hmamouche]

Yes this is the self computing time of the function, but when used as the integrand of another integral it gives a huge total computing time (i.e., more than 4 hours).

Answer 1

[Walter Roberson]

By going back and replacing that 1000 limit with infinity, and adding the assumption that all values are non-negative, and that hb < 2*hmax, Maple is able to combine terms to remove one integral. The symbolic form is

-exp(lambdab * int( int( u * (exp(1/2 * (1/4 * lmax * wmax * pi + (lmax+wmax) * ...
    (u^2 - 2 * cos(theta) * u * x + x^2)^(1/2)) * mu * (hb-2*hmax) / pi / hmax)-1), u, 0, y), ...
    theta, 0, 2*pi)) * lambdab * y * int(-1 + exp(1/8 * ((4 * lmax + 4 * wmax) * ...
    (x^2 - 2 * cos(theta) * x * y + y^2)^(1/2) + lmax * wmax * pi) * (hb-2*hmax) * mu / pi / hmax), ...
    theta, 0, 2*pi)

This should be more efficient

[Yassine Hmamouche]

This solution seems very interesting. I will even try to simplify the first expression (i.e., quad2d(S1,0,y,0,2*pi))

Let me check it and come back to you.

best,

[Yassine Hmamouche]

Please share the resulted simplification. I did not get anything by using MATHEMATICA.

[Walter Roberson]

What I posted here is the simplified version of funct. funct included one integral() and two quad2d for a total of 5 integrations. The simplified version (posted above) has fewer integrations.

[Yassine Hmamouche]

I think you did not consider in MAPLE the expression of the second quad2d (expression of the denominator, or f3 in the script of isakson).

Either way, MAPLE has only splitted quad2D into int(int(.)). So the proposed script will contain: 5 int(.) instead of 1 int(.) + 2 quad2D, knowing that quad2D is largely efficient than int(int(.)).

I summarize here validated actions that might save some mseconds of delay:

• Reduce the interval of integration from infinity to 1000 or even 100 (I prefer rapid results with acceptable errors than very delayed accurate results).
• Compute constants (e.g., A,B,C1) outside the function "funct".
• upgrade at least to R2018 to be able to change quad2D by a possibly more efficient integral2.

Best regards,

[Walter Roberson]

you would be incorrect. I converted the quad2d calls into nested int() calls, giving an original expression that had 5 int() calls in it. Then by adding some reasonable assumptions, and changing the 1000 to infinity, I was able to get maple to reduce the 5 int() calls into 3 int() calls.

[Yassine Hmamouche]

I think the integral of the second quad2D seems removed/simplified in https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#answer_438033

because it hasn't been considered in Maple.

[Walter Roberson]

You can see in the below that I input all 5 levels of integration, and that Maple believes that it can be reduced to 3 integrations.

[Yassine Hmamouche]

Dear Walter and isakson,

Thank you so much for your feedback and effort.

You have demonstrated in fact that we have a great symbolic simplification. The function f3 in

https://fr.mathworks.com/matlabcentral/answers/531213-how-to-reduce-computation-time#answer_437408

equals 1 all the time.

The processing time reduces then as

 tic
for jj = 1 : 10000
    out1 = f_Y0_X0_bis_bis(4,3,-4,-3,1,20,60,10,4,1); % Old function
end
toc
%%
tic
for jj = 1 : 10000
    out2 = f_Y0_X0_bis_bis2(4,3,-4,-3,1,20,60,10,4,1);% Function after simplification
end
toc
disp(out2-out1)
%%%%%%%%%%%%%%%
Elapsed time is 55.496311 seconds.
Elapsed time is 39.585332 seconds.

Let's summarize here the lessons learned from these discussions:

• Use a recent version of Matlab and a new vanilla PC.
• Doing steps 0, 1 and 2 manually or through powerful symbolic tools like Maple and Mathematica before considering vectorial analysis through Matlab.
• Avoid computing unecessary tasks inside loops, e.g., constants, wide integral ranges.

Many thanks and best regards.

[Walter Roberson]

Also, it helps to know any constraints that might apply. Even just knowing that a particular variable is real-valued might be enough to permit a calculation to simplify.

Answer 2

[per isakson]

On R2018b and a new vanilla PC I found a significant decrease in computation time by replacing quad2d() by integral2() (with the same signature, i.e. defaults are used). With your input data this replacement didn't affect the result, out. (Caveat, I do this because I'm curious, not because I know much about the topic.) There are two posts on Double Integration in MATLAB ( and here ) at Loren on the Art of MATLAB. I don't find backing for the difference in speed that I've seen, which makes me a bit nervous. What's the catch? There ought to exist a discussion on how to chose between the two functions, but I fail to find it.

Have you tried to simplify the expressions with the help of the Symbolic Toolbox? See https://se.mathworks.com/matlabcentral/answers/33399-problem-with-double-integral-dblquad-and-quad2d#answer_42126 it's a bit old but shows the idea. Maybe that toolbox sees something that I don't see.

In response to comments

The difference in speed between integral2(), introduced in R2012a, and quad2d(), introduced in R2009a, that I see is irking.

I've compared your function, funct(), to funct_int2(), in which I replaced quad2d() by integral2(). (And split the long statement into three to make the code and the result of profile() easier to read.)

>> cssm % first call to cssm after starting Matlab
Elapsed time is 0.364935 seconds.
Elapsed time is 0.140186 seconds.
     0
>> cssm
Elapsed time is 0.104307 seconds.
Elapsed time is 0.021930 seconds.
     0
>> cssm
Elapsed time is 0.088423 seconds.
Elapsed time is 0.016356 seconds.
     0
...
...    
>> cssm
Elapsed time is 0.091209 seconds.
Elapsed time is 0.017180 seconds.
     0     

where cssm.m is a script containing

%%
mu=1; lmax=4; wmax=6; hmax=30; lambdab=1; hu=40; hb=10; y=4; x=2;
%%
tic
for jj = 1 : 10
    out1 = funct(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x);
end
toc
%%
tic
for jj = 1 : 10
    out2 = funct_int2(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x);
end
toc
disp( out2-out1 )

and

function [out] = funct_int2(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x)
    nx=length(x);
    out=zeros(1,nx);
    A=(mu./pi).*(wmax+lmax);
    B=(mu.*wmax.*lmax)./4;
    C1=(hb-(2.*hmax))./(2*hmax);
    C2=(hu+hb-(2.*hmax))./(2*hmax);         %#ok<NASGU>
    for ii=1:nx
        S1=@(u,theta) u.*(1-exp(C1.*((A.*sqrt(u.^2+x(ii).^2-(2.*x(ii).*u.*cos(theta))))+B)));
        S2=@(theta) (1-exp(C1.*((A.*sqrt(y.^2+x(ii).^2-(2.*x(ii).*y.*cos(theta))))+B)));
        f1 = lambdab.*y.*integral(S2,0,2*pi,'ArrayValued',true);
        f2 = exp(-lambdab.*integral2(S1,0,y,0,2*pi));
        f3 = (1-exp(-lambdab.*integral2(S1,0,1000,0,2*pi)));
        out(ii) = f1.*f2./f3;
    end
end

Comments

• Fine print; don't blame me for mistakes.
• In this comparison both integral2() and quad2d() use the "tiled method"
• The fact that the two return exactly the same result makes me believe that deep inside the same calculation is made.

[Yassine Hmamouche]

Thank you so much for the feedback.

Yes I have tried to simplify the expressions, but without much success. The expressions are already very compact.

I also compared the two options: i) using Integral(integral) or ii) quad2d. I found that the second option is by far the fastest.

I think that one important solution is trying to vectorize the loop. Please any idea to do that !

function [out] = funct(mu,lmax,wmax,hmax,lambdab,hu,hb,y,x)
nx=length(x);
out=zeros(1,nx);
A=(mu./pi).*(wmax+lmax);
B=(mu.*wmax.*lmax)./4;
C1=(hb-(2.*hmax))./(2*hmax);
C2=(hu+hb-(2.*hmax))./(2*hmax);
for i=1:nx
   S1=@(u,theta) u.*(1-exp(C1.*((A.*sqrt(u.^2+x(i).^2-(2.*x(i).*u.*cos(theta))))+B)));
   S2=@(theta) (1-exp(C1.*((A.*sqrt(y.^2+x(i).^2-(2.*x(i).*y.*cos(theta))))+B)));
   out(i)=lambdab.*y.*integral(S2,0,2*pi,'ArrayValued',true).*exp(-lambdab.*quad2d(S1,0,y,0,2*pi))./(1-exp(-lambdab.*quad2d(S1,0,1000,0,2*pi)));
end
end

Best regards.

[Walter Roberson]

Vectorized loop would move from a 2d integral to 3d, but there is no quad3d. Also the integral() would have to be integral2() but integral2() does not have ArrayValued

[Yassine Hmamouche]

Thank you so much for the feedback.

Yes, I tried to make vectorization but it does not work with quad2.

Alternatively, I tried to gain some mseconds in the processing delay by reducing the interval of integration. For instance, 0 to infinity to be replaced by 0 to 1000 or even by 0 to 100. Also, I tried to make constants A, B, C1, et C2 outside the function "funct" and then as inputs of it rather than computed inside it.

Best regards

[per isakson]

@Yassine, which release of Matlab are you using?

I've added to my answer.

[Yassine Hmamouche]

@per isakson, I am using Matlab R2013a.

[per isakson]

@Yassine Hmamouche, Then you have integral2() and should be able to rum "my" funct_int2. See integral2, Numerically evaluate double integral. I'm curious to hear whether you can reproduce my result.

[Yassine Hmamouche]

I tried to use integral2() but it doesn't work.

R2013a has typically the following integral functions:

• For one fold integral: integral, quad, ...,etc.
• For two fold integral: quad2D.
• For three fold integral: integral3.