Why is the POLYFIT function in MATLAB unable to find a fit over my data values?
39 views (last 30 days)
MathWorks Support Team on 24 Jan 2013
I am using the POLYFIT function to fit a second order polynomial over my data values as follows.
However, I receive the following warning message
ERROR: Warning: Polynomial is badly conditioned. Add points with
distinct X values, reduce the degree of the polynomial, or try
centering and scaling as described in HELP POLYFIT.
MathWorks Support Team on 30 Nov 2022
Edited: MathWorks Support Team on 5 Dec 2022
The example fitting will throw a warning in R2009a but will not in later releases due to changes in how POLYFIT throws warnings. The numerical properties of the problem have not changed, it is still ill-conditioned.
The warning message occurs is due to the scaling of the data. It produces an ill-conditioned least squares problem which POLYFIT flags. One way to address this is to scale the data, such that each "X" in the original input is scaled into "X_new" as follows:
X_new = (X – mean(X))/std(X)
This scales the input values so that x' is contained in [-1 1]. The third calling syntax of POLYFIT does this scaling before fitting the polynomial and returns a third output:
mu = [mean(X) std(X)]
This will return a different set of coefficients than as if the data were to be fitted with no scaling.
you will receive:
p*(X_new) = A_new*(X_new)^2 +B_new*(X_new) + C_new
%For notational purposes, the unscaled polynomial is
p(X) = A*(X)^2 +B*(X) + C
These new coefficients are correct in the sense that p(X_new) will produce the same Y value as its corresponding X and the fitting can be used to estimate any query point if it is scaled properly. These scaled coefficients can be transformed back to their non-scaled form with some caveats which I will address later. Since these two equations must be equal for corresponding X and X_new values, there is a relationship between the coefficients for a scaled (A_new,B_new,C_new) and un-scaled (A,B,C).
A_new*(X_new)^2 +B_new *(X_new) + C_new = A*(X)^2 + B*(X) + C
By writing X_new in terms of X as above, this relationship can be manipulated to generate closed form expressions for A,B, and C:
[p,S,mu] = polyfit(X,y,2)
A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
A = A_new/(mu2)^2
B = (B_new*mu2)/mu2 – (2*A_new*mu1)/(mu2)^2
However, this approach can be generalized for higher order fittings using the binomial theorem. (https://en.wikipedia.org/wiki/Binomial_theorem ). Consider the closed form equations above, but in matrix form, with A, b, and C being the sums of the columns. \n\n
M = [nchoosek(2,0) nchoosek(2,1)*(-mu1) nchoosek(2,2)*(-mu1)^2;
0, nchoosek(1,0), nchoosek(1,1)*(-mu1);
0, 0, choosek(0,0)];
T = ( p.' ./ [mu2^2 mu2 1].' ) .* M
This should hopefully make it a bit clearer that the rows of the matrix M are the coefficients of the expansion of:
(X - mu1)^k
The general form of which is given by the Binomial Theorem.
Lastly, a couple of caveats; This is not a recommended workflow to manually convert from scaled to non-scaled coefficients. Manually transforming the data is largely is redundant as you can produce the un-scaled coefficients by simply calling POLYFIT with the 2 output syntax in addition to the 3 output syntax. However, please note that these results can differ slightly from the manually transformed coefficients. See the well-conditioned example from the documentation below using the generalized matrix form of the transformation to demonstrate.
>> load census
>> p = polyfit(cdate,pop,2);
>> [phat,S,mu] = polyfit(cdate,pop,2);
>> T = ( phat.' ./ [mu(2)^2 mu(2) 1].' ) .* [nchoosek(2,0) nchoosek(2,1)*(-mu(1)) nchoosek(2,2)*(-mu(1))^2; 0 nchoosek(1,0) nchoosek(1,1)*(-mu(1)); 0 0 nchoosek(0,0)]
0.0000 -0.0025 2.3366
0 0.0001 -0.2298
0 0 0.0062
>> pT = sum(T) % sum over cols
0.0000 -0.0024 2.1130
0.0000 -0.0024 2.1130
>> p - pT
-0.0000 0.0008 -0.7785
Note that even in this relatively nice example there are some numerical differences since the coefficients have been calculated differently. Now let's turn our attention back to the original example
>> x = -100000:10000:200000;
f = @(x) 1000*x.^2+ 3*x + 0.023;
y = f(x);
>> p = polyfit(x,y,2);
[phat,S,mu] = polyfit(x,y,2);
T = ( phat.' ./ [mu(2)^2 mu(2) 1].' ) .* [nchoosek(2,0) nchoosek(2,1)*(-mu(1)) nchoosek(2,2)*(-mu(1))^2; 0 nchoosek(1,0) nchoosek(1,1)*(-mu(1)); 0 0 nchoosek(0,0)]
0.0000 -0.0001 2.5000
0 0.0001 -5.0000
0 0 2.5000
>> pT = sum(T)
1.0000 0.0030 0.0000
1.0000 0.0030 0.0000
>> p - pT
0.0000 -0.0000 0.0025
In this case you will see that the numerical difference has increased by nearly 10^7. For similar situations, it is likely that the newly calculated non-scaled coefficients will not be an accurate fit for your data, and all you have done is essentially circumvented the warnings for an ill-conditioned problem.
More Answers (2)
Shaik Pasha on 27 Oct 2014
Hi, the provided information for 2nd order polynomial is very much useful. But my query is what would be the transformation equations for the 3rd and 4th order polynomial equations, wherein i had to extract [A,B,C,D] and [A,B,C,D,E], respectively.
Hope for the earliest reply.
Christoph on 31 Oct 2014
Edited: MathWorks Support Team on 6 Apr 2018
@Shaik: Basically, you have to exploit the binomial theorem http://en.wikipedia.org/wiki/Binomial_theorem to get a formula for arbitrary degrees of your polynomial. A working MATLAB solution can be found here: https://groups.google.com/forum/#!topic/comp.soft-sys.matlab/D4JPQhkylkE - the polyfit_convert function.