polyfitn

Dennis Wouters,
Thanks for providing this fix. It works!

Despite my having written polyfitn, high order polynomials (especially in multiple dimensions) are rarely a good substitute for a nonlinear model that is chosen to have the proper behavior. That is, IF you have such a model. Those polynomials tend to do nasty things to you, especially in extrapolation. And unless you have sufficient data, then the model may behave strangely even inside your data. I call this the problem of interpolation.

Basically, you need to be very careful with high order polynomial fits.

By the way, nlinfit DOES handle multiple independent variables. You suggest it does not, but there is no constraint on that. The problem is choosing an intelligent model in multiple dimensions. This is EXTREMELY difficult in general. One good reason is that a functional form like w = f(x,y,z) is so difficult to visualize.

Personally, my preference has always been to move to a spline-like model for these problems. I am often unable to pull a model out of thin air, and I would prefer to avoid a high order polynomial. My own gridfit is one simple example, with 2 independent variables.

Hi Karan,

You've fallen into a common fallacy in regression modeling. It goes like this: "If n is good, then n+1 is better, and n+2 must be really wonderful."

The fallacy starts out easily. We try a linear model on our data. It works, all right, but not very well. After all, must problems are only pretty locally linear. So we try a quadratic model. It fits better. The lack of fit is reduced. But there are still problems in the model. It is inadequate to fit the shape of our data.

Since we did so well with the quadratic, try out a third order model. Yep, it does better still. Wow. Keep pushing things, until we have a 6th order model, or 10th order. A polynomial model is just like a Taylor series approximation, so more terms must give a better result. Eventually though, it all just turns to numerical crap before our eyes. So what happened?

Here you have 3 variables. A first order model has a constant term, as well as x, y, and z terms. So 4 terms in that model. You have 74 data points, and surely had no problem estimating the corresponding coefficients.

A quadratic model adds in 6 more terms, x^2, y^2, z^2, but also the cross product terms x*y, x*z, and y*z, so a total of 10 coefficients to estimate.

The curse of our fallacy is that as we move up in order, the number of terms grows faster and faster. At 3rd order, we have 20 terms to estimate, 35 at 4th order, 56 at 5th order, and 84 at 6th order. A 10th order model would have 286 terms.

How about your data? You have only 74 data points. (In my early days when I did some statistical consulting, my clients would always ask how much data they needed. My answer was almost always "More than you have.")

Every distinct data point is at most one additional piece of information. (Replicate points do not count here. All they do is reduce the uncertainty in your estimates. So replicates are good, but not quite as valuable as distinct points in some respects.)

Worse, new data points do not always provide new information. For example, two points determine a line. A third point on that line adds no new information towards the slope or intercept of that line. Again, like a replicate, all it does is reduce the uncertainty in your estimate of the coefficients.

So what happened when you tried to estimate a 6th order model with only 74 data points? You had too little information in your data to support such a large model. The message that comes back is from MATLAB, telling you that your matrix is numerically singular (rank deficient.)

Having said all this, things can get worse yet, making the problem yet more difficult. Some sets of data are worse than others. For example, if you add 1e6 to all of your independent variables, a polynomial of some given order will still be estimable mathematically (in theory), yet it may result in a matrix that is numerically singular. This is because the linear algebra is written using double precision floating point arithmetic. The numbers may simply lie over too wide of a dynamic range for the algebra to succeed. (What happens when you raise numbers that are on the order of 1e6 to successively higher powers?)

So when you have more coefficients to estimate than data points, you will always have a rank deficient matrix, but depending on the data itself, you may still have rank deficiency problems.

Some polynomial models will become more tractable to estimate if you center and scale the variables. Thus make each variable so it has a mean of zero, and lie roughly in the interval [-1,1]. See that raising numbers that are on the order of +/-1 to higher powers does NOT create a dynamic range issue. This will force you to transform your problem, but it often allows you to gain viable estimates for an otherwise intractable problem. (The use of orthogonal polynomials is another approach that often helps here.)

I hope this helps. Really, this is a subject that I've seen taught as a semester long course, so I have only touched the surface.

John

John,

Great code as usual!

It seems to fail when given only one data point.

n = 1;
x = rand(n,2);
y = exp(sum(x,2)) + randn(n,1)/100;
p = polyfitn(x,y,3);

While failing might be desirable in some cases, it should really throw an error instead. And your code does supply an n-param output for n-1 inputs, so it would not be inconsistent to handle this case.

I think it is in the logic after line 120 - checking the size of the inputs and transposing accordingly.

Thanks!

John - There are two possible questions here. First, if you have a problem with what is known as multiple right hand sides, then yes it is very vectorizable. By this I mean a problem with the same sets of independent variables, but different sets of dependent variables. Then backslash is the way to go, allowing you to factorize the equations once for a huge time savings. You would need to build the design matrix of course, but this is not a difficult thing to do.

If however, you have multiple problems, each with different sets of independent variables, then vectorization would require you to build a sparse block diagonal design matrix. As such, you can gain a bit in this effort, but not as significant of a gain as the first case. You need to build the design matrix as a sparse matrix. The trick here is to supply the blocks to blkdiag as individually sparse matrices. Thus, create the individual (sparse) design matrices in a cell array, then pass it to blkdiag as a comma separated list. Can you do this efficiently enough to make the backslash solve more efficient than a loop around the calls? Yes, but as I said, it will take some effort.

You can probably glean from my response that it will not be a trivial action to vectorize the result from this code, but that you could do things more efficiently with some time invested.

This should do it...

Look for these two lines below:

[Q,R,E] = qr(M,0);

polymodel.Coefficients(E) = R\(Q'*depvar);

Assume a column vector of weights as W, one for each data point. Change the above lines to read:

[Q,R,E] = qr(bsxfun(@times,W(:),M),0);

polymodel.Coefficients(E) = R\(Q'*(W(:).*depvar));

This should give you weighted estimates of the parameters. I think the statistics on the parameters will work too, as well as R^2, etc. Of course, you will need to pass in the weight vector too as an argument.

The above change assumes you are using a version that has bsxfun available. Its been around for a while so that should not be a problem. Even if not, a call to repmat will replace it.

A linear regression (which is polyfitn) uses backslash at its core. For example, suppose one wished to use polyfit to estimate a quadratic model in one variable. The model would be

y = a1*x^2 + a2*x + a3 (+noise)

polyfit is able to estimate the three unknown coefficients [a1,a2,a3] by solving a linear (but overdetermined) system of equations in a way that minimizes the sum of squares of the residuals. Again, backslash is the engine underneath here.

Now, suppose one wished to estimate a more complex model, with two independent variables, z(x,y). I'll write it as a full two variable quadratic model, so there will be six unknown coefficients.

z = a1*x^2 + a2*x*y + a3*y^2 + a4*x + a5*y + a6 + noise

As long as you have sufficient data to estimate the six unknown coefficients, the problem is essentially the same. See that they enter the problem linearly in the same way as for the one variable problem we looked at before. Backslash again can solve the overdetermined problem to estimate the coefficients.

In fact, both polyfit and polyfitn use a QR factorization of the design matrix to achieve their goals. This is because they are both designed to also generate covariance matrix estimates on the parameters. A QR factorization is the easiest way to serve both purposes, but its all just linear algebra.