Speed up process of passing each row of large array to external function

I have a very large array consisting of 10,000 rows and 1,927 columns. Each row is a photoluminescence spectrum that I am passing from AppDesigner to an external function to fit multiple Gaussians. I have tried using Image Analyst's fit_multiple_gaussians.m and Tom O'Haver's findpeakfits.m, which both work very well. However, to no one's surprise, calling either of these functions 10,000 times in a for-loop is very slow for spectra that contain 1,927 data points. Because each row of the array is an independent spectrum, is there an efficient way to parallelize the function call, or possibly perform these operations on a GPU without having to decrease the number of data points in each spectrum? Below is the code snippet of the for-loop:

x = lambdas % 1 x 1927 element vector
gaussians = zeros(length(int(:,1)),numGaussians,length(lambda)); %pre-allocate size of gaussians array

%Loop over the length of the array and pass each row to the function
for i = 1 : length(int(:,1))
    
    %int is a 10,000 x 1,927 array containining intensity data. 10,000
    %independent spectra with 1,927 data points per spectrum.
    y = int(i,:);
    
    FitResults = []; %clear previous result. FitResults is a small Nx5 array, where N is determined by the number of Gaussians needed to fit the spectrum
    FitResults = findpeaksfit(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype,peakshape,extra,NumTrials,autozero,fixedparameters,plots);
   
    %Extract Gaussian peaks and widths from function output
    for j = 1 : size(FitResults,1)
        peakPosition = FitResults(j,2);
        width = FitResults(j,4);
        %pass Gaussian parameters to another function to rebuild Gaussian
        gaussians(i,j,:) = gaussian(app,lambda,peakPosition,width);
    end
    
end

function g = gaussian(app, x, peakPosition, width)
    g = exp(-((x - peakPosition) ./ (0.60056120439323 .* width)) .^ 2);
end

[Dyuman Joshi]

On a cursory glance -

There is no need of this line of code, get rid of it.

FitResults = []; %clear previous result. FitResults is a small Nx5 array, where N is determined by the number of Gaussians needed to fit the spectrum

Preallocate the gaussians array.

You have defined 4 inputs to the function gaussian, but the input app is not used. What is the correct function definition then?

Also, do you have the symbolic toolbox?

[eherrm]

Hi Dyuman, thank you for taking the time. I edited the code snippet above just to show that gaussians has already been pre-allocated. While pre-allocation increases speed, the total run-time of the code is still very long.

Also, the call to gaussian() is defined within AppDesigner. I get errors if I do not include app in the function call. However, removing this call to gaussian() entirely does not significantly increase speeds. A bare-bones code snippet is below, which is still very slow. After trying multiple fitting functions, I am more curious to know if there is a way to split up this array so I can run multiple processes at the same time or on a gpu.

I do not have the symbolic toolbox.

x = lambdas % 1 x 1927 element vector

%Loop over the length of the array and pass each row to the function
for i = 1 : length(int(:,1))
    
    %int is a 10,000 x 1,927 array containining intensity data. 10,000
    %independent spectra with 1,927 data points per spectrum.
    y = int(i,:);
    
    FitResults = findpeaksfit(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype,peakshape,extra,NumTrials,autozero,fixedparameters,plots);
   
    %Extract Gaussian peaks and widths from function output
    for j = 1 : size(FitResults,1)
        peakPosition = FitResults(j,2);
        width = FitResults(j,4);
    end
    
end

[dpb]

"...the call to gaussian() is defined within AppDesigner."

Ergo, we can presume it is a callback function, then???

Have you compared IA's solution to that of fit in the Curve Fitting TB? It has builtin model to support up to eight (8) term Gaussian. Whether it might be better optimized than IA's using the builtin basic fminsearch I don't know. I doubt IA really thought too much about trying to do 10,000 spectra at a time... :)

I've never messed with the Parallel TB nor tried seriously to use the GPU so I dunno about attacking the problem in those directions. Undoubtedly it should be possible if you've the hardware to support it...

[dpb]

How many peaks do you typically have in the model?

A local comparison of the demo in @Image Analyst's script showed that without setting initial guesses, the 'Gauss6' builtin model was about twice as fast as the computation section in the script, but the fit wasn't quite as good for the last peak.

Setting the options object to fit to also include starting values for the mean, width, parameters, however, caused the run time to almost quadruple or roughly twice as long so that's a loser.

I'd presume your spectra are better behaved than the demo there, but that's purely a guess...but, it's probably worth investigating given the number of these you apparently have to do.

[eherrm]

I actually haven't tried the builtin fit model...probably should have done that first!

The data shows 0-4 guassian peaks, so I have been using a 4-gaussian model in @Image Analyst's script. With @Catalytic's suggestion of using parfor, it looks like the fitting should finish after a total time of 1.5 hours (it's currently running), which is significantly faster. Once this run finishes, I will attempt your suggestion of using the builtin 'Gauss4' model, in addition to parfor.

The clarity of the spectra contained in one dataset varies dramatically, with varying numbers of peaks, intensities, noise levels, etc. I'm building this app to fit multiple datasets as I collect them, so I'll definitely investigate further to increase efficiency.

Thank you for taking the time, I really appreciate it! I'll update you once I run the 'Gauss4' model.

[dpb]

"The clarity of the spectra contained in one dataset varies dramatically, with varying numbers of peaks, intensities, noise levels, etc. ..."

I can imagine -- I developed alogorithms for gamma spec for an online coal elemental S analyzer some 40 years ago (in Visual Basic, believe it or not!) without any such helping features...implemented own Marquardt solver in conjunction with the linear quadratic background correction...

There, however, it took some 5 minutes to collect enough counts so that it was only having to do that calc that often. Unfortunately, with the technology then and no cryogenics, that was using NaI crystal and so the resolution was so poor the analyzer itself was a device ahead of its time for power plant applications. In the time since, it's night and day although in the end for combustion control, simply the elemental S isn't that useful because the coal firing qualities are so highly dependent upon the particular sulfur compounds present more than just the amount of total S.

Answer 1

[Catalytic]

This reduces it to one loop, which should work with parfor -

x = lambdas; % 1 x 1927 element vector


intt=int';

logg=nan(numGaussians,length(lambda), size(int,1));  %pre-allocate log of gaussians


%Loop over the array and pass each COLUMN to the function
for i = 1 : size(int,1)

   FitResults = findpeaksfit(x,intt(:,i),   SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype,peakshape,extra,NumTrials,autozero,fixedparameters,plots);
 
   peakPosition=FitResults(:,2);
   width=FitResults(:,4);
   
   logg(:,:,i)=(x - peakPosition) ./ width ; 
    
end

gaussians=exp(-(logg/0.60056120439323).^2);

gaussians=permute(gaussians,[3,1,2]);

[eherrm]

This speeds up the calculation significantly, thank you!