speed-up the given code

The following code is taking a most of the time in my script and I have to run the script many times Is there a way to reduce the time for its execution.

t=0:200000;
lambda=[1,2,3;4,5,6;7,8,9;10,11,12];
for trans=1:size(lambda,1)
  uh_fun{trans}=expconv(lambda(trans,:));
end

uh=zeros(length(t),1);
for i=1:length(uh_fun)
  temp_uh_fun=uh_fun{i};
  f_gamma=zeros(length(t),1);
  for j=1:length(temp_uh_fun)
      f_gamma=f_gamma+temp_uh_fun{j}(t)';   % this line take a lot of time
  end
  uh=uh+p(i)*f_gamma;  % p is an array of scalars
end

I have added expconv function for reference. I need some suggestions to improve the performance. I have also attached a pdf which contains profiler analysis of anonymous function in expconv.

[Geoff Hayes]

Abhinav - without knowing anything about your uh_fun functions (I couldn't load them into my version of MATLAB) it may be difficult to suggest any improvements. When you say that your code is taking a lot of time, what do you mean? Does it take minutes or hours?

[Pieter Hamming]

I get this error when opening the .mat:

 D:\Research\Thesis_work\Structural_uncertainty\MatLab_codes\20180222\expconv.m>@(t)exp(-lambda(i)*t)*C(i) 
 Warning: Could not find appropriate function on path loading function handle

[Pieter Hamming]

As the .mat isn't opening we can't do much, but I still may have a suggestion:

tall arrays are Matlab's way of handling large arrays without causing memory issues. As f_gamma has 200001 rows, this may be memory related. Consider changing into tall arrays.

[Abhinav]

@Hayes I have added more code so that uh_fun can be accessed. However, this is just an example code. I have to run this code for about 8000 times which 5-6 hours.

[Abhinav]

@Peter I have added a code to the questions so that uh_fun can be accessed. Thanks for the suggestion. I will let you know if it is working.

[OCDER]

Why do you need to create all these function handles?

%for example
lambda = 1:10;
C = 1:10;
for j = 1:10
  fh{j} = @(t)exp(-lambda(i)*t)*C(i)
end
%will create something like...
@(t)exp(-1*t) * 1
@(t)exp(-2*t) * 2
@(t)exp(-3*t) * 3

[Abhinav]

@OCDER it is necessary to keep them separate because of density obtained by different rows in lambda is weighted by p(i).

[Are Mjaavatten]

What is the problem? Since p is not given I add a line:

p = rand(size(lambda,1),1);

in your code. I then save it as abhinav.m and time it:

>> tic;abhinav;toc
Elapsed time is 0.023084 seconds.

[Abhinav]

@Are My script calls this line 2016028 times, which makes it very costly. Profiler shows that the line

 f_gamma=f_gamma+temp_uh_fun{j}(t)

is the bottleneck. That's why I want to make it faster.

[Geoff Hayes]

Abhinav - is the temp_uh_fun exactly like what you have shown us in the attached code (expconv.m) or is it something else? It is probably expensive to call exp so many times especially when your t is an array of 200001 elements. Can you reduce the number of calls to exp? For example, suppose you do

 expT = exp(t);

Then to get exp(-1*t) then this should be the same as

 expTA = 1 ./ expT;

and similarly for exp(-2*t) and exp(-3*t), these would be the same as

 expTB = expTA .* expTA;
 expTC = expTB .* expTA;

I'm not sure if that will be much faster, but you have reduced the number of calls to exp by two-thirds...

[Abhinav]

@Hayes elements in lambda can be any positive values. Therefore, it is not possible to reduce the calls to exp your way.

[Geoff Hayes]

have you posted the actual function that you are using or just an example? We need to see the actual code so that we can provide some hints on how to improve the performance.

[Abhinav]

This is the actual function. Just the lambda values are arbitrary. I think that you are right about calls to exp taking a lot of time.

Answer 1

[OCDER]

This might be a variant of the xy-problem that @Stephen points out. http://xyproblem.info/

Solve X: you want to add the values for a function for all P, C, lambda values. Attempt with using function handles, called Solution Y.

Solve Y: you make a lot of function handles for all P, C, lambda sets. But this is too slow. You want a solution to make Solution Y faster, not seeing a faster solution for X.

%This is just an EXAMPLE. Check to make sure you get the right results
C = rand(1000, 1);
P = rand(1000, 1);
lambda = rand(1000, 1);
t = 1:10000;

%Solution Y approach  : use handles
Ysum = zeros(1, length(t));
for i = 1:length(C)
  fh = @(t) P(i)*exp(-lambda(i)*t) + C(i);
  Ysum = Ysum + fh(t);
end
t2 = toc

%Solution X approach  : use vectorized math
tic
Ysum = sum(bsxfun(@times, exp(-lambda*t), P) + C, 1);
t1 = toc  % ~5 times faster

I notice you like using repmat and find. Note that use of repmat is slower than bsxfun. Also, if possible, avoid using find if dealing with logical indexing. For instance:

function uh_fun=expconv(lambda)

      k=length(lambda);
  %     Lambda=repmat(lambda,k,1);
  %     temp_den=Lambda'-Lambda;   
      temp_den = lambda' - lambda; %faster
  %     ind=find(temp_den==0);
      ind = temp_den == 0; %stick with logical indexing

      if sum(ind(:)) > k % length(ind)>k %corrected this
          error('atleast two of the pdfs to be convolved are identical')
      else
  %         temp_den(temp_den==0)=1; 
          temp_den(ind) = 1; %no need to re-find temp_den == 0;
          den=prod(temp_den);
          nume=prod(lambda);
          C=nume./den;
          for i=1:k
              uh_fun{i}=@(t)exp(-lambda(i)*t)*C(i); %Do you reallly have to make all these handles?
          end
      end
  end

If this is still too slow, then you could try using parfor or mex functions.

[Abhinav]

Thanks for the suggestions! I will need some time to check-out all these.

[Abhinav]

@OCDER I have been able to pin-point the problem making a lot of calls to 'exp' function which seems unavoidable. As per your question if I really need to use a lot of function handle: I can actually avoid using function handles by supplying one a more argument 't' to the 'expconv' function and get the required results at desired values of t directly. Is that going to save some time? I will still be calling 'exp' function same number of times.

[OCDER]

If you look at the example, you could get away by calling exp once, but on a 2D matrix.

exp(-lambda*t), where lambda is a Mx1 matrix, and t is 1xN matrix.
so lambda * t is a MxN matrix.

This could save some time instead of trying to do a for loop and multip exp calls for each column of data.

The other question is what other function is calling your slow code "2,016,028" times? That's a lot... Do you often have the same input and output values? If so, you can cache your function outputs via memoize to prevent matlab from having to recalculate the output values.

https://www.mathworks.com/help/matlab/ref/memoize.html

This is what it does, assuming it's applied to a function myFun:

 B = myFun(1, 2, 3, 2, 3)  --> computing value B... --> B = 3;
 B = myFun(1, 2, 1, 2, 1)  --> computing value B... --> B = 1;
 B = myFun(1, 2, 3, 2, 3)  --> B = 3; %Same input/output as before

[Abhinav]

I have another 'unitHydrograph' in which the slow code is written. unitHydrograph has changing inputs. I have declared other useful variables which are calculated by unitHydrograph and remain fixed as 'persistent' so that I am not recomputing all those. Using matrix computation in lambda has not helped though.

Thanks a lot for your help and suggestions for improving my code!!