How can i speed up this part of my code?

I have this part of code that i apply for different indns (which are different vectors each time):

    phi1=phi(indn1);
    p1=p(indn1,indn1);
    rtn1=(phi1'*(p1*phi1));
    Ctn1=((rtn1*ep2)/(st*(1+rtn1)));
    Ctn1=max(Ctn1,0);
    Sn1=gammainc(Ctn1*0.5,0.5) ;
    Sn1=1-Sn1;

    phi2=phi(indn2);
    ....
    ....

I tried cellfun and structfun but it was slower. Is there something more i can do to improve the speed?

[Cedric]

EDITED on 04.19

Have you used the profiler to see what is slow? If not, type

 profile viewer

and call your script from there to get a report.

Just one point, ep2 and st are scalars (?)

[amanita]

yes,ep2 and st are scalars.Basically, phi's are vectors and p's are matrices. all the rest end up into scalars. I ran the profiler and i have done changes (i wanted to use chi2cdf but using gammainc was faster). I believe i have done all necessary changes, it cant go faster. I am just wondering if i could vectorize this part of the code. Cant use parfor loops since i am already using all available workers (this is a function called as a cost function to a GA)

[Cedric]

It seems to me that it would be difficult to improve the efficiency of this part of the code and that the approach would have to be changed at a larger scale to improve the efficiency (if it is even possible). I write this comment to bump up your thread and someone else might see it and have some idea.

[amanita]

thanx Cedric

[amanita]

The problem was solved by creating mex files and calling them through the GA btw.

Answer 1

[Roger Stafford]

In your code the multiplication part of doing

 phi1'*p1*phi1

is being repeated over and over again as you select different index combinations. Perhaps you could speed up things, even in your mex method, by doing things in the following way. (I am assuming phi is a column vector: one column, with the same number of rows as there are rows and columns in p.)

First, compute, one time only, the matrix

 A = (phi*phi').*p;

Then for each new indn vector as it comes along do only this

 rtn = sum(sum(A(indn,indn)));

or

 rtn = sum(reshape(A(indn,indn),[],1));

whichever is faster. This step requires only addition operations rather than the multiplications plus additions required for matrix multiplication, and gives the same result.

Also in gammainc use the 'upper' option instead of doing Sn1=1-Sn1.

[amanita]

Thanx Roger! This really helped.

Sum sum was faster than reshape, the 'upper' though i think makes it a little bit slower.