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

106 views (last 30 days)
Justin on 9 Aug 2012
Edited: Jan on 9 Apr 2019
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 on 9 Aug 2012
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 on 14 Aug 2012
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

Teja Muppirala on 10 Aug 2012
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 on 14 Aug 2012
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 on 18 Dec 2017
+1, Teja, what a powerful solution!

Martin Offterdinger on 9 Apr 2019
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

#### 1 Comment

Jan on 9 Apr 2019