how to speed up gaussian fitting

Greetings,

I am working on a piece of code that takes in some data and fits a gaussian with a constant offset to it. I first wrote a quick way to fit the gaussian that wasn't a true least squares method but it ran very fast ( could fit all my data in a matter of 10s of seconds). I now would like to go back and take the parameters that the quick fit gavem and use them to initialize a true least squares fit. What I've got it as follows :

           fo=fitoptions('Method','NonlinearLeastSquares',...
                       'Lower',[0 0 0 0],...
                       'Upper',[Inf Inf Inf Inf],...
                       'StartPoint', [ Ai Bi Ci Di]); % - setup the fitting options
            
           ft=fittype( 'a+b*exp(-(x-c)*(x-c)/(2*d*d))','options',fo); % make the fit type
                  
    
           [CURVE,GOF] = fit(x,y,ft); % - do the fit!
           

and I do this for each set of data points I want to fit. I have roughly 10000-20000 fits that I would like to do and this code is taking forever ( as in hours!) am I missing a way to speed this up? I could constrainthe upper and lower bounds slightly more, would that make a big difference?

For completeness I have the entire function below :

           function [xc, resid, yfit, A,diff] = recenterg(x, y, xc, width)
    % Fit y = c + m* x + A * exp(- (x-xc)^2 / (2 * width^2)) )
    %           c -> average background intensity
    %           m -> slope of background intensity
    %           A -> Amplitude of peak
    %           xc -> peak center
    %           width -> width of peak (essential psf)
    %
    % Input -> x,y -> vector of positions
    %           xc -> initial guess for peak position
    %           width ->  width of psf
    %
    %
    % Outputs ->
    %           xc -> refined peak center position
    %           resid -> rms of fitted intensities
    %           yfit -> fitted values of y
    %           A  -> fitted peak amplitude
    %
%     
    % Pick initial fit values so loop excutes at least once
    
    Niter = 0;
    Niter_max = 5; % Limit to 5 iterations
    xc_cut = 0.01; % Limit to fit precision of 0.01 pixels
    xc_min = min(x) + width; % Limits on center position
    xc_max = max(x) - width; 
    
    while (Niter<=Niter_max)  
        
        Niter = Niter + 1;
        
        % Perform linear regression which should be fast
        % of y = c + m* x + A * f + (A * xc_err) * df/dxc
        %        c = cbest(1); m = cbest(2) ; A = cbest(3); A * xc_err =     cbest(4)
        
        f = exp( - (x - xc).^2 / (2*width^2))/ width; % Peak function
        M = [ones(size(x)), x, f, f.*(x-xc)/width^2];
        cbest = (M'*M) \ M' * y;
        rmsd = std(y - M*cbest); % RMS error
        xc_err = cbest(4)/cbest(3);
    
        if (cbest(3)> 2*rmsd)&(abs(xc_err)<width)  
            % Peak signal to noise > 2 sigma
            % Correction in peak position < width of peak
            xc = xc + xc_err; % Update center position
            yfit = M*cbest;
            resid = sqrt(mean((y - yfit).^2));
            A = cbest(3);

            % End fit if convert to within precision of xc_cut
            if (abs(xc_err)<xc_cut)    
                Niter = Niter_max + 1;
            end
        else
            %xc = xc + 0.3 * xc_err;
            Niter = Niter_max +1 ; % Abort fit
            xc = NaN;
            resid = NaN; 
            yfit = NaN;
            A = cbest(3);
        end
    
% ok ,ran with the quick and dirty f it, now take those quick
        %values and use them to initialize a true Gaussian fit
        
        %****** START HERE GAUSSIAN FITTING PROCEDURE
        % first check if the value was NaN or 0 if it was then leave it as
        % it
        if (isnan(xc)) % when the quick fit doesnt work the gaussin fit won't either so give them NaN 
            
            xc = NaN;
            resid = NaN; 
            yfit = NaN;
            A = cbest(3);
            diff=NaN;
        else
            Ai=  (sum(y(1:4))+sum(y(end-3:end)))./8;  % a rough guess of the background 
            Bi=   A; % - guess for the amplitude
            Ci=   xc; % - from the quick fit
            Di=   1.9 ; % the width as a first guess
            
           fo=fitoptions('Method','NonlinearLeastSquares',...
                       'Lower',[0 0 0 0],...
                       'Upper',[Inf Inf Inf Inf],...
                       'StartPoint', [ Ai Bi Ci Di]); % - setup the fitting options
            
           ft=fittype( 'a+b*exp(-(x-c)*(x-c)/(2*d*d))','options',fo); % make the fit type
                  
           [CURVE,GOF] = fit(x,y,ft); % - do the fit!
           diff=abs( xc - CURVE.c); % - to keep track of the difference between quick method and real fit
           xc = CURVE.c; % - the center position of the fit   
           yfit = CURVE.b; %- don't really need this but might be useful later
           resid = CURVE.d; % - the std of fit
           A = CURVE.b ; %- the fit amplitude
          
            
        end

Answer 1

[Matt J]

https://www.mathworks.com/matlabcentral/answers/431697-make-curve-fitting-faster#answer_348931