Floating point inverses and stability
Let f be a monotone, strictly convex function on a real interval I and let g be its inverse. For example, we could have f(x) = ex and g(x) = log x.
Now suppose we round our results to N digits. That is, instead of working with f and g we actually work with fN and gN where
fN(x) = round(f(x), N)
and
gN(x) = round(g(x), N)
and round(y, N) is the number y rounded to N significant figures [1].
This is what happens when we implement our functions f and g in floating point arithmetic. We don't actually get the values of f and g but the values of fN and gN.
We assumed that f and g are inverses, but in general fN and gN will not be exact inverses. And yet in some sense the functions fN and gN are like inverses. Harold Diamond [2] proved that if go back and forth applying fN and gN two times, after two round trips the values quit changing.
To make this more precise, define
hN(x) = gN( fN(x)).
In general, hN(x) does not equal x, but we do have the following:
hN( hN( hN(x) ) ) = hN( hN(x) ).
The value we get out of hN(x) might not equal x, but after we've applied hN twice, the value stays the same if we apply hN more times.
Connection to connectionsDiamond's stability theorem looks a lot like a theorem about Galois connections. My first reaction was that Diamond's theorem was simply a special case of a more general theorem about Galois connections, but it cannot.
A pair of monotone functions F and G form a Galois connection if for all a in the domain of F and for all b in the domain of G,
F(a) b a G(b).
Let F and G form a Galois connection and define
H(x) = G( F(x) ).
Then
H( H(x) ) = H(x).
This result is analogous to Diamond's result, and stronger. It says we get stability after just one round trip rather than two.
The hitch is that although the functions f and g form a Galois connection, the functions fN and gN do not. Nevertheless, Diamond proved that fN and gN form some sort of weaker numerical analog of a Galois connection.
ExampleThe following example comes from [2]. Note that the example rounds to two significant figures, not two decimal places.
from decimal import getcontext, Decimal # round to two significant figures getcontext().prec = 2 def round(x): return float( Decimal(x) + 0 ) def f(x): return 115 - 35/(x - 97) def f2(x): return round(f(x)) def g(x): return 97 + 35/(115 - x) def g2(x): return round(g(x)) def h2(x): return g2(f2(x)) N = 110 print(h2(N), h2(h2(N)), h2(h2(h2(N))))
This prints
100.0 99.0 99.0
showing that It shows that the function h2 satisfies Diamond's theorem, but it does not satisfy the identify above for Galois compositions. That is, we stabilize after two round trips but not after just one round trip.
Related posts[1] Our digits" here need not be base 10 digits. The stability theorem applies in any radix b provided bN >= 3.
[2] Harold G. Diamond. Stability of Rounded Off Inverses Under Iteration. Mathematics of Computation, Volume 32, Number 141. January 1978, pp. 227-232.
The post Floating point inverses and stability first appeared on John D. Cook.