3

I have a function f(x), that is positive and decreasing for x<c, and is zero for all x>=c. How can I find c, the threshold where the function hits zero (to within a tolerance)?

Here's an example:

zer = function(x){
    ifelse(x>5, rep(0,length(x)), 5 - x)
}

> x=-5:15
> plot(x,zer(x))

enter image description here

You can use uniroot to find where a function crosses zero, but that relies on the function being negative and positive on either side of the crossing, so I can't use that here. In the above case it evaluates zer at the upper interval (15), finds a zero, and returns that.

I wrote a bisection algorithm that starts with the interval and moves left if f(midpoint) == 0 or right if f(midpoint) > 0. This works, but I'm wondering if I've missed an implementation in a package or elsewhere that does this better, or if I've missed "one simple trick that lets you use uniroot to solve this".

The best I can find in the docs is the Numerical Task View's cryptic "There are implementations of the bisection algorithm in several contributed packages."

Note I don't have the gradient of f(x) so can't use Newton's method or anything that needs gradient evaluations at x values.

Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • Is the discontinuity always due to an `ifelse`. If so, you could just find the root of the decreasing function (in this case `uniroot(function(x) 5-x, range(x))`) – dww Nov 25 '17 at 15:51
  • No the discontinuity could come from anything. The motivating example for me was the area remaining in a district after removing anything within a buffer distance of the roads in the district. At some buffer distance the buffer covers the whole district, and for any value beyond that the area is zero. Computing the buffer is expensive so I wanted a solution with as few evaluations as possible, which is a general requirement for optimisation and root finding. – Spacedman Nov 25 '17 at 17:26

3 Answers3

2

One possibility - instead of returning 0 for f(x)==0, return a small constant negative number:

zer2 = function(x){
    y = zer(x)
    y[y==0]=-1e-6
    y
}

this gives a solution that can be found with uniroot:

> uniroot(zer2, c(-5,15))
$root
[1] 5.000043

The size of the small negative number might be important though.

Also, I'm not sure how well-behaved uniroot is if half the function is a constant value of -1 - it seems to cope in this case and its probably robust enough.

Spacedman
  • 92,590
  • 12
  • 140
  • 224
2

This problem seems to be well suited to a bisection method. We can do it something like this for example:

findroot = function(f, x, tol = 1e-9) {
  if (f(x[2]) > 0) stop("No root in range")
  if (f(x[1]) <= tol) return (x[1])
  else {
    xmid = mean(x)
    if(f(xmid) == 0) x = c(x[1], xmid) else x = c(xmid, x[2])
    return (findroot(f, x, tol))
  }
}

findroot(zer, range(x))
# [1] 5
dww
  • 30,425
  • 5
  • 68
  • 111
  • As I said, I have implemented a bisection method for this, but I was after something better. `uniroot` uses a different method (possibly a secant method) which gives better efficiency (needs fewer function evals). – Spacedman Nov 25 '17 at 17:23
  • Didn't read the question carefully enough to notice you'd already tried a bisection method. Secant & false position methods may not converge for these discontinuous functions. In fact, I tried implementing the false position method for your example, and in this case it does not converge. Secant method will suffer the same fate, as it uses the same method to find a bifurcation point as does false position. – dww Nov 25 '17 at 17:55
  • 1
    N.B. Speed of convergenge is only one desireable propoerty for optimizations. Robustness is also important, which bisection will give you here. Many methods will not be robust on discontinuous functions. – dww Nov 25 '17 at 18:02
0

I'm not sure this answers the question, but I've tried to follow your problem description to the letter:

  1. f(x) is positive and decreasing for x < c, and is zero for all x >= c;
  2. find its zero to within a tolerance.

So I start with a little helper function, a variant of a function I generally use in order to avoid FAQ 7.31. Then it's a simple matter of keeping the first TRUE it returns.

is.zero <- function(x, tol = .Machine$double.eps^0.5) abs(x) < tol

i <- which(is.zero(zer(x)))[1]
x[i]
#[1] 5

Off-topic.
The helper function mentioned is

is.equal <- function(x, y, tol = .Machine$double.eps^0.5) abs(x - y) < tol

This function is not equivalent to all.equal since it is vectorized and the values of the shorter argument x or y recycle.

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • 2
    That finds the first value in x where f(x) is zero, but that depends on the values in `x`. The `x` in my example is only meant for plotting (perhaps I should have plotted as a line) but f is defined at all continuous x values. To find x such that f(x) = 0 using this method to a precision of `x` requires an `x` vector with a separation equal to the tolerance, eg seq(-5,15,by=0.00001), which means a lot of function evaluations. – Spacedman Nov 25 '17 at 14:49
  • @Spacedman You are right, I just tested it with your new `seq` and my solution is 1200 times slower. Seemed so promising... – Rui Barradas Nov 25 '17 at 16:49