Home > Archive > Fortran > April 2007 > error analysis for brute force polynomial eval?
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
error analysis for brute force polynomial eval?
|
|
| Paul van Delst 2007-04-18, 8:04 am |
| Hi all,
Anyone have any links, refs, etc for error analysis in polynomial evaluations? I've got
Higham's book that provides the error analysis algorithm for Horner's method, but I would
like to do something similar for code that looks like:
x2=x*x
x3=x2*x
x4=x3*x
x5=x4*x
y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
I don't got the mathematical moxie to roll my own error bound analysis.
The above snippet is an example of code I've received and when I refactored to use
Horner's method I got quite different results, e.g. for one case:
y(horner) = 81.61805014616661
y(brute) = 81.61953466293143
(All double precision variables and coeffs too.)
Thanks for any info.
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
| |
| Beliavsky 2007-04-18, 7:05 pm |
| On Apr 18, 8:51 am, Paul van Delst <Paul.vanDe...@noaa.gov> wrote:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial evaluations? I've got
> Higham's book that provides the error analysis algorithm for Horner's method, but I would
> like to do something similar for code that looks like:
Probably someone in comp.lang.fortran can answer, but the place for
numerical analysis questions is sci.math.num-analysis .
| |
| Gordon Sande 2007-04-18, 7:05 pm |
| On 2007-04-18 09:51:07 -0300, Paul van Delst <Paul.vanDelst@noaa.gov> said:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial
> evaluations? I've got Higham's book that provides the error analysis
> algorithm for Horner's method, but I would like to do something similar
> for code that looks like:
>
> x2=x*x
> x3=x2*x
> x4=x3*x
> x5=x4*x
> y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
>
> I don't got the mathematical moxie to roll my own error bound analysis.
>
> The above snippet is an example of code I've received and when I
> refactored to use Horner's method I got quite different results, e.g.
> for one case:
>
> y(horner) = 81.61805014616661
> y(brute) = 81.61953466293143
>
> (All double precision variables and coeffs too.)
>
> Thanks for any info.
>
> cheers,
>
> paulv
Think of it as an inner product of the vector of coefficients with
the vector of powers. I would be surprised if there is no analysis
of error bounds for inner products. The only bother is that you will
need to derive the errors in the powers. But that should be an even
earlier result.
It requires patience rather than moxie to work your way back through
the several levels.
Usually this sort of request would draw snide comments about doing
your own homework. But I think you are past the homework stage.
| |
| Tim Prince 2007-04-18, 7:05 pm |
| Paul van Delst wrote:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial
> evaluations? I've got Higham's book that provides the error analysis
> algorithm for Horner's method, but I would like to do something similar
> for code that looks like:
>
> x2=x*x
> x3=x2*x
> x4=x3*x
> x5=x4*x
> y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
>
> I don't got the mathematical moxie to roll my own error bound analysis.
>
> The above snippet is an example of code I've received and when I
> refactored to use Horner's method I got quite different results, e.g.
> for one case:
>
> y(horner) = 81.61805014616661
> y(brute) = 81.61953466293143
>
> (All double precision variables and coeffs too.)
>
As optimizations are permitted by Fortran standard, you can't count on
the evaluation order in your "brute force" expression. If you are at
all interested in predictable numerical behavior, use something more
like Horner.
| |
| Paul van Delst 2007-04-18, 7:05 pm |
| Gordon Sande wrote:
> On 2007-04-18 09:51:07 -0300, Paul van Delst <Paul.vanDelst@noaa.gov> said:
>
>
> Think of it as an inner product of the vector of coefficients with
> the vector of powers. I would be surprised if there is no analysis
> of error bounds for inner products. The only bother is that you will
> need to derive the errors in the powers. But that should be an even
> earlier result.
Thanks, I'll continue digging.
> It requires patience rather than moxie to work your way back through
> the several levels.
>
> Usually this sort of request would draw snide comments about doing
> your own homework. But I think you are past the homework stage.
Cripes, I hope so! :o) I will admit to a modicum of laziness in asking the question,
though. The refs I do have are full of proofs and lemmas and what not.... not my forte.
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
| |
| Paul van Delst 2007-04-18, 7:05 pm |
| Tim Prince wrote:
> Paul van Delst wrote:
>
> As optimizations are permitted by Fortran standard, you can't count on
> the evaluation order in your "brute force" expression. If you are at
> all interested in predictable numerical behavior, use something more
> like Horner.
Hi,
I do use Horner by default (and not just because the adjoint form is ridiculously simple
compared to the same for the brute force method). I was looking for some numerical
ammunition so that I could eventually put in a request for changing the original code
(which is used to compute ocean permittivities that are used in a microwave ocean
emissivity model that is itself used in radiative transfer model to compute satellite
radiances so satobs can be used in a data assimilation system for NWP. Phewph.) Changing
code for a "solved" problem can engender a bit of resistance unless there's a well defined
reason.
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
| |
| Beliavsky 2007-04-18, 7:05 pm |
| On Apr 18, 9:35 am, Gordon Sande <g.sa...@worldnet.att.net> wrote:
> On 2007-04-18 09:51:07 -0300, Paul van Delst <Paul.vanDe...@noaa.gov> said:
>
>
>
>
>
>
>
>
>
>
>
>
> Think of it as an inner product of the vector of coefficients with
> the vector of powers. I would be surprised if there is no analysis
> of error bounds for inner products. The only bother is that you will
> need to derive the errors in the powers. But that should be an even
> earlier result.
>
> It requires patience rather than moxie to work your way back through
> the several levels.
Maybe the SEA ("Simple Error Arithmetic in Fortran 90/95") program at
http://shum.huji.ac.il/~agay/err.f90/download.html can do this.
Some compilers offer quadruple precision and/or REAL*10, and there
exist packages that implement multiple precision in software. The
question is whether evaluating the polynomials using these methods
slows down the program excessively.
| |
| Gordon Sande 2007-04-18, 7:05 pm |
| On 2007-04-18 11:27:44 -0300, Paul van Delst <Paul.vanDelst@noaa.gov> said:
> Tim Prince wrote:
>
> Hi,
>
> I do use Horner by default (and not just because the adjoint form is
> ridiculously simple compared to the same for the brute force method). I
> was looking for some numerical ammunition so that I could eventually
> put in a request for changing the original code (which is used to
> compute ocean permittivities that are used in a microwave ocean
> emissivity model that is itself used in radiative transfer model to
> compute satellite radiances so satobs can be used in a data
> assimilation system for NWP. Phewph.) Changing code for a "solved"
> problem can engender a bit of resistance unless there's a well defined
> reason.
>
> cheers,
>
> paulv
The monomials are about the worst possible basis if you are interested in
error bounds. If you want to go full press you could expand in terms of
Tchebychev (or any of the several dozen alternate spellings) polynomials
after you have shifted the range of definition.
What you want are smallish coefficents that do not depend on cancellation.
That says you are looking for an orthogonal basis. You want the cancellation
to be early and exact. (Like voting in Chicago - early and often!)
More likely sources of error are the determination, or even reporting,
of the coefficients rather than the numerical error in evaluation of the
polynomials.
A quick experiment is to drop/change/etc the last digit of the coefficients
so the cancellation is less accurrate. Watch the resulting error bounce around!
| |
| Ron Shepard 2007-04-18, 7:05 pm |
| In article <f0547t$p4h$1@news.nems.noaa.gov>,
Paul van Delst <Paul.vanDelst@noaa.gov> wrote:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial
> evaluations? I've got
> Higham's book that provides the error analysis algorithm for Horner's method,
> but I would
> like to do something similar for code that looks like:
>
> x2=x*x
> x3=x2*x
> x4=x3*x
> x5=x4*x
> y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
>
> I don't got the mathematical moxie to roll my own error bound analysis.
>
> The above snippet is an example of code I've received and when I refactored
> to use
> Horner's method I got quite different results, e.g. for one case:
>
> y(horner) = 81.61805014616661
> y(brute) = 81.61953466293143
>
> (All double precision variables and coeffs too.)
This is probably due simply to the order of operations. If you
change your expression to
y = a5*x5 + a4*x4 + a3*x3 + a2*x2 + a1*x + a0
or
y = a0 + (a1*x + (a2*x2 + (a3*x3 + (a4*x4 + a5*x5))))
then you will probably get the same result as horner's method. Your
difference is probably due to the fact that some of the terms are
positive and some are negative, and you get different truncation
errors when you add the terms in a different order.
Another possible reason is that your compiler may be using extended
precision with one of you expressions but not with the other.
However, count the total number of operations. In your original
expression there are four initial multiplications followed by five
more multiplications and five additions to evaluate y. If you use
horner's method directly
y = a0 + x * (a1 + x * (a2 + x * (a3 + (x * (a4 + x * a5)))))
there are only five multiplications and five additions. Of course,
the compiler may recognize some of this and rearrange things for
you, but why write complicated inefficient code in the first place.
Horner's method (which is just factoring out the common expressions
from the above expression for y) is both simple and efficient. If
you were to evaluate the expression by hand, you would naturally use
horner's method, right, even if you didn't know its name?
There are other ways to evaluate polynomials, which have different
errors, but the above is probably enough for you to understand what
is causing the difference you see. It might be interesting to
evaluate your expression at higher precision (either in fortran or
in some other language, such as mathematica) just to see which
result is the most accurate.
$.02 -Ron Shepard
| |
| robert.corbett@sun.com 2007-04-19, 4:08 am |
| On Apr 18, 5:51 am, Paul van Delst <Paul.vanDe...@noaa.gov> wrote:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial evaluations? I've got
> Higham's book that provides the error analysis algorithm for Horner's method, but I would
> like to do something similar for code that looks like:
>
> x2=x*x
> x3=x2*x
> x4=x3*x
> x5=x4*x
> y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
>
> I don't got the mathematical moxie to roll my own error bound analysis.
>
> The above snippet is an example of code I've received and when I refactored to use
> Horner's method I got quite different results, e.g. for one case:
>
> y(horner) = 81.61805014616661
> y(brute) = 81.61953466293143
>
> (All double precision variables and coeffs too.)
>
> Thanks for any info.
What are the values of X and the coefficients?
If you are interested in finding bounds for this particular
computation,
you could try using interval arithmetic.
Bob Corbett
| |
| Arjen Markus 2007-04-19, 4:08 am |
| On 18 apr, 14:51, Paul van Delst <Paul.vanDe...@noaa.gov> wrote:
> Hi all,
>
> Anyone have any links, refs, etc for error analysis in polynomial evaluations? I've got
> Higham's book that provides the error analysis algorithm for Horner's method, but I would
> like to do something similar for code that looks like:
>
> x2=x*x
> x3=x2*x
> x4=x3*x
> x5=x4*x
> y = a0 + a1*x + a2*x2 + a3*x3 +a4*x4 + a5*x5
>
> I don't got the mathematical moxie to roll my own error bound analysis.
>
> The above snippet is an example of code I've received and when I refactored to use
> Horner's method I got quite different results, e.g. for one case:
>
> y(horner) = 81.61805014616661
> y(brute) = 81.61953466293143
>
> (All double precision variables and coeffs too.)
>
> Thanks for any info.
>
> cheers,
>
> paulv
>
> --
> Paul van Delst Ride lots.
> CIMSS @ NOAA/NCEP/EMC Eddy Merckx
Just a thought:
Have you tried with different compiler settings for floating-point
operations?
This might turn out to be one area where storing intermediate results
in
extended precision makes a significant difference.
Regards,
Arjen
| |
| Arjen Markus 2007-04-23, 8:05 am |
| One reference I found via Google notes that the problem of evaluating
a polynomial function is numerically unstable if you are trying to
evaluate
it near a root - you get a cancellation of terms in that case, making
it
very hard to get the "right" answer.
Of course, these roots need not to be real! Can you examine whether
this is the case for your problem? If it is, you have a serious
problem -
no matter what the numbers you get now.
This is demonstrated with a real root in the following program:
program polyn
implicit none
real(kind=kind(1.0d00)), dimension(5) :: a
real(kind=kind(1.0d00)) :: x, y1, y2, y3
integer :: i
!
! (1 + 2.1*x + 3.5*x**2 + 100.0**3) * (x-11.2)
!
a = (/ 0.0d0, 1.0d0, 2.1d0, 3.5d0, 100.0d0 /) + &
-11.2d0 * (/ 1.0d0, 2.1d0, 3.5d0, 100.0d0, 0.0d0 /)
do i = 1,10
x = 11.2 + (i-6)*0.01
y3 = x**4.0 * a(5)
y1 = a(1) + a(2) * x + a(3) * x**2 + a(4) * x**3 + a(5) * x**4
y3 = y3 + a(1) + x * a(2)
y2 = a(1) + x * ( a(2) + x * ( a(3) + x * ( a(4) + x *
a(5) ) ) )
y3 = y3 + x**2 * (a(3) + a(4) * x )
write(*,*) x, y1, y2, y3
enddo
endprogram
The results are not bad at first sight with double precision (though I
do
get differences in the last two digits), but dramatic with single
precision.
I did not tune this equation at all - just chose a few numbers.
Regards,
Arjen
| |
| Paul van Delst 2007-04-23, 7:05 pm |
| Arjen Markus wrote:
> One reference I found via Google notes that the problem of evaluating
> a polynomial function is numerically unstable if you are trying to
> evaluate
> it near a root - you get a cancellation of terms in that case, making
> it
> very hard to get the "right" answer.
>
> Of course, these roots need not to be real! Can you examine whether
> this is the case for your problem? If it is, you have a serious
> problem -
> no matter what the numbers you get now.
Well, the polynomials and coefficients are from a paper for computing ocean permittivities
(a complex quantity). Just looking at the coefficient values for the various polynomials,
they seem to be well behaved; i.e. each successive coeff is smaller (typically by an order
of magnitude) and of opposite sign to the previous coeff.
I have not, however, checked to see where the roots are for a typical range of input (sea
surface temperatures).
Thanks for the tip. I will check it out. (This stuff is like an onion! :o)
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
| |
| SimonG 2007-04-25, 8:04 am |
| On Apr 23, 8:25 am, Arjen Markus <arjen.mar...@wldelft.nl> wrote:
....
> a = (/ 0.0d0, 1.0d0, 2.1d0, 3.5d0, 100.0d0 /) + &
> -11.2d0 * (/ 1.0d0, 2.1d0, 3.5d0, 100.0d0, 0.0d0 /)
....
Doesn't compile with g95 because it's illegal to have + -11
Simon Geard
| |
| Arjen Markus 2007-04-27, 7:05 pm |
| On 25 apr, 11:52, SimonG <s...@whiteowl.co.uk> wrote:
> On Apr 23, 8:25 am, Arjen Markus <arjen.mar...@wldelft.nl> wrote:
>
> ...
>
>
> ...
>
> Doesn't compile with g95 because it's illegal to have + -11
>
> Simon Geard
Ah! That is why! CVF was perfectly happy, but I noticed the same error
when I tried
it on Linux. Thanks for noticing this.
Regards,
Arjen
| |
| jamesgiles@att.net 2007-04-27, 7:05 pm |
| On Apr 27, 7:37 am, Arjen Markus <arjen.mar...@wldelft.nl> wrote:
> On 25 apr, 11:52, SimonG <s...@whiteowl.co.uk> wrote:
....
....[color=darkred]
....[color=darkred]
> Ah! That is why! CVF was perfectly happy, but I noticed the same error
> when I tried
> it on Linux. Thanks for noticing this.
CVF has an option to turn this off (it's the one that enforces
standard
compliance). You have to watch this if you don't turn that feature
off. It will accept:
X = Y ** -J * Z
But it will parse that as:
X= Y ** (-J * Z)
This is almost never what the programmer expects.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
| |
| SimonG 2007-04-27, 7:05 pm |
| On Apr 27, 2:37 pm, Arjen Markus <arjen.mar...@wldelft.nl> wrote:
> On 25 apr, 11:52, SimonG <s...@whiteowl.co.uk> wrote:
>
>
>
>
>
>
>
> Ah! That is why! CVF was perfectly happy, but I noticed the same error
> when I tried
> it on Linux. Thanks for noticing this.
>
> Regards,
>
> Arjen
Actually it was noticed by Andy Vaught when I reported it as a g95
bug!
Simon
|
|
|
|
|