On Monday, February 24, 2014 3:34:22 AM UTC-8, Hans W Borchers wrote:
>
> The following function is one of the most difficult test functions for
> univariate root finding I know of,
>
> julia> function fn(x::Real)
> return log(x) + x^2/(2*exp(1)) - 2 * x/sqrt(exp(1)) + 1
> end
>
> whose exact root between 1.0 and 3.4 is sqrt(e) = 1.6487212707001282... .
> Applying your nice version of bisect_root to this function returns
>
> julia> bisect_root(fn, 1.0, 3.4)
> (1.648707085670482,1.6487070856704822)
>
> with an absolute error of about 1.5e-05 ! This innocent looking test
> function appears to be so flat that only an arbitrary precision solver will
> locate this root more exactly.
>
> Hans Werner
>
Functions with "flat" zeros are definitely a problem for numerical root
finding. As you suggest, pretty much any "black box" root finder that works
with local float evaluations will lose a bunch of precision here--the issue
isn't specific to bisect_root. Evaluating the output values of bisect_root
gives a clue about what's going on:
julia> fn(1.648707085670482), fn(1.6487070856704822)
(-2.220446049250313e-16,0.0)
The test function evaluated on the upper bound returns 0.0, so as far as
the algorithm can tell, it has found a root. As written, bisect_root will
always do this when you have a sequence of successive floating inputs that
all evaluate to 0.0: it will return the smallest value in the sequence that
evaluates to 0.0 as its upper bound, and the next smallest float as its
lower bound. You could rewrite it to give you the upper bound of the
numerically 0 part instead, or maybe even to find both the upper and lower
edges. The mean of these might be close to the exact root, but it isn't
always.
The example you give isn't really any more or less pathological than other
functions that have coincident 0's of the value and the first and second
derivative. `x -> exp(-x^3)-1` and `x -> (x^3 + 1) - 1` behave about the
same.
I'm actually planning to address exactly this issue in my next post, in the
context of finding extrema of functions instead of zeros. Most functions
are flat near their extrema--only a few are at their zeros.
BTW It would be nice to have Ridders' algorithm available in Julia, too,
> about which the "Numerical Recipes" say:
>
> "*In both reliability and speed, Ridders's method is generally
> competitive with the more highly developed and better established (but more
> complicated) method of Van Wijngaarden, Dekker, and Brent [...].*"
>
> My Julia version of Ridders' method would be as follows, perhaps it can be
> improved with your proposed x1 < (x1 + x2)/2 < x2 :
>
I don't know whether the stopping criterion in your implementation is fine
or not, but it makes me a little suspicious. There aren't any pairs of
floats much larger than 1.0 that are separated by less than eps(), so it
looks to me like if you try to isolate a root larger than 1.0, this
algorithm will run either for maxiter iterations, or until the upper and
lower bounds exactly coincide. Exact coincidence might be fine, or it might
not. What I like about bisection is that it's so simple that I can analyze
it confidently. The *only* arithmetic in bisect_root is a right shift in
the _middle subroutine.
I'd be curious to see ridder's method benchmarked against the bisection
code for a few different objective functions. For simple objective
functions evaluated over floats, it's harder to outperform bisection than a
lot of people think because bisection has such low per-iteration overhead.
The square root in here might hurt it some.
#-- --------------------------------------------------------------------
> function ridders(f::Function, a::Real, b::Real;
> maxiter::Integer = 1000, tol::Real = eps())
>
> x1 = a; x2 = b
> f1 = f(x1); f2 = f(x2)
> if f1 * f2 > 0
> error("f(a) and f(b) must have different signs.")
> elseif f1 == 0
> return x1
> elseif f2 == 0
> return f2
> end
>
> niter = 2
> while niter < maxiter && abs(x1 - x2) > tol
> xm = (x1 + x2)/2.0
> fm = f(xm)
> if fm == 0; return xm; end
>
> x3 = xm + (xm - x1) * sign(f1 - f2) * fm / sqrt(fm^2 - f1 * f2)
> f3 = f(x3)
> niter += 2
> if f3 == 0; return x3; end
>
> if fm * f3 < 0
> x1 = xm; f1 = fm
> x2 = x3; f2 = f3
> elseif f1 * f3 < 0
> x2 = x3; f2 = f3
> elseif f2 * f3 < 0
> x1 = x3; f1 = f3
> else
> error("Inform the maintainer: you should never get here.")
> end
> end
>
> return (x1 + x2) / 2.0
> end
> #-- --------------------------------------------------------------------
>
>
