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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20080203/2e61385b/attachment.html