3D Least squares polynomial fit in x and y

Fit a two dimensional f(x,y) polynomial to sampled x,y,z data triplets

Description

Often, measured data is comprised of N sampled values of z, evaluated at N locations (x,y). With this function, you can calculate the coefficients of the best-fit x,y polynomial using a linear least squares approximation.

You can use this function if you have a set of N data triplets x,y,z, and you want to find a polynomial f(x,y) of a specific form (i.e. you know the terms you want to include (e.g. x^2, xy^3, constant, x^-3, etc.) in your fitting polynomial.

• MATLAB release: MATLAB 7.7 (R2008b)

[John D'Errico]

Sorry, but this is poorly written, and uses poor numerical methods to solve the problem.

Why do I say it is poorly written? It uses loops, building its arrays one element at a time, without benefit of any preallocation.

The poor numerical methods part comes from its use of the normal equations to solve the linear least squares problem. When the system is poorly conditioned, this solver will often result in inaccurate estimates of the parameters, but you will never be informed of that fact. The user will just get crap for results, and have no idea why, or even that they got a poor set of estimates.

No system scaling is done to avoid numerical problems.

There is no error checking done at all here.

There is no H1 line, at least not a meaningful one. The H1 line is the VERY FIRST line if your help. It is what MATLAB uses when you use the lookfor command. The lookfor command allows you to find a function when you don't remember the name of the function that you downloaded last year. Sorry, but "least_square_polyfit_xyz.m" is not the most memorable name in the world. So the H1 line should be a single line of comment that includes some reasonable keywords to search for. The author's name and e-mail address do not qualify as such. They may be useful information, but put them somewhere else in the code!

Finally, this code does not return parameter variances, or R^2, or any other piece of information that might be remotely useful in your analysis.

Look for polyfitn, a code that does all this does and does it better.

http://www.mathworks.com/matlabcentral/fileexchange/10065

[Thomas]

John makes excellent points. Polyfitn has lots of error checking, and 300+ lines of what looks like well written code. By comparison, this function delivers the core 8 line algorithm only. Some of my colleagues have been using nonlinear LSQ solvers when fitting multidimensional polynomials, and I wanted to provide an example of the linear alternative. this might not be the right venue, and I encourage those looking for all bells and whistles to use polyfitn.

[John D'Errico]

But if all you want is to show how to build a simple code, it is still easily done without using the normal equations. For example, given column vectors x, y, z, here is a better solution using only ONE line of code:

coef = (bsxfun(@power,x,m(:,1)').*bsxfun(@power,y,m(:,2)'))\z;

Better is to include various error checks to make the code friendly, so that one need not worry about whether the user has actually provided a set of column vectors. At least, make it robust, so either row or column vectors will suffice. Even this can be done in one line!

coef=(bsxfun(@power,x(:),m(:,1)').*bsxfun(@power,y(:),m(:,2)'))\z(:);

Simple scalings can be easily added too, even returning parameter variance estimates in only a few lines.

So even without all the bells and whistles, this could/should have been done better. I'll repeat that the normal equations are a poor way to solve linear least squares problems if you are EVER worried that those equations may be poorly conditioned. And polynomial problems are a common source of ill-conditioned linear systems.

AVOID the normal equations. Too often they are taught as the way to solve linear least squares problems by someone who does not understand why they are bad. This propagates to their students. You will even find them in books, but just because something is in a book means nothing. You don't need to pass a test on numerical analysis to write a book.