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

I am currently working with large data cubes in the form of an MxNxP matrix. The P-dimension represents the signal at each (m,n) pixel. I must obtain a Z-order polynomial fit(where Z varies depending on the situation) for each signal in the data cube. Currently, I utilize a "for" loop to obtain the signal at each pixel, obtain the polynomial coefficients, then calculate the fitted curve. The code I use is fundamentally similar to the following:

dataCube = rand(1000,1000,300);
x = rand(300,1);

sizeCube = size(dataCube);

polyCube = zeros(sizeCube);

for ii = 1:sizeCube(1);
   for iii = 1:sizeCube(2);

        signal = squeeze(dataCube(ii,iii,:));

        a = polyfit(x,signal,z)
        y = polyval(a,x);

        polyCube(ii,iii,:) = y;

     end
  end

Because of the quantity of iterations in the for loop, this operation takes a considerable amount of time for each data cube. Is there a faster way to obtain the polynomial fitting, without having to resort to the iterative process I use here. Perhaps, something similar to the filter function where you can apply the filter to a specific dimension of a matrix, rather than extracting each signal?

filteredCube = filter(b,a,dataCube,[],3)

Thanks, Justin

[Walter Roberson]

Have you considered re-ordering your data, at least during the processing, so that the dimension you are fitting over is the first dimension? Access (and assignment) over the first dimension is faster.

[Justin]

Walter,

I attempted reordering the matrix and performing the process as above, over the 1st dimension. There was an incremental improvement in speed, but not nearly the type of advance I was hoping I could achieve.

Check out the answer by Teja below. Its exactly what I was hoping for and works perfectly.

Thanks for the input! Justin

Answer 1

[Teja Muppirala]

This can be accomplished in a fraction of the time with some matrix operations.

dataCube  = rand(100,100,300);
sizeCube = size(dataCube);
x = rand(300,1);
z = 3;

V = bsxfun(@power,x,0:z);
M = V*pinv(V);
polyCube = M*reshape(permute(dataCube,[3 1 2]),sizeCube(3),[]);
polyCube = reshape(polyCube,[sizeCube(3) sizeCube(1) sizeCube(2)]);
polyCube = permute(polyCube,[2 3 1]);

[Justin]

Teja,

Thank you! This is a slick solution that drastically reduces the processing time. It went from ~50 minutes to process a 1920x1040x301 image_cube (@ 0.0015 seconds per loop), to about 15 seconds!

Thank you for the help! Justin

[Jan]

+1, Teja, what a powerful solution!

[Pavel Psota]

+1, cool, indeed! Could the fitting coefficients be obtained from your solution?

[Jan]

A small improvement is to avoid the expensive power operation:

% Replace:
V = bsxfun(@power,x,0:z);
% by:
V = [ones(300, 1), cumprod(repmat(x, 1, z), 2)];

Answer 2

[Martin Offterdinger]

Dear Teja,

I am having a similar problem- actually a simpler one even. I have the same array, but I always need to fit a first-order polynom (linear, z=1 in your code). Is it possible to get the coefficients of the linear fit (p1,p2) from your solution as well?

Thanks,

Martin

[Jan]

@Martin: Please do not attach a new question in the section for answer. Open a new thread instead and remove this pseudo-answer. Including a link to this thread is a good idea. Thanks.

What's wrong with setting z=1? Which array is "the same" and why do you need to determine the fit multiple times for the same array? (Please explain this in your new question...)