CODEWITHSUNDEEP
performancematlabsymbolic-mathintegral
Paul Cai
Paul Cai

Reputation: 33

Speed up Matlab integral calculation using either symbolic or numeric methods?

I'm a Matlab beginner as of two months ago. I used it for my summer project to process MRI images. Recently, I wrote code for the integration shown below. However, both methods are extremely slow. It took a day to run them. How can I improve them to shorten the running time?

  1. Symbolic method:

    syms t;
T=t*ones(79,95,78);
RF1=(1/2)*double(int(((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-1/2)),t,0,inf));
RD1=(3/2)*double(int(((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-3/2)),t,0,inf));

  2. Numeric method:

    fun1=@(T) ((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-1/2));
RF1=(1/2)*integral(fun1,0,inf,'ArrayValued',true);
fun2=@(T) ((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-3/2));
RD1=(3/2)*integral(fun2,0,inf,'ArrayValued',true);

Where x1, y1, z are 79-by-95-by-78 real matrices.

Upvotes: 3

Views: 3978

Answers (1)

horchler
horchler

Reputation: 18484

You are trying to calculate 555,750 integrals, so it will take some amount of time. Here are a few suggestions that might help.

1. First, you can express your equations slightly more efficiently (this will potentially effect the numeric precision of your result depending on the size of your values and the stability of the problem):

fun1 = @(T)1./sqrt((T+x1).*(T+y1).*(T+z));
fun2 = @(T)1./sqrt((T+x1).*(T+y1).*(T+z).*(T+z).*(T+z));

An identical transformation can be performed for the symbolic case. In a test just using random values for smaller (4,524 element arrays) x1, y1, and z arrays, both of your sets of integrals were calculated over four times faster on my machine.

2. You may also be able to get a small increase in speed by turning your anonymous functions into regular M-File functions. In separate files:

function RF1 = fun1(T,x1,y1,z)
RF1 = 1./sqrt((T+x1).*(T+y1).*(T+z));

and

function RD1 = fun2(T,x1,y1,z)
RD1 = 1./sqrt((T+x1).*(T+y1).*(T+z).*(T+z).*(T+z));

Then call them passing x1, y1, and z as parameters:

RF1 = 0.5*integral(@(T)fun1(T,x1,x2,z),0,Inf,'ArrayValued',true);
RD1 = 1.5*integral(@(T)fun2(T,x1,x2,z),0,Inf,'ArrayValued',true);

I see a small, but non-negligible, difference on my machine: about 10% faster for the data I tried. This is likely due to the M-File being JIT compiled separately and because regular functions may be treated slightly differently than anonymous ones.

3. Of course you should also play with the absolute and relative tolerances for integral. You could also see if the vectorized quadv is faster and gives accurate results (note, however, that this function will be removed in a future version of Matlab). Depending on the nature of your data and the functions, it may even be possible to discretize the problem (use a large value for the upper bound instead of Inf) and use a method like trapz to perform the integration, though I don't know if this would be faster either.

For the symbolic case, you might try defining t as real (i.e., syms t real;) and possibly specifying the 'IgnoreAnalyticConstraints' option as true to see if either of these improve the speed. Symbolic math can also be inefficient with memory, so you might also try breaking up the problem and integrating in chunks (say row by row). There are other options besides int that may be faster, but without your actual data values it's hard to say (and pure numeric methods will be much faster if you don't need the precision).

4. Lastly, if you don't need much precision, and care only about speed, you could implement the "fast inverse square root trick" in C mex.

Upvotes: 6

Related Questions