Weighted Polynomial Surface for 3D Points

Hi

I have an mx3 array of point cloud data (X,Y,Z) and a vector of weights for each point (mx1). I'm trying to create a Polynomial Surface (poly22) of the form:

f(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2

From this I need the value of the first coefficient, p00. At the moment I'm using the fit command as follows:

x = XYZ(:,1); % x-column

y = XYZ(:,2); % y-column

Z = XYZ(:,3); % elevation

fitsurface = fit([x,y],z,fittype('poly22'),'Weight',weight);

p00 = fitsurface(0,0);

This works well but is rather slow, considering that this is part of a function I'm working on which acts one-by-one on a point cloud of nearly 2 million points. If the points had no weight, I would use something like:

A = [ones(size(x)) x y x.^2 x.*y y.^2] \ z;

p00 = A(1);

This is much quicker than the 'fit' command; however, I have no idea how to make this work with weights involved. I have also tried the polyfitn function from the File Exchange - around 10% quicker for my dataset.

Any help finding a faster solution to this would be greatly appreciated.

Cheers, Jack

Answer 1

[John D'Errico]

Weighted fitting is simpler than it might seem. Of course, since I wrote polyfitn, why would you EVER want to use anything else? :)

% I'll pick some random weights here. Your weights
% probably make more sense than mine.
W = rand(size(x));

M = [ones(size(x)) x y x.^2 x.*y y.^2];
A = (diag(w)*M)\(w.*z);
p00 = A(1);

The idea is you simply multiply every line of the least squares problem by the corresponding weight. That scales the i'th residual by w(i).

I used diag to build a matrix to scale the rows of M there. If you had a HUGE number of points, that multiply will be less efficient. If you wanted speed, it would be better to use spdiags to build the diagonal matrix with the weights as sparse diagonal. Or, I could have used bsxfun here.

This will in fact be the fastest way to compute the vector of coefficients A, I think:

A = bsxfun(@times,w,M)\(w.*z);

It will certainly be faster than polyfitn, because this does no error checking, no model building, no computation of statistics on the coefficients, or things like R^2.

[Jack Williams]

John,

Many thanks for your swift response! This is around 5x quicker for my data than 'fit' and, as you note, also quicker than polyfitn. The use of bsxfun here is indeed the fastest option for my data.

One small query to follow up on this if that's ok. The following commands return slightly different results:

% Create points

xyz = randn(10000,3);

x = xyz(:,1);

y = xyz(:,2);

z = xyz(:,3);

% Create a set of weights - random here

weight = abs(rand(length(x),1));

% Apply the fit function

fitsurface = fit([x,y],z,fittype('poly22'),'Weight',weight);

p00 = fitsurface(0,0)

% Apply John's suggestion

M = [ones(size(x)) x y x.^2 x.*y y.^2];B = bsxfun(@times,weight,M)\(weight.*z);

p00 = B(1)

Having said this, the second result is within the 95% confidence bound of the fit command's coefficient. Is there a particular reason for this? Please excuse any mathematical naivety I may have displayed!

Thanks again and thank you also to Image Analyst for your response,

Jack

[John D'Errico]

For a random data set like yours, I get:

p00 =
    0.0078365

w = sqrt(weight);
M = [ones(size(x)) x y x.^2 x.*y y.^2];
B = bsxfun(@times,w,M)\(w.*z)
B =
  0.0078365
 -0.0082185
 -0.0050725
 -0.0096802
 -0.0053211
 -0.0014558

So B(1) matches here now, exactly so. The difference is in how to treat the weights.

By taking the square root of the weight vector, when the regression implicitly forms the sum of squares of the residuals, it multiplies each data point by that weight in the sum of squares. So we needed to take the square root of the weights to be consistent with how fit works.

Arguably the scheme that fit uses (with the square root of the weights) makes sense. A weight of 2 means that point has twice the importance of another point with a weight of 1, in terms of the weighted sum of squares of residuals. So a weight of 2 is equivalent to replicating that point in the set.

It is all in how you choose to treat the weights.

[Jack Williams]

John - many thanks indeed for your explanation. This makes sense to me now.

All the best,

Jack

[Chad Greene]

Note: Since John posted his answer, Matlab now has implicit expansion, meaning

A = (diag(w)*M)\(w.*z);

can be written simply as

A = (w.*M)\(w.*z);

and memory is no longer an issue by this method, even with a long vector of weights.

One point I'm not clear about, though, is whether perhaps it should actually be the square root of weights. I would love an explanation of whether we should or shouldn't take the square root of formal weights when calculating least squares polynomial fits.

Answer 2

[Image Analyst]

You might want to check this one out also: http://www.mathworks.com/matlabcentral/fileexchange/13719-2d-weighted-polynomial-fitting-and-evaluation. I haven't used it myself though.