Fit - Nonlinear Regression - can this run faster? Python >20x faster?

Hi,

I'm looking for a way to speed up fitting an exponential model to a stack of eleven double images. I've been at it for days but can't seem to find a way to speed things up. Even using parfor this takes 20 minutes on a downsampled image stack (100x100), which is WAY slowwer than single thread of python (6 minutes per 1024x1024 pixel stack!). I was wondering if somebody could help me out speeding things up a bit - I guess I must clearly be missing something...

I load all 11 images into an array and permute it in order to be able to access the fitting dimension right.

imgdata2 = permute(imgdata,[3,1,2]);
size(imgdata2)
ans =

        11        1024          1024

The model I want to fit is:

fun = @(j, x)  j(1).*exp(x./0.0257)+j(2).*exp(x./(2*0.0257))- j(3);
pguess = [100e-15, 10e-9, 0.04];

Y is given and the same for all fits. A simple fit now looks like this:

   x = imgdata2(:,i,k)./1000;
  out = fitnlm(x,y,fun, pguess);

Now the fastest way I could figure out to fit all 1024*1024 pixels of the image stacks was using parfor (yes, I pre-start the workers). For this I had to create a dummy cell array and then unwrap the results afterwards in another loop.

%%preallocate
    fitresults_wrapper = cell(1,(size(imgdata2,2)*size(imgdata2,3)));

%%fitting
parfor i = 1:(size(imgdata2,2)*size(imgdata2,3));   % for all 1024 by 1024 'pixel-stacks'
        x = imgdata2(:,i)./1000;
        if sum(isnan(x))== 0
                fitresults_wrapper{i} = nlinfit(x,y,fun, pguess);
        end;
end;

Now. Is there no way to speed things up? If this always takes a couple of minutes it's pretty much useles. A simple python loop can do the full images in 6 minutes using just one thread (curve_fit(func, x, y, [50e-12,1e-6, 0.04])) I was unable to find a way to vectorize the whole thing, which is usually the first thing I do in matlab in order to speed things up. Is there a way to vectorize fitting here?

Any help or hints appreciated Arnold

Answer 1

[Ilya]

I am no expert on image processing, but based on your code samples, you don't need to use fitnlm. When you do this

x = imgdata2(:,i,k)./1000;

or this

x = imgdata2(:,i)./1000;

you get a vector. Just form a matrix with two columns by putting

X = [exp(x./0.0257) exp(x./(2*0.0257))];

and pass this matrix to a linear utility such as say fitlm to fit for the two linear coefficients, j(1) and j(2), and intercept, -j(3). This is going to be much faster.