Binary search on real space

Question

I want to perform a binary search on a continuous unimodal function ƒ(x)=y, where x and y are real numbers. I'm not looking up values in an array, and so I don't have clean integer inputs that I'm stepping along.

My original attempt at a function (written in JavaScript) has the problem that when the input value you are looking for happens to sit on a boundary chosen by the algorithm, the algorithm continues to run until the numeric precision is exhausted:

/**
 * y: the target output value to find a matching input
 * minX: the smallest input value
 * maxX: the largest input value
 * ƒ:    a function that returns a value given an X
 * ε:    the threshold to compare outputs versus the target (default:0)
 */
function binarySearch(y, minX, maxX, ƒ, ε) {
    if (ε===undefined) ε=0;
    let m=minX, n=maxX, k, v, Δ;
    while (m<=n) {
        k = (n+m)/2;
        v = ƒ(k);
        Δ = y-v;
        if (Math.abs(Δ)<=ε) return k;
        if (Δ>0) m = k;
        else     n = k;
    }
    if (Math.abs(y-ƒ(m))<=ε) return m;
    if (Math.abs(y-ƒ(n))<=ε) return n;
}

With the above, the call binarySearch( 0, 0, 10, n=>n ) will run 1077 iterations until m=0 and n=5e-324 before (m+n)/2 is finally so close to 0 that even with ε=0 the JavaScript interpreter cannot tell the difference.

A hack I was going to use is to provide a minimum step function that modifies each boundary by a fixed amount (similar to array index +/- 1). This forces the boundary to move faster, but also requires an ε>0 in case the boundary overshoots the value. It feels gross:

// step: a minimum amount to move each boundary each time
function binarySearch(y, minX, maxX, ƒ, ε, step) {
    if (ε===undefined) ε=0;
    if (step===undefined) step=0;
    let m=minX, n=maxX, k, v, Δ;
    while (m<=n) {
        k = (n+m)/2;
        v = ƒ(k);
        Δ = y-v;
        if (Math.abs(Δ)<=ε) return k;
        if (Δ>0) m = k+step;
                 n = k-step;
    }
    if (Math.abs(y-ƒ(m))<=ε) return m;
    if (Math.abs(y-ƒ(n))<=ε) return n;
}

A ~clean fix is to check the boundaries on every pass by moving the final two if statements into the while loop. This calls ƒ() three times as often each pass, and so seems inelegant.

/**
 * y: the target output value to find a matching input
 * minX: the smallest input value
 * maxX: the largest input value
 * ƒ:    a function that returns a value given an X
 * ε:    the threshold to compare outputs versus the target (default:0)
 */
function binarySearch(y, minX, maxX, ƒ, ε) {
    if (ε===undefined) ε=0;
    let m=minX, n=maxX, k, v, Δ;
    while (m<=n) {
        k = (n+m)/2;
        v = ƒ(k);
        Δ = y-v;
        if (Math.abs(y-ƒ(m))<=ε) return m;
        if (Math.abs(y-ƒ(n))<=ε) return n;
        if (Math.abs(Δ)<=ε) return k;
        if (Δ>0) m = k;
        else     n = k;
    }
}

It smells to me like there ought to be an elegant solution, some sort of fencepost I'm not thinking of, that fixes this efficiently and elegantly.

Answer 1

I realize now that the only time this is a problem is when the value I'm looking for is either of the starting boundaries. In all other cases, a boundary is only moved after it is checked.

So, all I had to do was check both boundaries once, before the loop starts. A fixed, slightly terser version of the original function:

/**
 * y: the target output value to find a matching input
 * minX: the smallest input value
 * maxX: the largest input value
 * ƒ:    a function that returns a value `y` given an `x`
 * ε:    the threshold to compare outputs versus the target (default:0)
 */
function binarySearch(y, minX, maxX, ƒ, ε) {
    if (ε===undefined) ε=0;
    let m=minX, n=maxX, k, Δ;
    if (Math.abs(y-ƒ(m))<=ε) return m;
    if (Math.abs(y-ƒ(n))<=ε) return n;
    while (m<n) {
        Δ = y - ƒ( k = (n+m)/2 );
        if (Math.abs(Δ)<=ε) return k;
        if (Δ>0) m=k;
        else     n=k;
    }
}

Note that this function also fixes a ~bug in the while loop by changing the condition to m<n instead of m<=n.

Answer 2

Something to consider is that the computer representation of the real field is not uniform. Half of all the possible unique values for a float (or double) occur in the interval 1 to -1. This usually doesn't cause a problem, but if your algorithm can choose a value close to the max or min values it can produce surprises.

I suspect the reason your initial implementation had problems was that your calculations were getting jitter in the lowest bit positions.