Floating point arithmetic (epsilon)

From: Steven E. Harris (seharris_at_raytheon.com)
Date: 01/19/04 

Date: Mon, 19 Jan 2004 12:29:39 -0800

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        :: seharris@raytheon.com
Raytheon                :: http://www.raytheon.com

Relevant Pages

  • Re: [PATCH] hrtimer: increase clock min delta threshold while interrupt hanging
    ... tick_dev_program_eventunless min delta has a very low value... ... more agressive one for the timer tick only? ... Since embedded systems could need high resolution ... Actually perhaps 10 iterations would be more reasonable.... ...
    (Linux-Kernel)
  • Re: numeric issue: 0.95 - 0.05 == 0.8499996
    ... Precisely I proceede as follows: ... float e = 1.0f; ... size (delta = 0.05) times the number of iterations from the original number: ...
    (comp.lang.java.programmer)
  • Re: Optimized Shadow Volume Rendering
    ... delta was employed, at a delta of 2,3 or 4 then the shadow started to ... Jack ... > (n frames below threshold could aggregate to significant distance). ... >> problems for high end systems that update every frame for a light that ...
    (microsoft.public.win32.programmer.directx.graphics)