Find simple binary decision boundary of heaviside step function using binary search

findDecisionBoundary.ts

export interface FinderOption {
  f: (x: number) => boolean
  p: number
  epsilon: number
}

// This is theoretically pure, but impractical since this may result stack overflow due to heavy recursion
export function findDecisionBoundary(option: FinderOption, lower: number, upper: number, n: number = 0): { 
  boundary?: number 
  n: number
} {
  const { f, p, epsilon } = option 
  const m = lower * p + upper * (1 - p)
  const fm = f(m)

  if (upper - lower <= epsilon) {
    if (fm) return { boundary: m, n }
    return { n }
  }
  if (fm) return findDecisionBoundary(option, lower, m + epsilon * p / 2, n + 1)
  return findDecisionBoundary(option, m, upper, n + 1)
}

export function findDecisionBoundaryIterative({ f, p, epsilon }: FinderOption, lower: number, upper: number): { 
  boundary?: number
  n: number
} {
  let n = 0
  let prev = 1e+10
  do {
    const m = lower * p + upper * (1 - p)
    const fm = f(m)
    const diff = upper - lower

    // convergence detection
    if (diff <= epsilon || prev === diff) {
      if (fm) return { boundary: m, n }
      return { n }
    }
    prev = diff

    if (fm) {
      upper = m + epsilon * p / 2
    } else {
      lower = m
    }
    n += 1
  } while (n <= 1024) // prevent infinite loop
  
  return { boundary: lower * p + upper * (1 - p), n }
}

This is an algorithm searching for one dimensional simple binary decision boundary.

Let f be the function between R and {0, 1} where f(x) = x >= B ? 1 : 0.
Here B is unknown so that we want to find it efficiently.

Our computer has limitation of discretization, so it'll try to find an open set (L, U) includes B and approximates B using linear interpolation of L and U. Here the objective U - L ≤ ε.

By using above algorithm with L = inf(dom(f)) and U = sup(dom(f)), I observed that empirical average complexity goes O(-lg ε) with respect to precision ε.
Maybe it is due to similarity between this and classical binary search.

Note that p / 2 in the upper = m + epsilon * p / 2 must be preserved. This guarantees that algorithm never fall into infinite loop.
Theoretically any positive non zero number below p is ok. (explanation in Korean. I'm still working in progress)

Now the optimal p is 0.5 which is pretty intuitive, but I can't prove this as I'm not an expert computer scientist. But by conducting a simple experiment, I observed that p = 0.5 makes the lowest iterations in average. Roughly speaking, using binary search dividing not equally has fairly bad edge case (which the wanted value is placed on the opposite side), and same thing works with real numbers.

function getAverageExecutions(p: number, epsilon: number): number {
  let sum = 0
  let cnt = 0
  for (let E = 0.1; E <= 0.9; E += 0.005) {
    sum += findDecisionBoundaryIterative({ f: (x) => x >= 100 + E * 100, p, epsilon }, 100, 200).n
    cnt += 1
  }
  return sum / cnt
}

function format(x: number): string {
  return `${x}`.slice(0, 4)
}

const epsilon = 2 ** -24
for (let p = 0.001; p < 1; p += 0.001) {
  const n = getAverageExecutions(p, epsilon)
  console.log(`\t${p}\t${n}`)
}