Polyfit with scaling

[Thomas Weber-8]

I'd like to propose the attached patch for polyfit.m
It introduces an optional fourth bool parameter 'scale'. If not set, or
if it's set to 'false', nothing changes.

If it's set to true, the x values are scaled by their maximal absolute
value prior to calculating the polynomial coefficients. This has a
stabilizing effect on the calculation.

A test case showing the difference is included.

There are three FIXME's in the code. Two are about the usage of repmat()
and reshape, respectively (speed reasons). The third one deals with the
second parameter from chol(), which is undocumented. So I don't know
what the code should do in case this return value is non-zero (the
current code ignores this).

        Thomas

---

diff -r b280628de44e scripts/polynomial/polyfit.m
--- a/scripts/polynomial/polyfit.m Mon Jan 28 21:18:03 2008 +0100
+++ b/scripts/polynomial/polyfit.m Wed Jan 30 18:13:39 2008 +0100
@@ -51,17 +51,31 @@
 ## @item yf
 ## The values of the polynomial for each value of @var{x}.
 ## @end table
+##
+## If a fourth bool parameter @var{scale} is given and set to
+## @code{true}, the values of @var{x} are scaled by their maximum absolute
+## value before doing the actual calculation. This helps in numercial difficult
+## cases. The returned coefficients @var{p} are valid for the unscaled
+## values of @var{x}, while the values in the structure @var{S} are
+## valid for the scaled polynomial used for the actual calculation.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@...>
 ## Created: 13 December 1994
 ## Adapted-By: jwe
 
-function [p, s, mu] = polyfit (x, y, n)
+function [p, s, mu] = polyfit (x, y, n, scale)
 
+  if (nargin < 3 || nargin > 4)
+    print_usage ();
+  endif
 
-  if (nargin != 3)
-    print_usage ();
+  if (nargin == 3)
+    scale = false;
+  else
+    if ( ! isbool(scale) || ! isscalar(scale))
+      error ("polyfit: scale must be a bool scalar");
+    endif
   endif
 
   if (! (isvector (x) && isvector (y) && size_equal (x, y)))
@@ -77,8 +91,16 @@ function [p, s, mu] = polyfit (x, y, n)
   l = length (x);
   x = reshape (x, l, 1);
   y = reshape (y, l, 1);
+  ## FIXME: is reshape() faster than (:) ?
+
+  ## scale the x values
+  if (scale)
+    xmax = max(abs(x));
+    x = x/xmax;
+  endif
 
   X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0));
+  # FIXME: maybe use repmat() for the ones() part, but this here seems faster
 
   p = X \ y;
 
@@ -93,16 +115,22 @@ function [p, s, mu] = polyfit (x, y, n)
     endif
 
     [s.R, dummy] = chol (X'*X);
+    ## FIXME: what if dummy is !0 ?
     s.X = X;
     s.df = l - n - 1;
     s.normr = norm (yf - y);
 
   endif
 
-  ## Return value should be a row vector.
+  ## Scale the polynomial back and return a row vector
+  if (scale)
+    power = length(p)-1:-1:0;
+    power = (xmax .^ power)(:);
+    p = (p ./ power);
+  endif
 
-  p = p.';
-
+  ## Return a row vector
+  p = p(:)';
 endfunction
 
 %!test
@@ -123,3 +151,20 @@ endfunction
 %! x = [-2, -1, 0, 1, 2];
 %! fail("polyfit (x, x.^2+x+1, [])");
 
+%% test difficult case where scaling is really needed
+%% checks also usage of second output argument
+%!shared s1, s2
+%!  x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, \
+%!        -1186.8, -1185.6, -1184.4, -1183.2, -1182 ];
+%!  y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966,    \
+%!        315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201,    \
+%!        315600.7143, 315602.9508, 315605.1765 ];
+%!  [p1,s1] = polyfit(x,y,10);
+%!  [p2,s2] = polyfit(x,y,10, true);
+%!
+%!  # fails without scaling (s1.norm ~6.33)
+%!error assert(s1.normr, 0.08,  0.01);
+%!
+%!  # but works with scaling
+%!assert(s2.normr, 0.08,  0.01);
+

---

[Dmitri A. Sergatskov]

On Jan 30, 2008 11:19 AM, Thomas Weber <thomas.weber.mail@...> wrote:
> I'd like to propose the attached patch for polyfit.m
> It introduces an optional fourth bool parameter 'scale'. If not set, or
> if it's set to 'false', nothing changes.
>
> If it's set to true, the x values are scaled by their maximal absolute
> value prior to calculating the polynomial coefficients. This has a
> stabilizing effect on the calculation.
>

I do not understand why this change had not been made a long time ago
-- we had few reports of faulty fits. I think the better
scaling is:

x = (x-mean(x))/std(x).

wpolyfit from octave-forge does the scaling like that.


>         Thomas
>


Dmitri.
--

[Thomas Weber-8]

Am Mittwoch, den 30.01.2008, 17:48 -0600 schrieb Dmitri A. Sergatskov:

> On Jan 30, 2008 11:19 AM, Thomas Weber <thomas.weber.mail@...> wrote:
> > I'd like to propose the attached patch for polyfit.m
> > It introduces an optional fourth bool parameter 'scale'. If not set, or
> > if it's set to 'false', nothing changes.
> >
> > If it's set to true, the x values are scaled by their maximal absolute
> > value prior to calculating the polynomial coefficients. This has a
> > stabilizing effect on the calculation.
> >
>
> I do not understand why this change had not been made a long time ago
> -- we had few reports of faulty fits. I think the better
> scaling is:
>
> x = (x-mean(x))/std(x).
>
> wpolyfit from octave-forge does the scaling like that.

Tarball with a test directory attached. Output of the script in it:

norm difference with normal polyfit: 6.330024
norm difference with wpolyfit: 0.465359
norm difference with new polyfit: 0.082466

There are warnings from dgelsd, though. The data is taken from the unit
test, which in turn was taken from a bug report in Debian.

I don't use polynomial approximations often, so it would be nice if a
few people tested it with their data.

        Thomas

[Ben Abbott]

On Jan 31, 2008, at 3:40 AM, Thomas Weber wrote:

Am Mittwoch, den 30.01.2008, 17:48 -0600 schrieb Dmitri A. Sergatskov:

On Jan 30, 2008 11:19 AM, Thomas Weber <thomas.weber.mail@...> wrote:

I'd like to propose the attached patch for polyfit.m

It introduces an optional fourth bool parameter 'scale'. If not set, or

if it's set to 'false', nothing changes.

If it's set to true, the x values are scaled by their maximal absolute

value prior to calculating the polynomial coefficients. This has a

stabilizing effect on the calculation.

I do not understand why this change had not been made a long time ago

-- we had few reports of faulty fits. I think the better

scaling is:

x = (x-mean(x))/std(x).

wpolyfit from octave-forge does the scaling like that.

Tarball with a test directory attached. Output of the script in it:

norm difference with normal polyfit: 6.330024
norm difference with wpolyfit: 0.465359
norm difference with new polyfit: 0.082466

There are warnings from dgelsd, though. The data is taken from the unit
test, which in turn was taken from a bug report in Debian.

I don't use polynomial approximations often, so it would be nice if a
few people tested it with their data.

Thomas

<polyfits.tar.gz>

The test script uses a set of data for which there is a very small variation in both x and y. As a result, if I make the non-sensical change below;

  ## scale the x values
  if (scale)
 #  xmax = max(abs(x));
    xmax = min(abs(x));
    x = x/xmax;
  endif

I still get impressive results

warning: dgelsd: rank deficient 13x11 matrix, rank = 3

warning: dgelsd: rank deficient 13x11 matrix, rank = 7

norm difference with normal polyfit: 6.330024

norm difference with wpolyfit: 0.150155

norm difference with new polyfit: 0.082631

It doesn't appear that the example has much comparative value. 

I suspect that using max(abs(x)) to normalized the x-data will also result in difficulties when that value is substantially large than the other x-data values.

I suggest some tests be designed to break both wpolyfit as well as your new version, and see how the three fair comparatively.

I'll have some time over the next few days, and will take a shot at writing a few tests myself.

Ben

[Thomas Weber-8]

On 31/01/08 07:47 -0500, Ben Abbott wrote:
> It doesn't appear that the example has much comparative value.

Honestly, I don't care what we change. However, what shouldn't happen is
that things just stay the same. If we have a better implementation in
wpolyfit, we should take that for polyfit. Having a better algorithm
lying around and not using is just a waste (actually, it's even worse:
we not even have a better algorithm, but it's also already implemented).

        Thomas

[Ben Abbott]

On Jan 31, 2008, at 3:08 PM, Thomas Weber wrote:

On 31/01/08 07:47 -0500, Ben Abbott wrote:

It doesn't appear that the example has much comparative value.

Honestly, I don't care what we change. However, what shouldn't happen is
that things just stay the same. If we have a better implementation in
wpolyfit, we should take that for polyfit. Having a better algorithm
lying around and not using is just a waste (actually, it's even worse:
we not even have a better algorithm, but it's also already implemented).

Thomas

? hmmm ... I agree we  should use the better algorithm.

In any event, I've modified polyfit.m to optionally use the method employed by wpolyfit. To do this; I wrote two new scripts. One to offset a polynomial's dependent variable (polyshift), and one to scale a polynomial's dependent variable (polyscale).

Below I compared the present version of polyfit to the one attached to this email.

octave:1> x = 1:4;

octave:2> y = polyval ([1, 1, 1], x);

octave:3> polyfit (x, y, 2)

ans =   1.00000   1.00000   1.00000

octave:4> polyfit (x, y, 2, 1)

ans =   1.00000   1.00000   1.00000

octave:5> polyfit (x, y, 3)

ans =   0   1   1   1

octave:6> polyfit (x, y, 3, 1)

ans =   5.2633e-16   1.0000e+00   1.0000e+00   1.0000e+00

octave:7> polyfit (x, y, 4)

ans =  -0.0020445   0.0204453   0.9284416   1.1022263   0.9509314

octave:8> p = polyfit (x, y, 4, 1)

p =    -5.0886    50.8856  -177.0996   255.4281  -121.1255

octave:9> polyval (p, x)

ans =    3.0000    7.0000   13.0000   21.0000

octave:10> polyval ([1, 1, 1], x)

ans =    3    7   13   21

Comparing to Matlab ...

>> x = 1:4;

>> y = polyval([1, 1, 1], x);

>> polyfit(x, y, 2)

ans =    1.0000    1.0000    1.0000

>> polyfit(x, y, 3)

ans =    0.0000    1.0000    1.0000    1.0000

>> polyfit(x, y, 4)

Warning: Polynomial is not unique; degree >= number of data points.

> In polyfit at 72

ans =    0.0200   -0.2000    1.7000    0    1.4800

The offset and scaling appears to significantly modify the result when the order of the fit exceeds to order of the original polynomial.

I am not an expert in such matters. How might we quantitatively determine which algorithm is best?

I've attached the modified polyfit.m as well as the two new scripts it relies upon.

Perhaps it is worth the effort for others to compare their results with mine, and to compare Matlab's results on a computer other than an Intel/ppc based Mac (each gave the same result).

Ben

[Ben Abbott]

On Jan 31, 2008, at 9:26 PM, Ben Abbott wrote:

On Jan 31, 2008, at 3:08 PM, Thomas Weber wrote:

On 31/01/08 07:47 -0500, Ben Abbott wrote:

It doesn't appear that the example has much comparative value.

Honestly, I don't care what we change. However, what shouldn't happen is
that things just stay the same. If we have a better implementation in
wpolyfit, we should take that for polyfit. Having a better algorithm
lying around and not using is just a waste (actually, it's even worse:
we not even have a better algorithm, but it's also already implemented).

Thomas

? hmmm ... I agree we  should use the better algorithm.

In any event, I've modified polyfit.m to optionally use the method employed by wpolyfit. To do this; I wrote two new scripts. One to offset a polynomial's dependent variable (polyshift), and one to scale a polynomial's dependent variable (polyscale).

Below I compared the present version of polyfit to the one attached to this email.

octave:1> x = 1:4;

octave:2> y = polyval ([1, 1, 1], x);

octave:3> polyfit (x, y, 2)

ans =   1.00000   1.00000   1.00000

octave:4> polyfit (x, y, 2, 1)

ans =   1.00000   1.00000   1.00000

octave:5> polyfit (x, y, 3)

ans =   0   1   1   1

octave:6> polyfit (x, y, 3, 1)

ans =   5.2633e-16   1.0000e+00   1.0000e+00   1.0000e+00

octave:7> polyfit (x, y, 4)

ans =  -0.0020445   0.0204453   0.9284416   1.1022263   0.9509314

octave:8> p = polyfit (x, y, 4, 1)

p =    -5.0886    50.8856  -177.0996   255.4281  -121.1255

octave:9> polyval (p, x)

ans =    3.0000    7.0000   13.0000   21.0000

octave:10> polyval ([1, 1, 1], x)

ans =    3    7   13   21

Comparing to Matlab ...

>> x = 1:4;

>> y = polyval([1, 1, 1], x);

>> polyfit(x, y, 2)

ans =    1.0000    1.0000    1.0000

>> polyfit(x, y, 3)

ans =    0.0000    1.0000    1.0000    1.0000

>> polyfit(x, y, 4)

Warning: Polynomial is not unique; degree >= number of data points.

> In polyfit at 72

ans =    0.0200   -0.2000    1.7000    0    1.4800

The offset and scaling appears to significantly modify the result when the order of the fit exceeds to order of the original polynomial.

I am not an expert in such matters. How might we quantitatively determine which algorithm is best?

I've attached the modified polyfit.m as well as the two new scripts it relies upon.

Perhaps it is worth the effort for others to compare their results with mine, and to compare Matlab's results on a computer other than an Intel/ppc based Mac (each gave the same result).

Ben

I should mention that nothing has been done to ensure that the second output respects the normalization of the dependent variable. 

I'll look into what needs to be done there, if it is possible ... or more to the point, if I'm up to the task ;-)

If anyone has experience or background in how the job might be done, please feel free to contribute or to point me in the right direction!

Ben

[Rolf Fabian]

Ben Abbott wrote:

On Jan 31, 2008, at 3:08 PM, Thomas Weber wrote:

> On 31/01/08 07:47 -0500, Ben Abbott wrote:
>> It doesn't appear that the example has much comparative value.
>
> Honestly, I don't care what we change. However, what shouldn't  
> happen is
> that things just stay the same. If we have a better implementation in
> wpolyfit, we should take that for polyfit. Having a better algorithm
> lying around and not using is just a waste (actually, it's even worse:
> we not even have a better algorithm, but it's also already  
> implemented).
>
> Thomas

? hmmm ... I agree we  should use the better algorithm.

In any event, I've modified polyfit.m to optionally use the method  
employed by wpolyfit. To do this; I wrote two new scripts. One to  
offset a polynomial's dependent variable (polyshift), and one to scale  
a polynomial's dependent variable (polyscale).

Below I compared the present version of polyfit to the one attached to  
this email.

octave:1> x = 1:4;
octave:2> y = polyval ([1, 1, 1], x);
octave:3> polyfit (x, y, 2)
ans =   1.00000   1.00000   1.00000
octave:4> polyfit (x, y, 2, 1)
ans =   1.00000   1.00000   1.00000
octave:5> polyfit (x, y, 3)
ans =   0   1   1   1
octave:6> polyfit (x, y, 3, 1)
ans =   5.2633e-16   1.0000e+00   1.0000e+00   1.0000e+00
octave:7> polyfit (x, y, 4)
ans =  -0.0020445   0.0204453   0.9284416   1.1022263   0.9509314
octave:8> p = polyfit (x, y, 4, 1)
p =    -5.0886    50.8856  -177.0996   255.4281  -121.1255
octave:9> polyval (p, x)
ans =    3.0000    7.0000   13.0000   21.0000
octave:10> polyval ([1, 1, 1], x)
ans =    3    7   13   21

Comparing to Matlab ...

 >> x = 1:4;
 >> y = polyval([1, 1, 1], x);
 >> polyfit(x, y, 2)
ans =    1.0000    1.0000    1.0000
 >> polyfit(x, y, 3)
ans =    0.0000    1.0000    1.0000    1.0000
 >> polyfit(x, y, 4)
Warning: Polynomial is not unique; degree >= number of data points.
 > In polyfit at 72
ans =    0.0200   -0.2000    1.7000    0    1.4800

The offset and scaling appears to significantly modify the result when  
the order of the fit exceeds to order of the original polynomial.

I am not an expert in such matters. How might we quantitatively  
determine which algorithm is best?

I've attached the modified polyfit.m as well as the two new scripts it  
relies upon.

Perhaps it is worth the effort for others to compare their results  
with mine, and to compare Matlab's results on a computer other than an  
Intel/ppc based Mac (each gave the same result).

Ben

Dear Ben

Without testing it, I looked shortly at the code of
your polyfit.m and found  close to the end:

  ## Return a row vector
  p = p(:)';

This is a bug because it converts complex column vectors
to complex conjugated polynomials ( row-vectors)

Rolf

[Ben Abbott]

On Feb 1, 2008, at 9:15 AM, Ben Abbott wrote:

>
> On Jan 31, 2008, at 9:26 PM, Ben Abbott wrote:
>
> I should mention that nothing has been done to ensure that the  
> second output respects the normalization of the dependent variable.
>
> I'll look into what needs to be done there, if it is possible ... or  
> more to the point, if I'm up to the task ;-)
>
> If anyone has experience or background in how the job might be done,  
> please feel free to contribute or to point me in the right direction!
>
> Ben
>

I've added some code to mock the 2nd output, and have fixed a bug.

In this version, the Cholesky factor and the Vandermonde matrix which  
are returned as fields on the 2nd output are not used in the  
solution.  The matrices that are used are representative of the scaled  
version of x.

While this is not a pretty solution, I don't know what else might be  
done. Any thoughts?

Ben

[Thomas Weber-8]

On 01/02/08 07:26 -0800, Rolf Fabian wrote:
> Without testing it, I looked shortly at the code of
> your polyfit.m and found  close to the end:
>
>   ## Return a row vector
>   p = p(:)';
>
> This is a bug because it converts complex column vectors
> to complex conjugated polynomials ( row-vectors)

I think Ben just took this line from me. Your observation regarding the
complex values is obviously correct. However, is polyfit supposed to
work with comples input values? Otherwise, I don't see how complex
values can show up.

        Thomas

[Ben Abbott]

On Feb 1, 2008, at 3:57 PM, Thomas Weber wrote:

> On 01/02/08 07:26 -0800, Rolf Fabian wrote:
>> Without testing it, I looked shortly at the code of
>> your polyfit.m and found  close to the end:
>>
>>  ## Return a row vector
>>  p = p(:)';
>>
>> This is a bug because it converts complex column vectors
>> to complex conjugated polynomials ( row-vectors)
>
> I think Ben just took this line from me. Your observation regarding  
> the
> complex values is obviously correct. However, is polyfit supposed to
> work with comples input values? Otherwise, I don't see how complex
> values can show up.
>
> Thomas

I don't believe Octave has any functions that are restricted to real  
values. So ... yes, polyfit accepts and returns complex values.

Ben

[Ben Abbott]

Using the same approach as applied in wpolyfit, it is possible to  
improve the numerical stability/accuracy of polyfit.

However, it is not clear how the second output should be handled.

[P, S] = polyfit (X, Y, N)

S.R: The Cholesky factor of the Vandermonde matrix used to compute the  
polynomial coefficients.
S.X: The Vandermonde matrix used to compute the polynomial coefficients.
S.df: The degrees of freedom.
S.normr: The norm of the residuals.
S.yf: The values of the polynomial for each value of X.

Essentially, when there is no normalization of X the outputs are

         L = numel (X);
        V = (X * ones (1, N+1)) .^ (ones (L, 1) * (N : -1 : 0));
        P = (V \ Y).';
        S.yf = polyval (P, X);
        S.df = numel (X) - N;
        S.X  = V;
        s.R  = chol (V'*V);

When X is normalized

        Xn = (X - mean (X)) / std (X - mean (X));
        V = (Xn * ones (1, n+1)) .^ (ones (L, 1) * (n : -1 : 0));
        P = (V \ Y).';
         P = polyscale (P, std (X - mean (X)));
        P = polyshift (P, -mean (X));

I'm not sure what should be done with the second output when  
normalization is in place. I think the proper treatment depends upon  
what the results are used for.

Does anyone have experience/insight into how/where S.S, and S.R are  
used?

Ben

[Dmitri A. Sergatskov]

Though I was partially responsible for all this discussions,
I am now rethinking this approach.
polyfit should be as stable numerically as we can make it,
but I do not think it should try to make a "better"
fit than the one it was asked for.

Polyfit does fit to a*x^n+b*x^(n-1)... polynomial.
If people want it to fit to a*(x-x0)^n+b*(x-x0)^(n-1)
they should pre-condition data accordingly prior to
passing it to polyfit.

So to make long story short -- I think the original
Thomas' suggestion (just to normalize the data to max(x))
was the good one.

On related note -- it appears that wpolifit does a
better job than polyfit on un-conditioned data as well.
It appears that the difference is in polyfit
uses straight "\",  while the wpolyfit does QR decomposition.
Perhaps it is worth to port this to polyfit as well.
Frankly I do not see why would not we copy wpolyfit to
octave as polyfit (with an additional change to scale x
unconditionally).

Sincerely,

Dmitri.
--

[Ben Abbott]

On Feb 2, 2008, at 3:00 PM, Dmitri A. Sergatskov wrote:

> Though I was partially responsible for all this discussions,
> I am now rethinking this approach.
> polyfit should be as stable numerically as we can make it,
> but I do not think it should try to make a "better"
> fit than the one it was asked for.
>
> Polyfit does fit to a*x^n+b*x^(n-1)... polynomial.
> If people want it to fit to a*(x-x0)^n+b*(x-x0)^(n-1)
> they should pre-condition data accordingly prior to
> passing it to polyfit.
>
> So to make long story short -- I think the original
> Thomas' suggestion (just to normalize the data to max(x))
> was the good one.
>
> On related note -- it appears that wpolifit does a
> better job than polyfit on un-conditioned data as well.
> It appears that the difference is in polyfit
> uses straight "\",  while the wpolyfit does QR decomposition.
> Perhaps it is worth to port this to polyfit as well.
> Frankly I do not see why would not we copy wpolyfit to
> octave as polyfit (with an additional change to scale x
> unconditionally).

Regarding "scale x unconditionally", do you refer to the scaling used  
by wpolyfit;

        (x - mean (x)) / std (x)

or to Thomas' suggestion to just scale the magnitude?

        x / max (x)

If you refer to Thomas' suggestion, the maximum value will result in  
as much trouble/benefit as the minimum value.

Perhaps a better solution would be (a) the geometric mean of the  
magnitudes, (b) the median of the magnitudes, (c) the mean of the  
magnitudes, (d) consider several normalization options and select the  
most numerically stable one.

In any event, what should be done about s.R and s.X? Are they to  
represent the scaled dependent variable?

Ben

[Dmitri A. Sergatskov]

On Feb 2, 2008 2:28 PM, Ben Abbott <bpabbott@...> wrote:

> Regarding "scale x unconditionally", do you refer to the scaling used
> by wpolyfit;
>
>         (x - mean (x)) / std (x)
>
> or to Thomas' suggestion to just scale the magnitude?
>
>         x / max (x)
>

I mean Thomas' suggestion. That is to be precise x / max(abs(x))

> If you refer to Thomas' suggestion, the maximum value will result in
> as much trouble/benefit as the minimum value.

No. If your data are well centered, the min(abs(x)) ~  eps, so scaling to
min does not work. In case of data having a large offset,
min(abs(x)) ~ max(abs(x)) ~ mean(abs(x)), so scaling to any of these
numbers would be equally helpful. But scaling to max(abs(x)) would guarantee
to make all the data in (-1,1) range in all cases and that should help
with numerical precision.

>
> Perhaps a better solution would be (a) the geometric mean of the
> magnitudes, (b) the median of the magnitudes, (c) the mean of the
> magnitudes, (d) consider several normalization options and select the
> most numerically stable one.

See above. I doubt that the fit will be that sensitive to the scaling
parameter, i.e. to say instead of max(abs(x)) you can use 0.5*max(abs(x))
and you probably will not see much of a difference.
Using some more sofisticated approximation is to make some assumption of
the data distribution and we do not want to do that in a generic function.

>
> In any event, what should be done about s.R and s.X? Are they to
> represent the scaled dependent variable?
>

I do not know.

> Ben
>

Regards,

Dmitri.
--

[Ben Abbott]

On Feb 2, 2008, at 3:57 PM, Dmitri A. Sergatskov wrote:

> On Feb 2, 2008 2:28 PM, Ben Abbott <bpabbott@...> wrote:
>
>> Regarding "scale x unconditionally", do you refer to the scaling used
>> by wpolyfit;
>>
>>        (x - mean (x)) / std (x)
>>
>> or to Thomas' suggestion to just scale the magnitude?
>>
>>        x / max (x)
> I mean Thomas' suggestion. That is to be precise x / max(abs(x))
>
>> If you refer to Thomas' suggestion, the maximum value will result in
>> as much trouble/benefit as the minimum value.
>
> No. If your data are well centered, the min(abs(x)) ~  eps, so  
> scaling to
> min does not work. In case of data having a large offset,
> min(abs(x)) ~ max(abs(x)) ~ mean(abs(x)), so scaling to any of these
> numbers would be equally helpful. But scaling to max(abs(x)) would  
> guarantee
> to make all the data in (-1,1) range in all cases and that should help
> with numerical precision.
>

It is been my understanding that the importance is how geometrically  
different the numbers are from 1.

Meaning that eps and 1/eps are equally bad.

Specifically, I thought it would be best if the geometric mean (or  
perhaps median) were at unity.

>>
>> Perhaps a better solution would be (a) the geometric mean of the
>> magnitudes, (b) the median of the magnitudes, (c) the mean of the
>> magnitudes, (d) consider several normalization options and select the
>> most numerically stable one.
>
> See above. I doubt that the fit will be that sensitive to the scaling
> parameter, i.e. to say instead of max(abs(x)) you can use  
> 0.5*max(abs(x))
> and you probably will not see much of a difference.
> Using some more sofisticated approximation is to make some  
> assumption of
> the data distribution and we do not want to do that in a generic  
> function.
>
>>
>> In any event, what should be done about s.R and s.X? Are they to
>> represent the scaled dependent variable?
>>
>
> I do not know.

This is the point of my greatest concern. Since I don't know how this  
output is used, I have no way to determine who it should behave when  
the dependent variable is shifted and/or scaled.

Ben

[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

[Dmitri A. Sergatskov]

On Feb 3, 2008 6:20 PM, Ben Abbott <bpabbott@...> wrote:

> 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 that case can we just leave things as they are and perhaps just add some
words to help file? (May be also pointing to the wpolyfit file).
I am curious to see how well Matlab handles your benchmark.

> 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
>
>

Dmitri.
--

[Ben Abbott]

On Feb 3, 2008, at 7:50 PM, Dmitri A. Sergatskov wrote:

On Feb 3, 2008 6:20 PM, Ben Abbott <bpabbott@...> wrote:

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 that case can we just leave things as they are and perhaps just add some
words to help file? (May be also pointing to the wpolyfit file).
I am curious to see how well Matlab handles your benchmark.

Regarding leaving things as they are, the normalization is not something I'd expect most users to desire to handle for themselves. 

However we might either provide for optional normalization, or provide the routines to do so.

Regarding how Matlab handles this benchmark here are the results.

                       Matlab

 N   M     norm(y-yfit)        norm(p-fit)

 1   2    2.542114973e-13    2.955859084e-12

 1   3    3.410605132e-13    1.693934712e-11

 2   3    1.028152284e-09    1.323761177e-06

 2   4    4.440674381e-10    3.319473118e-06

 2   5    4.939083810e-10    4.947078972e-06

 3   4    2.189135306e-06    5.694889319e+00

 3   5    1.099055345e-06    2.955133100e+00

 3   6    7.348554138e-07    3.296905699e+00

 3   7    1.217871795e-06    5.411753547e+00

 4   5    1.196039913e-03    2.780029492e+06

 4   6    2.441406250e-03    9.003165440e+05

 4   7    3.895778546e-03    4.993247502e+06

 4   8    3.054689407e-03    5.591453444e+05

 4   9    1.708984375e-03    3.790911374e+05

Matlab's results appear consistent with QR without normalization.

Ben

[Rolf Fabian]

Ben Abbott wrote:

I don't believe Octave has any functions that are restricted to real  
values. So ... yes, polyfit accepts and returns complex values.

Ben

Lets say that way: The majority of Octave's mapping functions is
able to handle complex arguments. But unfortunately some
important functions do not and need to be improved.

E.g.:
octave-3.0.0.exe:1> x = fsolve('x.^2+1',0.9999*i)
error: fsolve: iteration is not making good progress
        not returning expected x=+/-i

octave-3.0.0.exe:2> gamma (x=1+i)
error: gamma: unable to handle complex arguments
      Handling complex x requires the Lanczos approximation
      of the gamma function which is  unfortunately
      not implemented in Octave.    

octave-3.0.0.exe:3> erf(1+i)
error: erf: unable to handle complex arguments

octave-3.0.0.exe:4> erfc(1+i)
error: erfc: unable to handle complex arguments

Form 3 of the qz decomposition

Rolf