Faster method for polyfit along 3rd dimension of a large 3D matrix ?

Hello, I have the same problem as mentiond in the link http://www.mathworks.co.jp/matlabcentral/answers/45680-faster-method-for-polyfit-along-3rd-dimension-of-a-large-3d-matrix

I tried to understand the answer but I cant seem to understand the mathematics behind it, I dont know why does it give the right answer.

Also, in my program I do exponential fit along the 3rd dimension of a large 3d matrix by two for loops.

My code:

dataCube = rand(480,640,500);
x = rand(500,1);
sizeCube = size(dataCube);
exCube = zeros(sizeCube);
for ii = 1:sizeCube(1);
   for iii = 1:sizeCube(2);
        s(1,:) = dataCube(ii,iii,:);
        signal = s';
        a = fit(x,signal,'exp2')
        y = feval(a,x);
        exCube(ii,iii,:) = y';
     end
  end

This code takes several hours to finish running.

Is there a fast method for exp fit along 3d dimesnion for large 3d data matrix ?

Any help would be greatly appreciated.

Answer 1

[Matt J]

Because a polynomial has a linear dependence on the unknown coefficients, the solution for a polynomial fit is a closed-form linear function of the input data, and therefore can be expressed as a matrix multiplication. You can't accelerate the fitting for exp2 the same way as for a polynomial, because the least squares fit for exp2 is neither linear nor closed-form as a function of the input signal.

Nevertheless, there are tricks you can use, depending on your dataCube. For example, if neighboring signals in your dataCube are always very similar to each other, you can initialize the fit of one signal in the set with the solution from the previous signal. This reduces the number of required iterations in the fit.

[Matt J]

You can also try using FMINSPLEAS on the File Exchange instead of the Curve Fitting Toolbox. FMINSPLEAS can take advantage of the fact that the exp2 model is linear in 2 of the 4 unknown parameters. I don't know whether fit() does the same.

[John D'Errico]

(I'm not sure that fit does use a partitioned scheme like fminspleas. It should, and it is quite doable. In fact, I once wrote a nonlinear regression toolbox with a complete gui in front, years before the curvefitting toolbox existed. It used a variable partitioning scheme on all models.)

Anyway, you could accelerate this to make it run like blazes, with a semi-custom MATLAB code.

Use the fminspleas idea on each problem. It looks like exp2 is a two term exponential with no constant term, so there are two nonlinear parameters and two intrinsically linear parameters. So if you have n problems, there are 2*n nonlinear parameters, 2*n linear ones.

So solve this as a batched problem, with 2*n unknowns, using perhaps lsqnonlin or lsqcurvefit. Specify a Jacobian pattern as a block diagonal matrix, so that lsqcurvefit knows that the subproblems are disjoint in reality. If there are n sub-problems, with m data points in each sub-problem, then there will be n*m data points in total for the problem, with 2*n unknowns.

Essentially, the idea is to solve ONE massive problem, but help lsqcurvefit by letting it know that it is really a system of sub-problems.

The above scheme will be eminently doable, as long as you choose decent starting values for the nonlinear parameters. Since you need no starting value at all for the linear ones, the problem gets much simpler, and far more robust. My guess is you can expect perhaps a 100-1 time reduction for this scheme, and perhaps more, since you were solving n problems with 4 unknowns for each subproblems before, and you would have been doing it with a loop.

I discuss some of the issues of how to build a batched problem in my optimization tips and tricks document on the File Exchange.

[Lucia]

Thank you very much, I looked at the file but didn't quite understood it (lots of stuff there).

Can you please give me a coded example of how to use it in my code ?

Thanks in advance

[Matt J]

Here is an example exp2 fitting using FMINSPLEAS,

        x = linspace(0,3,100);
        y = 4*exp(-3*x) - 5*exp(-2*x);

        funlist = {@(cd,x) exp(cd(1)*x),...
                  @(cd,x) exp((cd(1)+cd(2))*x)};

        cd0 = [-1,-4]; %initial guess

        options = optimset('disp','iter','TolX',1e-20,'TolFun',1e-20);
        [cd,ab] = fminspleas(funlist,cd0,x,y,[],[],[],options);

        a=ab(1),
        b=ab(2),
        c=cd(1),
        d=cd(1)+cd(2),