Solving linear system with rank-deficient matrix in MATLAB

Question

I want to perform quartic regression on following data:

x =
     375
    1029
    2187
    3993
    6591
   10125
   14739

y=
   0.5608
   0.8522
   0.9075
   0.9994
   0.8668
   0.7162
   0.9143

To do that I've implemented following function:

function[str] = quarticRegression(x, y)

    A = [x.^4 x.^3 x.^2 x ones(length(x), 1)];
    res = A\y;

    str = sprintf('%fx^4 + (%f)x^3 + (%f)x^2 + (%f)x + %f',...
                  res(1), res(2), res(3), res(4), res(5));
end

Unfortunately, I'm getting following error:

Warning: Rank deficient, rank = 4, tol =  7.520684e+01.

When I try to change x vector for smaller numbers, error disappears.

What is the correct way to get the result?

I guess I should somehow normalize the x vector, but I don't know how to do that and not change the system solution.

Answer 1

The condition number of $A$ (the ratio of the largest singular value of A to the smallest) is very big $\left(\approx 10^{16}\right)$ so that Matlab iterprets $A$ as a singular matrix.

Even the Matlab function polyfit (that do the job you want) shows that the problem is badly conditioned. Indeed, notice that the graph lying the points $(x,y)$ is almost an horizontal line.

With your data, we have for a quadratic fit

>> p = polyfit(x,y,2)

p =

-0.0000    0.0000    0.7501

which is approximately a constant value.

Answer 2

Matlab solution: Load your data as x and y vector (check them in the workspace), I had NaN in the first values.

polyfit(x,y,2) % last argument is for the order of the polynomial

Excel solution (different solutions as Excel is using linear regression and I think Matlab is using nonlinear regression): This is not an answer for Matlab, bu if you just want to find a simple polynomial regression. You can always use Excel for this. Just plot the data in an x,y-plot and left-click >> add trend line >> select polynomial with degree 2 and check the box for view formula and check the box for display coefficient of determination.

To be quite honest I think a cubic polynomial would suit much better (loot at the plot).

Answer 3

To scale x, just multiply it by a suitable number. Thus:

z = 0.01 * x;

With $z$ instead of $x$, it seems to work. Note that $y = a_0 + a_1 z + a_2 z^2 + a_3 z^3 + a_4 z^4$ corresponds to $y = a_0 + 10^{-2} a_1 x + 10^{-4} a_2 x^2 + 10^{-6} a_3 x^3 + 10^{-8} a_4 x^4$.

A = [z.^4 z.^3 z.^2 z ones(length(z), 1)];
a = A\y;
ax = a .* [10^(-8),10^(-6),10^(-4),10^(-2),1]';
X = 300:100:16000;
R = ax(1) * X .^ 4 + ax(2) * X .^ 3 + ax(3) * X .^2 + ax(4) * X + ax(5);
plot(X,R)
hold
plot(x,y,'Color','red')