On Monday, February 24, 2014 8:16:28 PM UTC-8, Stefan Karpinski wrote:
>
> This sort of comparison seems inevitably bound to favor bisection for
> objective functions that are fast to evaluate and to favor "sophisticated"
> methods when the objective function is costly, no?
>
Yes, that's right. And also if you want to work to work to arbitrary
precision, you will definitely want a method with a higher order of
convergence. And if you are working in more than one dimension.
But there are definitely interesting questions to ask about objective
functions that cost less than 100 transcendental function evaluations, and
I think it's useful to know the ballpark of the crossover.
> On Mon, Feb 24, 2014 at 11:01 PM, Jason Merrill
> <jwme...@gmail.com<javascript:>
> > wrote:
>
>> On Monday, February 24, 2014 7:48:50 PM UTC-8, Jason Merrill wrote:
>>
>>> On Monday, February 24, 2014 3:34:22 AM UTC-8, Hans W Borchers wrote:
>>>
>>>> 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.
>>>
>>
>> Because I can't leave well enough alone...
>>
>> julia> function ridder_loop()
>> accum = 0.0
>> for i = 1:10000
>> accum += ridders(sin, 2.0, 4.0)
>> end
>> accum
>> end
>>
>> julia> function bisect_loop()
>> accum = 0.0
>> for i = 1:10000
>> accum += bisect_root(sin, 2.0, 4.0)[1]
>> end
>> accum
>> end
>>
>> julia> @time ridder_loop()
>> elapsed time: 10.430907316 seconds (1596537204 bytes allocated)
>> 31415.92653589205
>>
>> julia> @time ridder_loop()
>> elapsed time: 10.440722588 seconds (1596480048 bytes allocated)
>> 31415.92653589205
>>
>> julia> @time bisect_loop()
>> elapsed time: 0.152900673 seconds (25984456 bytes allocated)
>> 31415.92653589205
>>
>> julia> @time bisect_loop()
>> elapsed time: 0.144590689 seconds (25920064 bytes allocated)
>> 31415.92653589205
>>
>> So for this example, admittedly with an overly simplistic objective
>> function, bisection is winning by a factor of about 70.
>>
>>
>>
>>
>>> #-- --------------------------------------------------------------------
>>>> 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
>>>> #-- ------------------------------------------------------------
>>>> --------
>>>>
>>>>
>
