Multiple use of polyfit - could I get it faster?

I have a more general problem in the sense that I am supposed to do polyfit more than 10000 times on a given problem. Is there any way I can improve the speed of this instead of looping more than 10000 times through all of this? I do assume that the polynoms are all of the same degree.

Basically, I want to do matrixP = polyfit(matrixX, matrixY, n)

[Bruno Luong]

There might be some way if the numbers of data remain unchanged. Can we assume that?

Answer 1

[Jan]

If you are sure, that your X-values are properly scaled (perferred: mean(x)==0, std(x)==1) and do not need the 2nd and 3rd output of POLYFIT, you can omit the nice and secure error checks of POLYFIT:

P = polyfitM(matrixX, matrixY, n)

[sx1, sx2] = size(x);
P = zeros(sx1, n+1);

for Index = 1:sx1
  x = matrixX(Index, :);
  y = matrixY(Index, :);
  x = x(:);
  y = y(:);
  q = ones(sx1, 1);
  V(:, n+1) = q;
  for j = n:-1:1
    q = q .* x;
    v(:, j) = q;
  end
  [Q, R] = qr(V, 0);
  P(Index, :) = transpose(R \ (transpose(Q) * y));
end

If your X-values are equidistant, you could save 99.9% processing time by calculating the Vandermonde matrix V and the QR decomposition once only. But even the reduced error checks in the posted method can run a factor of 2 to 3 faster compared to calling POLYFIT in a loop - please check this on your computer.

Answer 2

[Matt Tearle]

In the manner of Jan's answer, just stripping down to a brute-force least-squares fit can give a lot of speed-up. I'm assuming that you have a matrix of x and a matrix of y, with each column representing a different data set. And also, therefore, assuming that each data set has a different set of x values. Here's a function to compare them:

function [c1,c2] = multpolyfit(x,y,n)

m = size(x,2);
c1 = zeros(n+1,m);
tic
for k=1:m
    c1(:,k) = polyfit(x(:,k),y(:,k),n)';
end
toc

c2 = zeros(n+1,m);
tic
for k = 1:m
    M = repmat(x(:,k),1,n+1);
    M = bsxfun(@power,M,0:n);
    c2(:,k) = M\y(:,k);
end
toc

My results:

>> x = rand(20,10000);
>> y = pi + x - 0.5*x.^2 + 0.1*rand(size(x));
>> [c1,c2] = multpolyfit(x,y,2);
Elapsed time is 9.531862 seconds.
Elapsed time is 0.610350 seconds.

[Ian]

Nice function. Might you also be able to modify it to output the R^2 value (ie the coefficient of determination)?

Would be most grateful for that addition. Can't seem to figure out how to add it without doing the calculation for every fit. Thanks.

[Jan]

@Ian: Do the R values created in my solution help?

[Matt Tearle]

@Ian: something like this?

function [c,R2] = multpolyfit(x,y,n)

m = size(x,2);

c = zeros(n+1,m);

r = zeros(size(y));

for k = 1:m

M = repmat(x(:,k),1,n+1);

M = bsxfun(@power,M,0:n);

c(:,k) = M\y(:,k);

r(:,k) = M*c(:,k)-y(:,k);

end

sserr = sum(r.^2);

sstot = sum(bsxfun(@minus,y,mean(y)).^2);

R2 = 1 - sserr./sstot;

It seems to be about 25% faster than doing the individual calculation every time through the loop. Caveat: I *think* the formula/implementation is correct... but I wouldn't bet my house on it.

Answer 3

[Roy]

Thanks for the tips - works nicely now with less than 5% of the CPU-time.

[Jan]

Does this mean, that you accept one of these answers?