Too much calculus time on 'lsqcurvefit' + 'fsolve'

Good afternoon everybody,

I am working with a code that combines 'fsolve' and 'lsqcurvefit' to calculate the parameters, k(1), k(2) and k(3), that better fits my function, called 'Iteor', to my experimental data, (texp, Iexp):

rteor=@(k,r) fsolve(@(r) arrayfun(@(R,T) FC.*integral(@(x) exp(-a./(x.*log(x))), 1, R)-(T-(.001.*k(1))),r,texp), inpts, options);
Iteor=@(k,r) ((.01.*k(2)./FC).*(rteor(k,texp)).*exp(a./(rteor(k,texp)).*log(rteor(k,texp))))+(Vo./(100.*k(3).*log(rteor(k,texp))));
k = lsqcurvefit(Iteor, x0, texp, Iexp,[],[],options);

The problem is my experimental data are 22000 points large, and the calculus takes too much time. I'm working on a Intel Core i7 4.00 GHz and 32.0 GB. With FORTRAN the process is much faster... Is there any way to reduce the calculus time without eliminating initial experimental points ?

Thank you so much !!

Answer 1

[OCDER]

Your Iteor is causing a lot of issues since it's calling rteor 4 times, meaning it's doing repeated fsolve and integral - both of which are computationally expensive.

To translate what's happening:

rteor = @(k, r) HEAVY_MATH_STUFF;
Iteor = @(k, r) HEAVY_MATH_STUFF .* HEAVY_MATH_STUFF + log(HEAVY_MATH_STUFF) ./ HEAVY_MATH_STUFF;

To fix, create a function that does the heavy math only once.

function y = Iteor(k, r, ...)
   b = rteor(k, r, ...);
   y = b .* b + log(b) ./ b;  %b was only computed once!

I wish I could give you the corrected code, but that equation is very complex with all the nested anonymous functions... I'll end up messing it up. Overall, prevent recalculating fsolve and integral. Hope this helps to speed things up a little.

[Sergio Quesada]

Thanks a lot, OCDER. I have already written an output function for the 'lsqcurvefit' routine to create and write up the results found. Something like this:

function stop=outfun(x,optimValues,state) 
persistent fID_I_vs_texp
stop = false;
      switch state    
        case 'init' 
          fID_I_vs_texp=fopen([...
        case 'iter' 
          disp({...)
              if optimValues.resnorm<initialerror; 
                 fprintf(...  
             initialerror=optimValues.resnorm; 
              end
       case 'done'
           fclose ('all');
    end

So, as you say, Should I write an additional order in the case 'iter' containing your commands, in order to reevaluate 'rteor' only once after each iteration? How can I do it?

Thanks!!

[OCDER]

I'm not quite sure where this fits in with the other codes. Where is the iteration?

[Sergio Quesada]

Outfun is an output function of my lsqcurvefit. The commands after "state 'iter'" are executed At the beggining of each iteration

[OCDER]

Sorry, I can't seem to figure out where state iter is executed in the following statement, or locate the call to outfun. Is this "iteration" the one used by lsqcurvefit internally? It'll be very slow then as every time it tries to fit your data, it's printing the results too - HDD read/write access is much slower than RAM access:

rteor=@(k,r) fsolve(@(r) arrayfun(@(R,T) FC.*integral(@(x) exp(-a./(x.*log(x))), 1, R)-(T-(.001.*k(1))),r,texp), inpts, options);
Iteor=@(k,r) ((.01.*k(2)./FC).*(rteor(k,texp)).*exp(a./(rteor(k,texp)).*log(rteor(k,texp))))+(Vo./(100.*k(3).*log(rteor(k,texp))));
k = lsqcurvefit(Iteor, x0, texp, Iexp,[],[],options); %%? Where is outfun?

Since this seems like a different question, you could try posting another question to the forum like:

How to save outputs of each iteration of lsqcurvefit efficiently?

Perhaps one of the other forum users might know a clever trick.

[Sergio Quesada]

Sorry, OCDER, I didn't explain myself properly. On my lsqcurvefit, y selected the output function like this:

options = optimoptions(@lsqcurvefit,'Algorithm','levenberg-marquardt','Display','iter-detailed','OutputFcn',@outfun_3);
k = lsqcurvefit(Iteor, x0, texp, Iexp,LB,UB,options);

, where outfun3 is my OutPutFcn containing the switch. I have written my outputfunction and the code containing lsqcurvefit on different .m files. Any idea now ?

Thanks again !

[OCDER]

I see. The output function you have is only going to save the results of iteration to a file and does not summon your slow code, rteor, so no issue there with fsolve or integral being summoned by your output function.

However, any function that does fprintf per iteration could be slow. This is why if speed is important, most codes disable fprintf. A workaround is to save the data you want to print into a memory-stored variable instead, and print it to file in blocks. An example:

function stop = outfun(x,optimValues,state)
persistent fID_I_vs_texp StoredData J;

stop = false;
switch state
    case 'init'
        fID_I_vs_texp = fopen('filename.txt');
        StoredData = cell(1000, 1); %Data to store your print data
        J = 0;
    case 'iter'
        J = J + 1;
        if J > 1000 
            fprintf(fID_I_vs_texp, '%s\n', StoredData{:}); %Save your data in 1000-line blocks
            J = 0;
        end
        StoredData{J} = sprintf('iter = %d, fval = %d', optimValues.iteration, optimValues.fval); 
    case 'done'
        fprintf(fID_I_vs_texp, '%s\n', StoredData{1:J}); %Save the last bit of data
        fclose ('all');
end

[OCDER]

And your Iteor function should be looking something like this. I can't test it out since I don't have license to use lsqcurvefit.

function F = Iteor(k, texp, R, T, FC, a, r, inpts, options)
G = rteor(k, texp); %Do this once
LogG = log(G);      %Do this once
F = ((.01.*k(2)./FC).*G.*exp(a./(G).*LogG))+(Vo./(100.*k(3).*LogG));

function G = rteor(k, texp, R, T, FC, a, r, inpts, options)
G = fsolve(@(r) arrayfun(@(R,T) FC.*integral(@(x) exp(-a./(x.*log(x))), 1, R)-(T-(.001.*k(1))),r,texp), inpts, options);

Use lsqcurvefit something like this:

 x = lsqcurvefit(@(k,texp) objfun(k, texp, R, T, FC, a, r, inpts), k0, texp0, ydata, ...);

WARNING! your options structure that has your output function is going to print for both fsolve AND lsqcurvefit iterations. Double printing, especially with a persistent variable, will give you improper results and also will be very slow. To see if this outfun is the slow function, try running your code without it first.

[Sergio Quesada]

Thank you so much, OCDER !

[OCDER]

You're welcome!