Fastest (parallelized, maybe) way to run exponential fit function in a big matrix of data

Dear experts,

I am trying to ensure the most optimized syntax to run a part of my code. In summary, I need to perform an exponential fit of a series with 6 points, but ~2 million times (independent samples). In other words, I have a reshaped matrix 6x2,000,000.

This is my code right now (each loop iteration takes 0.02 seconds what makes this process prohibitive):

Xax = [30;60;90;120;150]; % The constant X series.
f = @(b,Xax) b(1).*exp(b(2).*Xax);  % Exponential model
% LoopMap == my 6x2,000,000 matrix
for k = 1:size(LoopMap,2) %parfor? Any GPU model (if this is and optimal example)
    bet = fminsearch(@(b) norm(LoopMap(:,k) - f(b,Xax)),[0;0]); 
    CoefAtmp(1,k) = abs(1/bet(2)); 
end

My PC is a standard machine (8th gen Core i7, 32Gb ram) with a not exceptional graphic card.

Thank you all in advance.

[Matt J]

Do you have the Parallel Computing Toolbox? Is parfor an option?

Answer 1

[Matt J]

Additional speed-up should also be possible by reducing your problem to a 1-variable estimation and removing subsref operations from your objective function.

Xax = [30;60;90;120;150]; % The constant X series.
B2=nan(1,size(LoopMap,2));

Initial=Xa.^[0,1]\log(LoopMap);
Initial=Initial(2,:);

for k = 1:size(LoopMap,2) %parfor? Any GPU model (if this is and optimal example)
    
    y=LoopMap(:,k);
    
    B2(k) = fminsearch(   @(b2) objective(b2,Xax,y) , Initial(k)); 
    
end
CoefAtmp = abs(1./B2); 


function  cost=objective(b2,Xax,y)

  ex=exp(b2*Xa);
  
  b1=ex\y;
  
  cost=norm(b1*ex-y);


end

[Brunno Machado de Campos]

Thanks again, this solution is 200 times faster and just 7 times slower than the linear fit on the log transformed data, very reasonable!

Answer 2

[Matt J]

If it's acceptable to you, a log-linear least squares fit can be done very fast and without loops

A=Xax.^[0,1];

Bet=A\log(LoopMap);

CoefAtmp=1./abs(Bet(2,:));

[Brunno Machado de Campos]

Thank you Matt. This is exactly the code I am running in the last version of my code and indeed provide very similar results, but, the academic community expects the exponential fit. I am evaluating the possibility to implement just to avoid predictable and easy criticism. Thank you.

[Matt J]

Both methods are exponential fits.

Do you mean you must have a linear least squares objective? If so, intiailizing fminsearch with the log-linear solution will probably speed things up.

However, you should first check that the linear and the loglinear methods give significantly different results over a small but representative subset of your LoopMap(:,k). If not, you can use that to defend the loglinear method to the academic community.

[Brunno Machado de Campos]

Right! Sure, thank you. Since a simple exponential fits well the model, both option will give very similar answeres.

So (just to documment):

For an example were sig = [527.32;398.95;284.22;197.562;155.01];

1)

Xax = [30;60;90;120;150]; % The constant X series.
f = @(b,Xax) b(1).*exp(b(2).*Xax);
bet = fminsearch(@(b) norm(sig - f(b,Xax)),[0;0]); 
CoefAtmp1 = abs(1/bet(2)); 

0.02 seconds

2)

Xax = [30;60;90;120;150]; % The constant X series.
M = bsxfun(@power,Xax,0:1);
c2 = M\log(sig);
CoefAtmp2 = abs(1./c2(2));

0.000002 seconds

Nothing more to do here, thanks again!