Re: Polyfit with scaling

[Ben Abbott]

I invested some time looking at the stability of polyfit under different conditions. The four cases I considered were

(1) Using the present algorithm (which uses the backslash).

(2) Using the present algorithm, with normalized inputs.

(3) Using QR for the solution.

(4) Using QR for the solution, with normalized inputs.

The normalization used is

x = (x-mean(x)) / abs(std(x));

y = (y-mean(y)) / abs(std(y));

The short script I used permits the order of the polynomial and the length of the input vectors to be varied independently.

xoffset = 1000;

pshift = ones (1, N+1);

p = polyshift (pshift, -xoffset);

x = xoffset + linspace(-5,5,M);

y = polyval (pshift, x-xoffset);

[pfit, s] = polyfit (x, y, N);

yfit = polyval (pfit, x);

I ran several cases and compared the norm of the residues, as well as the norm of the difference of the original and fitted polynomial coefficients.

                            norm of (yfit-y)

             Present Algorithm         Solution using QR

   N   M     unnorm        norm        unnorm        norm

   1   2   1.1369e-13   0.0000e+00   1.1369e-13   1.2561e-15

   1   3   1.1369e-13   1.2561e-15   1.1369e-13   1.2610e-15

   2   3   2.6031e-10   0.0000e+00   1.3824e-09   8.7504e-15

   2   4   6.6134e-08   1.1409e-14   3.4041e-10   1.1580e-14

   2   5   6.9249e-08   8.2486e-15   5.4604e-10   1.1572e-14

   3   4   4.9642e+01   1.8310e-15   1.3767e-06   6.3815e-14

   3   5   5.9235e+01   1.4262e-13   1.0926e-06   6.4272e-14

   3   6   6.4335e+01   1.8519e-14   2.6008e-06   3.5718e-14

   3   7   6.7974e+01   1.0843e-13   1.3334e-06   6.1238e-14

   4   5   1.2522e+02   3.5527e-15   2.0280e-03   2.2775e-13

   4   6   1.5735e+02   2.2581e-12   1.9224e-03   3.3016e-13

   4   7   1.7629e+02   9.1410e-13   1.5611e-03   1.6262e-13

   4   8   1.8952e+02   3.8459e-13   3.6440e-03   4.6810e-13

   4   9   1.9985e+02   1.7892e-13   1.6603e-03   1.7505e-13

                            norm of (pfit-p)

             Present Algorithm         Solution using QR

   N   M     unnorm        norm        unnorm        norm

   1   2   6.7075e-12   0.0000e+00   2.8422e-12   0.0000e+00

   1   3   1.6939e-11   0.0000e+00   1.7167e-11   0.0000e+00

   2   3   2.8582e-06   2.3283e-10   1.2274e-05   4.6566e-10

   2   4   3.3198e-06   4.6566e-10   3.3194e-06   4.6566e-10

   2   5   4.9461e-06   4.6566e-10   4.9466e-06   4.6566e-10

   3   4   9.9901e+08   4.7684e-07   2.0040e+00   1.1921e-07

   3   5   9.9901e+08   4.7684e-07   2.9059e+00   2.3842e-07

   3   6   9.9901e+08   2.3842e-07   1.2420e+01   4.7684e-07

   3   7   9.9901e+08   7.1526e-07   7.3990e+00   4.7684e-07

   4   5   9.9901e+11   3.6621e-04   2.7385e+06   3.6621e-04

   4   6   9.9901e+11   5.9815e-03   9.0039e+05   8.5450e-04

   4   7   9.9901e+11   3.0518e-03   1.6692e+06   8.5450e-04

   4   8   9.9901e+11   2.4414e-04   4.4956e+05   1.2207e-03

   4   9   9.9901e+11   1.2207e-03   3.2162e+06   4.8828e-04

With no normalization, the QR algorithm performs better than the current implementation, Although, when using normalized variables that benefit is nearly absent.

I tried other polynomials as well. All I tried gave similar comparative results (including poles of multiplicity).

Give these results and the need to remain compatible with past expectations, and with Matlab, I'm inclined to make both QR and normalization optional.

In the event, the optional approach is to be used, I'm prepared to incorporate it in other core functions, such as residue.m.

Thoughts?

Ben