I am writing a small program that employs Newton's method to solve for
roots of an equation. The function to drive Newton's method is simple;
my trouble lies in finding a reasonable accuracy threshold to stop the
refining iterations.
The problem is floating-point arithmetic accuracy. When the value
"delta" in the function below gets too small, the (decf val delta)
form fails to change the value "val." First, take a look at the
current implementation:
(defun newtons-method (func func-prime initial-val max-iterations
&optional (threshold single-float-negative-epsilon))
(loop repeat max-iterations
with val = initial-val
do (let ((r (funcall func val)))
(when (< (abs r) threshold)
(return val))
(let ((dr (funcall func-prime val)))
(let ((delta (/ r dr)))
(when (< (abs delta) threshold)
(return val))
(decf val delta))))
finally return
(restart-case
(error 'convergence-failure :value val
:threshold threshold
:iterations max-iterations)
(use-final-value ()
:report (lambda (stream)
(format stream "Use final computed value ~A." val))
val)
(use-value (other-val)
:report (lambda (stream)
(format stream "Input another value instead of ~A." val))
:interactive (lambda ()
;; Should the stream be 't' instead?
(format *query-io* "Use instead of ~A: " val)
(multiple-value-list (eval (read))))
other-val))))
With some tracing in place, a sample invocation looks like this:
,----
| initial val for newton's method: -8.700001
| newton's method val: -8.700001
| newton's method r: 0.004221797
| newton's method dr: -6.665623
| newton's method delta: -6.333687E-4
|
| newton's method val: -8.699368
| newton's method r: 1.4066696E-5
| newton's method dr: -6.6235843
| newton's method delta: -2.1237288E-6
|
| newton's method val: -8.699366
| newton's method r: 1.4305115E-6
| newton's method dr: -6.623458
| newton's method delta: -2.1597653E-7
|
| newton's method val: -8.699366
| newton's method r: 1.4305115E-6
| newton's method dr: -6.623458
| newton's method delta: -2.1597653E-7
|
| [Accept the first restart...]
|
| final normal scale: -8.699366
`----
Note how all the values (val and delta in particular) settle in and
don't change through the iterations. I ran some tests (CLISP 2.32) to
see what's going on:
> (type-of -8.699366)
SINGLE-FLOAT
> (type-of -2.1597653E-7)
SINGLE-FLOAT
> (= -8.699366 (- -8.699366 -2.1597653E-7))
T
> (= -8.699366 (- -8.699366 single-float-negative-epsilon))
T
> (< (abs -2.1597653E-7) single-float-negative-epsilon)
NIL
> (< (abs -2.1597653E-7) single-float-epsilon)
NIL
So 2.1597653E-7, which we'll call "delta," is greater than
single-float(-negative)-epsilon, suggesting that delta is large enough
to do useful arithmetic with on single-floats. But subtracting delta
from our stuck "val" -8.699366 does nothing. I must be
misunderstanding how to interpret these epsilon constants.
I'm looking for how to express the threshold such that when delta
becomes so small that adding it to or subtracting it from some other
number produces no observable difference, we'll give up the iterative
refinement. That is, we want the smallest threshold such that
(not (= -8.699366 (- -8.699366 threshold)))
is true.
Advice and clarification would be most welcome.
--
Steven E. Harris :: ········@raytheon.com
Raytheon :: http://www.raytheon.com
>>>>> "Steven" == Steven E Harris <········@raytheon.com> writes:
>> (= -8.699366 (- -8.699366 -2.1597653E-7))
Steven> T
>> (= -8.699366 (- -8.699366 single-float-negative-epsilon))
Steven> T
>> (< (abs -2.1597653E-7) single-float-negative-epsilon)
Steven> NIL
>> (< (abs -2.1597653E-7) single-float-epsilon)
Steven> NIL
Steven> So 2.1597653E-7, which we'll call "delta," is greater than
Steven> single-float(-negative)-epsilon, suggesting that delta is large enough
Steven> to do useful arithmetic with on single-floats. But subtracting delta
Steven> from our stuck "val" -8.699366 does nothing. I must be
Steven> misunderstanding how to interpret these epsilon constants.
Yes, I think you are misinterpreting this epsilon constant.
single-float-epsilon is the smallest number such that, when added to
1.0, you get a different number when properly rounded. Your number is
not 1.0, but 8.699.... So adding some number larger than
single-float-epsilon may or may not change it.
A relative threshold may work for you, though. Something like if
delta/root > single-float-epsilon, continue refining.
But I think I'd choose 2 or 5 or 10 times single-float-epsilon as the
threshold. It's pretty hard to get all intermediate calculations of
delta accurate to that many digits.
But IANA numerical analyst.
Ray
On Mon, 19 Jan 2004 12:29:39 -0800
"Steven E. Harris" <········@raytheon.com> wrote:
> I am writing a small program that employs Newton's method to solve for
> roots of an equation. The function to drive Newton's method is simple;
> my trouble lies in finding a reasonable accuracy threshold to stop the
> refining iterations.
>
> The problem is floating-point arithmetic accuracy. When the value
> "delta" in the function below gets too small, the (decf val delta)
> form fails to change the value "val." First, take a look at the
> current implementation:
>
>
> (defun newtons-method (func func-prime initial-val max-iterations
> &optional (threshold
> single-float-negative-epsilon))
> (loop repeat max-iterations
> with val = initial-val
> do (let ((r (funcall func val)))
> (when (< (abs r) threshold)
> (return val))
> (let ((dr (funcall func-prime val)))
> (let ((delta (/ r dr)))
> (when (< (abs delta) threshold)
> (return val))
> (decf val delta))))
> finally return
> (restart-case
> (error 'convergence-failure :value val
> :threshold threshold
> :iterations max-iterations)
> (use-final-value ()
> :report (lambda (stream)
> (format stream "Use final computed
> value ~A." val))
> val)
> (use-value (other-val)
> :report (lambda (stream)
> (format stream "Input another value instead
> of ~A." val))
> :interactive (lambda ()
> ;; Should the stream be 't' instead?
> (format *query-io* "Use instead of ~A:
> " val)(multiple-value-list (eval
> (read))))
> other-val))))
>
>
>
> With some tracing in place, a sample invocation looks like this:
>
> ,----
> | initial val for newton's method: -8.700001
> | newton's method val: -8.700001
> | newton's method r: 0.004221797
> | newton's method dr: -6.665623
> | newton's method delta: -6.333687E-4
> |
> | newton's method val: -8.699368
> | newton's method r: 1.4066696E-5
> | newton's method dr: -6.6235843
> | newton's method delta: -2.1237288E-6
> |
> | newton's method val: -8.699366
> | newton's method r: 1.4305115E-6
> | newton's method dr: -6.623458
> | newton's method delta: -2.1597653E-7
> |
> | newton's method val: -8.699366
> | newton's method r: 1.4305115E-6
> | newton's method dr: -6.623458
> | newton's method delta: -2.1597653E-7
> |
> | [Accept the first restart...]
> |
> | final normal scale: -8.699366
> `----
>
> Note how all the values (val and delta in particular) settle in and
> don't change through the iterations. I ran some tests (CLISP 2.32) to
> see what's going on:
>
> > (type-of -8.699366)
> SINGLE-FLOAT
>
> > (type-of -2.1597653E-7)
> SINGLE-FLOAT
>
> > (= -8.699366 (- -8.699366 -2.1597653E-7))
>
>
> > (= -8.699366 (- -8.699366 single-float-negative-epsilon))
>
>
> > (< (abs -2.1597653E-7) single-float-negative-epsilon)
> NIL
>
> > (< (abs -2.1597653E-7) single-float-epsilon)
> NIL
>
> So 2.1597653E-7, which we'll call "delta," is greater than
> single-float(-negative)-epsilon, suggesting that delta is large enough
> to do useful arithmetic with on single-floats. But subtracting delta
> from our stuck "val" -8.699366 does nothing. I must be
> misunderstanding how to interpret these epsilon constants.
>
> I'm looking for how to express the threshold such that when delta
> becomes so small that adding it to or subtracting it from some other
> number produces no observable difference, we'll give up the iterative
> refinement. That is, we want the smallest threshold such that
>
> (not (= -8.699366 (- -8.699366 threshold)))
>
> is true.
>
> Advice and clarification would be most welcome.
As Raymond Toy suggests you are mistaking the meaning. You have about
seven digits of precision with an (IEEE-754) single precision float.
Imagine adding 1e-7 to 1e7. It's not going to change the value.
Steven E. Harris wrote:
> I am writing a small program that employs Newton's method to solve for
> roots of an equation. The function to drive Newton's method is simple;
> my trouble lies in finding a reasonable accuracy threshold to stop the
> refining iterations.
>
> The problem is floating-point arithmetic accuracy. When the value
> "delta" in the function below gets too small, the (decf val delta)
> form fails to change the value "val." First, take a look at the
> current implementation:
>
>
> (defun newtons-method (func func-prime initial-val max-iterations
> &optional (threshold single-float-negative-epsilon))
> (loop repeat max-iterations
> with val = initial-val
> do (let ((r (funcall func val)))
> (when (< (abs r) threshold)
> (return val))
> (let ((dr (funcall func-prime val)))
> (let ((delta (/ r dr)))
> (when (< (abs delta) threshold)
> (return val))
> (decf val delta))))
> finally return
> (restart-case
> (error 'convergence-failure :value val
> :threshold threshold
> :iterations max-iterations)
> (use-final-value ()
> :report (lambda (stream)
> (format stream "Use final computed value ~A." val))
> val)
> (use-value (other-val)
> :report (lambda (stream)
> (format stream "Input another value instead of ~A." val))
> :interactive (lambda ()
> ;; Should the stream be 't' instead?
> (format *query-io* "Use instead of ~A: " val)
> (multiple-value-list (eval (read))))
> other-val))))
>
>
>
> With some tracing in place, a sample invocation looks like this:
>
> ,----
> | initial val for newton's method: -8.700001
> | newton's method val: -8.700001
> | newton's method r: 0.004221797
> | newton's method dr: -6.665623
> | newton's method delta: -6.333687E-4
> |
> | newton's method val: -8.699368
> | newton's method r: 1.4066696E-5
> | newton's method dr: -6.6235843
> | newton's method delta: -2.1237288E-6
> |
> | newton's method val: -8.699366
> | newton's method r: 1.4305115E-6
> | newton's method dr: -6.623458
> | newton's method delta: -2.1597653E-7
> |
> | newton's method val: -8.699366
> | newton's method r: 1.4305115E-6
> | newton's method dr: -6.623458
> | newton's method delta: -2.1597653E-7
> |
> | [Accept the first restart...]
> |
> | final normal scale: -8.699366
> `----
>
> Note how all the values (val and delta in particular) settle in and
> don't change through the iterations. I ran some tests (CLISP 2.32) to
> see what's going on:
>
>
>>(type-of -8.699366)
>
> SINGLE-FLOAT
>
>
>>(type-of -2.1597653E-7)
>
> SINGLE-FLOAT
>
>
>>(= -8.699366 (- -8.699366 -2.1597653E-7))
>
> T
>
>
>>(= -8.699366 (- -8.699366 single-float-negative-epsilon))
>
> T
>
>
>>(< (abs -2.1597653E-7) single-float-negative-epsilon)
>
> NIL
>
>
>>(< (abs -2.1597653E-7) single-float-epsilon)
>
> NIL
>
> So 2.1597653E-7, which we'll call "delta," is greater than
> single-float(-negative)-epsilon, suggesting that delta is large enough
> to do useful arithmetic with on single-floats. But subtracting delta
> from our stuck "val" -8.699366 does nothing. I must be
> misunderstanding how to interpret these epsilon constants.
>
> I'm looking for how to express the threshold such that when delta
> becomes so small that adding it to or subtracting it from some other
> number produces no observable difference, we'll give up the iterative
> refinement. That is, we want the smallest threshold such that
>
> (not (= -8.699366 (- -8.699366 threshold)))
>
> is true.
>
> Advice and clarification would be most welcome.
Well, I do not know anything about this stuff, but that has never
stopped me before:
(loop with trials = 200
repeat trials
for base = (random 10.0)
when (= base (- base -2.1597653E-7))
count 1 into eqls
finally (print (/ eqls trials)))
...returns under ACL6.2/Win32 values between .5+ to .7+
I am thinking about the butterfly effect story. I think it depends on
what internal value you get from -8.699366, and how -8.699366 is to
being an integral multiple of the epsilon.
I dare say (* 1.5 epsilon) or twice epsilon would always work.
How's that for guesswork?
kenny
>
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
Kenny Tilton wrote:
>
>
> Steven E. Harris wrote:
>
>> I am writing a small program that employs Newton's method to solve for
>> roots of an equation. The function to drive Newton's method is simple;
>> my trouble lies in finding a reasonable accuracy threshold to stop the
>> refining iterations.
>>
>> The problem is floating-point arithmetic accuracy. When the value
>> "delta" in the function below gets too small, the (decf val delta)
>> form fails to change the value "val." First, take a look at the
>> current implementation:
>>
>>
>> (defun newtons-method (func func-prime initial-val max-iterations
>> &optional (threshold
>> single-float-negative-epsilon))
>> (loop repeat max-iterations
>> with val = initial-val
>> do (let ((r (funcall func val)))
>> (when (< (abs r) threshold)
>> (return val))
>> (let ((dr (funcall func-prime val)))
>> (let ((delta (/ r dr)))
>> (when (< (abs delta) threshold)
>> (return val))
>> (decf val delta))))
>> finally return
>> (restart-case
>> (error 'convergence-failure :value val
>> :threshold threshold
>> :iterations max-iterations)
>> (use-final-value ()
>> :report (lambda (stream)
>> (format stream "Use final computed
>> value ~A." val))
>> val)
>> (use-value (other-val)
>> :report (lambda (stream)
>> (format stream "Input another value instead
>> of ~A." val))
>> :interactive (lambda ()
>> ;; Should the stream be 't' instead?
>> (format *query-io* "Use instead of ~A:
>> " val)
>> (multiple-value-list (eval (read))))
>> other-val))))
>>
>>
>>
>> With some tracing in place, a sample invocation looks like this:
>>
>> ,----
>> | initial val for newton's method: -8.700001
>> | newton's method val: -8.700001
>> | newton's method r: 0.004221797
>> | newton's method dr: -6.665623
>> | newton's method delta: -6.333687E-4
>> | | newton's method val: -8.699368
>> | newton's method r: 1.4066696E-5
>> | newton's method dr: -6.6235843
>> | newton's method delta: -2.1237288E-6
>> | | newton's method val: -8.699366
>> | newton's method r: 1.4305115E-6
>> | newton's method dr: -6.623458
>> | newton's method delta: -2.1597653E-7
>> | | newton's method val: -8.699366
>> | newton's method r: 1.4305115E-6
>> | newton's method dr: -6.623458
>> | newton's method delta: -2.1597653E-7
>> | | [Accept the first restart...]
>> | | final normal scale: -8.699366
>> `----
>>
>> Note how all the values (val and delta in particular) settle in and
>> don't change through the iterations. I ran some tests (CLISP 2.32) to
>> see what's going on:
>>
>>
>>> (type-of -8.699366)
>>
>>
>> SINGLE-FLOAT
>>
>>
>>> (type-of -2.1597653E-7)
>>
>>
>> SINGLE-FLOAT
>>
>>
>>> (= -8.699366 (- -8.699366 -2.1597653E-7))
>>
>>
>> T
>>
>>
>>> (= -8.699366 (- -8.699366 single-float-negative-epsilon))
>>
>>
>> T
>>
>>
>>> (< (abs -2.1597653E-7) single-float-negative-epsilon)
>>
>>
>> NIL
>>
>>
>>> (< (abs -2.1597653E-7) single-float-epsilon)
>>
>>
>> NIL
>>
>> So 2.1597653E-7, which we'll call "delta," is greater than
>> single-float(-negative)-epsilon, suggesting that delta is large enough
>> to do useful arithmetic with on single-floats. But subtracting delta
>> from our stuck "val" -8.699366 does nothing. I must be
>> misunderstanding how to interpret these epsilon constants.
>>
>> I'm looking for how to express the threshold such that when delta
>> becomes so small that adding it to or subtracting it from some other
>> number produces no observable difference, we'll give up the iterative
>> refinement. That is, we want the smallest threshold such that
>>
>> (not (= -8.699366 (- -8.699366 threshold)))
>>
>> is true.
>>
>> Advice and clarification would be most welcome.
>
>
> Well, I do not know anything about this stuff, but that has never
> stopped me before:
>
> (loop with trials = 200
> repeat trials
> for base = (random 10.0)
> when (= base (- base -2.1597653E-7))
> count 1 into eqls
> finally (print (/ eqls trials)))
>
> ...returns under ACL6.2/Win32 values between .5+ to .7+
>
> I am thinking about the butterfly effect story. I think it depends on
> what internal value you get from -8.699366, and how -8.699366 is to
> being an integral multiple of the epsilon.
>
> I dare say (* 1.5 epsilon) or twice epsilon would always work.
>
> How's that for guesswork?
Not too good. Seems like 8 (not 7.999) time the epsilon is needed to get
reliable results. A power of two, I like that! But I still need a new
theory...
:)
kenny
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
Kenny Tilton wrote:
> But I still need a new
> theory...
..or to read the doc:
> The value of each of the constants short-float-epsilon,
> single-float-epsilon, double-float-epsilon, and long-float-epsilon is
> the smallest positive float of the given format, such that the
> following expression is true when evaluated:
(not (= (float 1 > epsilon) (+ (float 1 epsilon) epsilon)))
:(
kt
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
From: Joe Marshall
Subject: Re: Floating point arithmetic (epsilon)
Date:
Message-ID: <smib77ow.fsf@comcast.net>
"Steven E. Harris" <········@raytheon.com> writes:
> I'm looking for how to express the threshold such that when delta
> becomes so small that adding it to or subtracting it from some other
> number produces no observable difference, we'll give up the iterative
> refinement. That is, we want the smallest threshold such that
>
> (not (= -8.699366 (- -8.699366 threshold)))
>
> is true.
Why not just check if the value changes?
If you don't want to do that, though, your threshold will have to be a
ratio rather than an absolute value because it will depend upon your
answer. Another alternative would be testing
(< (abs delta) (* threshold val))
--
~jrm
Joe Marshall <·············@comcast.net> writes:
> Why not just check if the value changes?
I hadn't thought of that. I just implemented it and it works fine,
irrespective of threshold. But it feels a little weird reading code
that subtracts one non-zero number from another and then tests if
they're the same. Actually, I guess it's really testing whether
"non-zero" is sufficiently so.
> If you don't want to do that, though, your threshold will have to be
> a ratio rather than an absolute value because it will depend upon
> your answer. Another alternative would be testing
>
> (< (abs delta) (* threshold val))
I tried that, and it worked to a point, but I still found that I had
to scale single-float-negative-epsilon by 4 or so to avoid the trap,
per Raymond Toy's suggestion. First I tried no scaling and got
stuck. So I tried a factor of 2 and got stuck again. I couldn't make
it stick with a factor of 4. If someone asked me where the magic
factor came from, all I could say is, "I just kept trying examples
until they didn't break any more." I'm not comfortable with that
approach to my work.
I've been reading up on search hits for "floating point math epsilon"
and keep finding the same simple advice: Equality tests can't be
trusted; instead you want something like
(defun sufficiently-equal (a b)
(< (abs (- a b)) +some-epsilon+))
which is pretty much the same as my original "delta" test. The problem
-- largely unaddressed -- is coming up with the right epsilon.
--
Steven E. Harris :: ········@raytheon.com
Raytheon :: http://www.raytheon.com
In article <···············@l75001820.us.ray.com>,
"Steven E. Harris" <········@raytheon.com> writes:
> ...
> I've been reading up on search hits for "floating point math epsilon"
> and keep finding the same simple advice: Equality tests can't be
> trusted; instead you want something like
>
> (defun sufficiently-equal (a b)
> (< (abs (- a b)) +some-epsilon+))
>
> which is pretty much the same as my original "delta" test. The problem
> -- largely unaddressed -- is coming up with the right epsilon.
if you have an idea about the range your results are in, and if this
range isn't too wide, it should be relatively easy. with a very wide
range you run, for any epsilon, the risk that at one end the numbers
are so small that thos test succeeds where it shouldn't, while at the
other end they are so large that the epsilon doesn't effect the
values. you probably would be better off with something like
(defun sufficiently-equal (a b)
(< (abs (- 1 (/ a b)) +some-epsilon+)))
where the epsilon would depend on the precision used for the
calculation (you would have to make some provisions for b being 0)
hs
--
Not everything that counts can be counted,
not everything that can be counted counts.
A. Einstein
* Steven E. Harris
| So 2.1597653E-7, which we'll call "delta," is greater than
| single-float(-negative)-epsilon, suggesting that delta is large
| enough to do useful arithmetic with on single-floats.
No, this is very wrong. 1.0�epsilon is the smallest representable
number different from 1.0. To calculate the smallest representable
number different from any other value, you have to scale epsilon.
This should be instructive:
(rational (- 1.0 single-float-negative-epsilon))
=> 16777215/16777216
(rational (+ 1.0 single-float-epsilon))
=> 8388609/8388608
To find the epsilon of any given floating-point number, you have to
find the number that would be a change to the least significant bit
away from the number you look at. This will be a scaled +epsilon when
the exponent does not change, and a scaled -epsilon when it does.
The function DECODE-FLOAT returns a number of type (float 1/b (1.0)),
where b is the floating point representation base, which you can
always find with (/ (decode-float 1.0)), which means that in order to
scale the epsilon correctly, you need to compute the scaling factor:
(defun scaled-epsilon (float)
(let ((scale (- (nth-value 1 (decode-float float))
(nth-value 1 (decode-float 1.0)))))
(scale-float single-float-epsilon scale)))
To show how this works, consider the following dialogue:
(scaled-epsilon 8.5)
=> 4.768372e-7
(integer-decode-float (- 8.5 *))
=> 8912895
=> -20
=> 1
(integer-decode-float (+ 8.5 **))
=> 8912897
=> -20
=> 1
(integer-decode-float 8.5)
=> 8912896
=> -20
=> 1
This takes care of +epsilon, but -epsilon is a different matter, since
it applies only when the value in question is exactly representable
with a single bit and the operation yields a smaller value. Let me
therefore consider the full range:
(defun scaled-epsilon (float &optional (operation '+))
"Return the smallest number that would return a value different from
FLOAT if OPERATION were applied to FLOAT and this number. OPERATION
should be either + or -, and defauls to +."
(multiple-value-bind (significand exponent)
(decode-float float)
(multiple-value-bind (1.0-significand 1.0-exponent)
(decode-float (float 1.0 float))
(if (and (eq operation '-)
(= significand 1.0-significand))
(scale-float (typecase float
(short-float short-float-negative-epsilon)
(single-float single-float-negative-epsilon)
(double-float double-float-negative-epsilon)
(long-float long-float-negative-epsilon))
(- exponent 1.0-exponent))
(scale-float (typecase float
(short-float short-float-epsilon)
(single-float single-float-epsilon)
(double-float double-float-epsilon)
(long-float long-float-epsilon))
(- exponent 1.0-exponent))))))
| I must be misunderstanding how to interpret these epsilon constants.
I hope this helps.
--
Erik Naggum | Oslo, Norway
Act from reason, and failure makes you rethink and study harder.
Act from faith, and failure makes you blame someone and push harder.
Erik Naggum <····@naggum.no> writes:
> I hope this helps.
Definitely. Thank you very much. I spent the morning experimenting
with your examples and reading about decode-float and scale-float, and
learned more about floating-point representation in the process.
I integrated scaled-epsilon into my Newton's method function and it
works well to trap "delta" when it gets too small. Now I have a
different problem, and I suspect this one has little to do with Lisp,
and more to do with Newton's method. I'll ask here as a clarifying
process.
First, consider that almost every description of Newton's method I've
read mentions iterating until the maximum iteration count is met or
the solution gets "close enough." Now, seeing that the whole point of
using Newton's method is to find the solution, "close enough" can't
mean the distance between the candidate root and the real root since
we don't know the real root ahead of time. Instead, it usually means
that F(t), or "r" in my code below, is "close enough" to zero. Unlike
my previous question about an appropriate epsilon with respect to a
given number, I don't see an obvious way to determine what "close
enough to zero" means.
I tried to just avoid this problem with the following reasoning. With
Erik's epsilon solution in hand, we have a deterministic way to know
when "delta" -- the horizontal adjustment of the candidate root -- is
too small to make any significant difference to the answer. Any actual
solution should also see "delta" decrease below epsilon, possibly to
zero. So our deterministic "delta" check should be enough. The only
trade-off is that we wind up doing a few extra calculations and
comparisons on the final iteration.
Here's the current version of newtons-method:
(defun newtons-method (func func-prime initial-val max-iterations)
(loop repeat max-iterations
with val = initial-val
do (let ((r (funcall func val)))
(let ((threshold (scaled-epsilon val '-)))
(when (< (abs r) threshold) ; *** Shouldn't be
(return val)) ; *** necessary.
(let ((dr (funcall func-prime val)))
(let ((delta (/ r dr)))
(when (< (abs delta) threshold)
(return val))
(decf val delta)))))
finally return
(restart-case
(error 'convergence-failure :value val :threshold threshold :iterations max-iterations)
(use-final-value ()
:report (lambda (stream)
(format stream "Use final computed value ~A." val))
val)
(use-value (other-val)
:report (lambda (stream)
(format stream "Input another value instead of ~A." val))
:interactive (lambda ()
(format *query-io* "Use instead of ~A: " val)
(multiple-value-list (eval (read))))
other-val))))
Note the two lines marked with the "***" comments. Per my argument
above, this check for a small enough "r" or F(t) should not be
necessary; the delta v. threshold comparison lower down should be
sufficient. Also, note that I'm somewhat arbitrarily comparing "r" to
the scaled-epsilon threshold. The reasoning behind this particular
threshold (as a function of the current candidate root, or "val" here)
doesn't apply to "r."
But I'm finding test cases where the "r" check is necessary. Consider
this run (with the "***" lines above disabled), where delta does not
grow small enough below epsilon, but is large enough to keep pushing
us past the solution and bouncing back again. That is, the solution
must lay between 23.32404 and 23.324041, but some error among r, dr,
and the resultant delta prevents convergence. Note that the function
is not of the cubic variety where loops like these are expected.
,----
| initial val for newton's method: 8.00001
| newton's method val: 8.00001
| newton's method threshold: 4.768372E-7
| newton's method r: 3.1204147
| newton's method dr: -0.5893294
| newton's method delta: -5.2948565
|
| newton's method val: 13.294867
| newton's method threshold: 4.768372E-7
| newton's method r: 1.2106586
| newton's method dr: -0.21821134
| newton's method delta: -5.548101
|
| newton's method val: 18.842968
| newton's method threshold: 9.536744E-7
| newton's method r: 0.36898565
| newton's method dr: -0.10429378
| newton's method delta: -3.537945
|
| newton's method val: 22.380913
| newton's method threshold: 9.536744E-7
| newton's method r: 0.06392646
| newton's method dr: -0.07103838
| newton's method delta: -0.8998862
|
| newton's method val: 23.2808
| newton's method threshold: 9.536744E-7
| newton's method r: 0.0028018951
| newton's method dr: -0.064934544
| newton's method delta: -0.04314953
|
| newton's method val: 23.32395
| newton's method threshold: 9.536744E-7
| newton's method r: 5.841255E-6
| newton's method dr: -0.06465997
| newton's method delta: -9.033805E-5
|
| newton's method val: 23.32404
| newton's method threshold: 9.536744E-7
| newton's method r: 1.1920929E-7
| newton's method dr: -0.06465941
| newton's method delta: -1.8436496E-6
|
| newton's method val: 23.324041
| newton's method threshold: 9.536744E-7
| newton's method r: -1.1920929E-7
| newton's method dr: -0.06465939
| newton's method delta: 1.8436501E-6
|
| newton's method val: 23.32404
| newton's method threshold: 9.536744E-7
| newton's method r: 1.1920929E-7
| newton's method dr: -0.06465941
| newton's method delta: -1.8436496E-6
|
| [Accept the first restart...]
|
| final normal scale: 23.32404
`----
Again, note that delta never gets small enough for the epsilon
trap. "r" has fallen below our epsilon threshold, but that seems
incidental; the epsilon threshold is computed with respect to "val,"
while we have no criteria to decide what a "small enough" or "too
small" r is. Every implementation example I've found online compares r
to some constant EPSILON, which here might mean
single-float-epsilon. Even that check would still fail.
By inspection, we can see this looping with no progress. That suggests
that one approach would be to keep the second-to-last "val" around and
compare it to the current "val." If they match, we may be in one of
these reflective looping scenarios.
I interpret the loop above to mean, geometrically, that the slope dr
is too shallow and the projected tangent line winds up cutting through
the curve before hitting the F(t) axis. Cutting off "r" below some
(arbitrary) threshold works to hide this problem, but lacks the
reasoning Erik brought to bear with the scaled delta epsilon.
Again, any suggestions and guidance would be welcome.
--
Steven E. Harris :: ········@raytheon.com
Raytheon :: http://www.raytheon.com
>>>>> "Steven" == Steven E Harris <········@raytheon.com> writes:
Steven> First, consider that almost every description of Newton's method I've
Steven> read mentions iterating until the maximum iteration count is met or
Steven> the solution gets "close enough." Now, seeing that the whole point of
Steven> using Newton's method is to find the solution, "close enough" can't
Steven> mean the distance between the candidate root and the real root since
Steven> we don't know the real root ahead of time. Instead, it usually means
Steven> that F(t), or "r" in my code below, is "close enough" to zero. Unlike
Steven> my previous question about an appropriate epsilon with respect to a
Steven> given number, I don't see an obvious way to determine what "close
Steven> enough to zero" means.
I think you should look in netlib.org for some root-finding algorithms
to see how numerical analysts do this.
But consider finding the root of (x-pi)^50. Or (x-pi)^50+1d-20. I
think it's pretty hard in general to know you're at a root, unless you
actually know something about the function you're working with.
Ray
+ "Steven E. Harris" <········@raytheon.com>:
| Unlike my previous question about an appropriate epsilon with
| respect to a given number, I don't see an obvious way to determine
| what "close enough to zero" means.
With iterative schemes, when I want to iterate to machine accuracy (or
until numerical flaws in my code limit the accuracy), I often find it
useful to just iterate until the residual stops decreasing. It's
simple and effective, though perhaps not totally kosher. Remember
that Newton's method, in the usual case, roughly doubles the number of
correct bits with each iteration, so convergence is extremely fast
once you get reasonably close to a root. (Maybe a forum dedicated to
numerical methods would be a better place to ask?)
By the way, rather than using Newton's method, I find a method based
on inverse quadratic interpolation easy to use and quite robust. It
requires no derivative. It works this way: Given three points
x[n-2], x[n-1], x[n] with corresponding function values f(x[i]),
determine the coefficients of a quadratic polynomial approximating the
/inverse/ of f by interpolating at the given points, then evaluate
this at zero to get x[n+1]. However you may need to use linear
interpolation in some cases. I found some C code (zeroin.c) in netlib
years ago which explains it. (Hmm. Is it called Brent's algorithm?)
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- Debating gives most of us much more psychological satisfaction
than thinking does: but it deprives us of whatever chance there is
of getting closer to the truth. -- C.P. Snow
In article <···············@l75001820.us.ray.com>,
"Steven E. Harris" <········@raytheon.com> writes:
> Erik Naggum <····@naggum.no> writes:
>
>> I hope this helps.
>
> Definitely. Thank you very much. I spent the morning experimenting
> with your examples and reading about decode-float and scale-float, and
> learned more about floating-point representation in the process.
>
> I integrated scaled-epsilon into my Newton's method function and it
> works well to trap "delta" when it gets too small. Now I have a
> different problem, and I suspect this one has little to do with Lisp,
> and more to do with Newton's method. I'll ask here as a clarifying
> process.
Newton's method is pretty tricky when your function has multiple
roots. it is good when all roots are single. i once used the trick
exploiting the fact that polynomials with multiple roots share roots
with their derivatives. by dividing the polynomial by its derivative
i ended up with one that was guaranteed to have only single roots, so
that newton could be applied without concerns
> ...
hs
--
Not everything that counts can be counted,
not everything that can be counted counts.
A. Einstein
From: Joe Marshall
Subject: Re: Floating point arithmetic (epsilon)
Date:
Message-ID: <vfn5tkdz.fsf@ccs.neu.edu>
"Steven E. Harris" <········@raytheon.com> writes:
> Again, any suggestions and guidance would be welcome.
You are entering into the territory of heavy-duty numerical analysis.
Contrary to popular belief, it is more tedious than difficult.
Floating point operations are allowed to give slightly erroneous
results if the mathematically correct result is not representable.
IEEE 754 specifies that results for individual floating point
operations should appear to have been computed to infinite precision
and then rounded to the nearest representable number. The floating
point result for any single operation should not be more than 1 unit
in the least significant position (ULP).
Presumably, the function of which you are finding the roots involves a
few floating point operations. You can analyze the error that may be
produced, but for the sake of argument, let's assume that the bottom 2
bits may be erroneous.
I do not know where FUNC-PRIME is coming from, but let us assume that
it too may have errors in the bottom two bits. (If you are doing
this:
(defun derivative (func)
(lambda (x dx)
(/ (- (funcall func (+ x dx)) (funcall func x)) dx)))
and (funcall func (+ x dx)) is erroneous by being too high and
(funcall func x) is erroneous by being too low, then
(- (funcall func (+ x dx)) (funcall func x)) will have twice the error
of either. If (funcall func (+ x dx)) returns f(x+dx)+e0 (where e0
is the error), and (funcall func x) returns f(x)+e1, the the
difference will be (f(x+dx)-f(x))+(e0+e1). As dx approaches 0,
f(x+dx) - f(x) will approach 0, but (e0+e1) will stay the same. If
the difference term gets small enough, the error term will be large
enough to make a difference.)
In any case, your value of R has an error term and your value of dR
has an error term, so your delta will have an error term. (The value
of the error term will depend on R and dR.)
> But I'm finding test cases where the "r" check is necessary. Consider
> this run (with the "***" lines above disabled), where delta does not
> grow small enough below epsilon, but is large enough to keep pushing
> us past the solution and bouncing back again. That is, the solution
> must lay between 23.32404 and 23.324041, but some error among r, dr,
> and the resultant delta prevents convergence.
The mathematical delta gets below epsilon, but the actual delta (which
includes the error term) does not.
Fixing this is an interesting problem. There is no solution that is
independent of FUNC and FUNC-PRIME, but there are solutions that may
be good enough for your purposes. Single-precision floating point
(such as you are using) gives you only about 7 decimal digits in the
best case. Erik's epsilon scaling gives you the value of 1 binary
digit in the least significant place. If you multiply this by a power
of two, you can get your algorithm to disregard the least significant
digits for convergence. The appropriate power of two will depend on
FUNC and FUNC-PRIME.
> I interpret the loop above to mean, geometrically, that the slope dr
> is too shallow and the projected tangent line winds up cutting through
> the curve before hitting the F(t) axis. Cutting off "r" below some
> (arbitrary) threshold works to hide this problem, but lacks the
> reasoning Erik brought to bear with the scaled delta epsilon.
This is a different (and also interesting) possibility. Since you are
using floating point, you cannot compute F at an arbitrary point, and
F(t) takes on discrete values that approximate the mathematical
value. The mathematical function F may cross the axis where the
floating point F does not.
Note, too, that as dR gets near zero, the error term in dR becomes
more significant.
If you need a truly accurate result, you will need to tediously
analyze the function and the errors from start through finish. If you
don't want to do this, simply raise the threshold until it works.
Additionally, you should use double precision floating-point where
possible because it will reduce the error terms by several orders of
magnitude.
In article <···············@L75001820.us.ray.com>,
Steven E. Harris <········@raytheon.com> wrote:
>Erik Naggum <····@naggum.no> writes:
>
>> I hope this helps.
>
>Definitely. Thank you very much. I spent the morning experimenting
>with your examples and reading about decode-float and scale-float, and
>learned more about floating-point representation in the process.
Hats to hearts for DEC and IBM-mainframe style floating point
representation.
*Brother Lupo plays taps on the pipes*
Say, why are you using clisp and single precision anyway? There is very
little overhead in double precision code which fits in the cache like
Newton's method should. In fact, I can't think of too many reasons to use
a single precision float *anywhere* in the 21st century. And there is very
little point in using byte compiled code for numerics when you have access
to nice compilers like cmu-cl, unless you just plain like clisp and don't
care about speed.
In fact, what are you trying to do anyway? Is this just a learning
exercise?
>First, consider that almost every description of Newton's method I've
>read mentions iterating until the maximum iteration count is met or
>the solution gets "close enough." Now, seeing that the whole point of
>using Newton's method is to find the solution, "close enough" can't
>mean the distance between the candidate root and the real root since
>we don't know the real root ahead of time. Instead, it usually means
>that F(t), or "r" in my code below, is "close enough" to zero. Unlike
>my previous question about an appropriate epsilon with respect to a
>given number, I don't see an obvious way to determine what "close
>enough to zero" means.
It sort of depends what the problem is, doesn't it? Close enough to zero
means close enough to zero for you, personally, but not so close to zero
that you can't represent it in your machine without being awfully clever.
Things you should know:
Newton-Raphson method doesn't work unless you have a good first guess;
i.e. it doesn't converge globally; for that you need to be more clever.
If you calculate a naieve numeric derivative, aka something like
(defun func-prime (func x dx)
(\ (- (func (+ x dx)) (func x)) dx))
-it's also extremely susceptable to numeric instability. I suspect this
was the case; no? SICP makes this mistake, BTW, so you are in good
company.
The way in which you calculate the function in question can also get you
into trouble with numeric stability.
Newton-Raphson is also non optimal for certain kinds of multivariate
problem, though it is kinda kewl in that it can be used in more general
cases.
>I tried to just avoid this problem with the following reasoning. With
>Erik's epsilon solution in hand, we have a deterministic way to know
>when "delta" -- the horizontal adjustment of the candidate root -- is
>too small to make any significant difference to the answer. Any actual
>solution should also see "delta" decrease below epsilon, possibly to
>zero.
Think not of zero when you use floating point code. In many
FP representations, including IEEE, which is probably what your computer
uses, *there is no zero.*
I'm overstating this a little, but just a little (actually, zero can be
positive or negative, which is, IMO an abomination).
http://en.wikipedia.org/wiki/IEEE_Floating_Point_Standard
A common conceptual mistake, one implicitly made by people who should know
better (like, actual computer scientists or applied math guys), is that
computer floating point numbers are the same thing as real numbers, in the
same sense that computer integers represent mathematical integers (which
they do, perfectly, assuming the integers fit in a register somewhere). FP
numbers are no such thing. FP is a convenient fiction that may or may not
actually represent a real number in a way that is useful.
You should really get and read a copy of "Numerical Recipes in Fortran" or
whatever the modern analog of Press, Teukolsky, Vetterling and Flannery's
cookbook is. While it is no guarantee you'll do the right thing, it's not
a bad place to start. Someone at Raytheon must have a copy of Horowitz and
Hill as well -they have useful words on FP representations.
You could also probably use a previously canned root-finder from netlib
if having a sophisticated working piece of code is a priority.
Particularly if you can get the Fortran-Lisp converter to work (if so,
let me know, as I could not easily get it to work). There is no dishonor
in using canned floating point code, and a lot of advantages, especially
if you aren't familiar with its pitfalls.
-Lupo
"Always be willing to speak your mind, and a base man will avoid you."
-William Blake <··@io.com>
From: Ray Dillinger
Subject: Re: Floating point arithmetic (epsilon)
Date:
Message-ID: <400F539B.F5DB8170@sonic.net>
Lupo LeBoucher wrote:
>
> Say, why are you using clisp and single precision anyway? There is very
> little overhead in double precision code which fits in the cache like
> Newton's method should. In fact, I can't think of too many reasons to use
> a single precision float *anywhere* in the 21st century.
I have little use for designations like "single", "float", "double",
"long double", etc. I think it's better style to use macros to
refer directly to a number of bits and get a float type that has
*at least* that many bits.
> >I tried to just avoid this problem with the following reasoning. With
> >Erik's epsilon solution in hand, we have a deterministic way to know
> >when "delta" -- the horizontal adjustment of the candidate root -- is
> >too small to make any significant difference to the answer. Any actual
> >solution should also see "delta" decrease below epsilon, possibly to
> >zero.
>
> Think not of zero when you use floating point code. In many
> FP representations, including IEEE, which is probably what your computer
> uses, *there is no zero.*
> I'm overstating this a little, but just a little (actually, zero can be
> positive or negative, which is, IMO an abomination).
> http://en.wikipedia.org/wiki/IEEE_Floating_Point_Standard
I understand about branch cuts and how impractical it is to get
multiple answers and do something useful with them, but signed
zeros and picking a single value for equations with multiple
solutions has always seemed like a copout to me. And in my own
mind I distinguish exact zero (which is neither positive nor
negative) from inexact zero (which will have one sign or
another). Hmmmm.
I wonder... How stupid would it get if all the basic math
functions were defined with range and domain including lists
of numbers (including potentially complex numbers) as well
as single numbers? You'd run the same code, but you'd get
lists of possible values back, having split wherever there
was an ambiguity.
Bear
> I wonder... How stupid would it get if all the basic math
> functions were defined with range and domain including lists
> of numbers (including potentially complex numbers) as well
> as single numbers? You'd run the same code, but you'd get
> lists of possible values back, having split wherever there
> was an ambiguity.
>
> Bear
What would you do with function that could return infinitely many
values, like inverse sin? Use streams rather than lists?
Cheers,
rif
··@io.com (Lupo LeBoucher) writes:
> better (like, actual computer scientists or applied math guys), is that
> computer floating point numbers are the same thing as real numbers, in the
> same sense that computer integers represent mathematical integers (which
> they do, perfectly, assuming the integers fit in a register somewhere). FP
> numbers are no such thing. FP is a convenient fiction that may or may not
> actually represent a real number in a way that is useful.
Indeed. We must keep in mind that "floating point" is actually a
representation of a subset of rational numbers!
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/