Finding multiple roots of univariate functions in R

In this post I describe some methods for root finding in R and the limitations of these when there is more than one root. As an example consider the following function

f <- function(x) cos(x)^4 - 4*cos(x)^3 + 8*cos(x)^2 - 5*cos(x) + 1/2

Here the function is plotted and we can see it has two roots on this interval.

Perhaps the most widely used root finding function in R is uniroot. However, this function has major limitations if there is more than one root (as implied by its name!).

uniroot(f, c(0,2))

Error in uniroot(f, c(0, 2)) : 
 f() values at end points not of opposite sign

In the case when there are an odd number of roots on the interval, uniroot will only find and return one of them. This is a serious problem arising from the underlying mathematical difficulty of finding multiple roots. If you need to find all the roots, you must call uniroot multiple times and each time specify an interval containing only one root.

c(uniroot(f, c(0,1))$root, uniroot(f, c(1,2))$root)

[1] 0.6293433 1.4478550

The package rootSolve contains a function uniroot.all which attempts to automate this procedure. It divides the search interval into many subintervals (the default is 100) and uses uniroot on those intervals which have a sign change.

uniroot.all(f, c(0,2))

[1] 0.6293342 1.4478556

However, this approach is not guaranteed to work unless the number of subintervals is sufficiently large (which depends on how close together or to the boundary the roots are). I was recently working on an algorithm that required finding multiple roots of many functions. I could not manually plot them or check to see if uniroot.all was using a sufficiently fine grid to capture all the roots. Luckily for me, I was able to transform the functions I was working with into trigonometric polynomials. And this brings us to another one of R’s great root finding functions, polyroot.

Note that the “poly” here does not stand for multiple, but rather for polynomial. Our running example happens to be a trigonometric polynomial (what a coincidence!). If we call polyroot with the coefficients (in order of increasing degree) it will return all roots, including complex valued ones.

coefs <- c(1/2, -5, 8, -4, 1)
roots <- polyroot(coefs)
roots

[1] 0.1226314+0.000000i 0.8084142-0.000000i
    1.5344772-1.639788i 1.5344772+1.639788i

Since the polynomial had degree 4, there are 4 roots. The first two are real. The complex roots correspond to roots of the polynomial which are not achieved by the trigonometric version because . As long as we carefully select the correct polynomial roots and transform back to the trigonometric scale, we’ll get all the roots of the original function:

acos(Re(roots[2:1]))

[1] 0.6293434 1.4478554

Another word of caution regarding substitutions: they might be numerically unstable. Careful readers may have noticed the above roots have some different digits depending on which function computed them. This error can be much larger if the substitution rule has a large derivative near one of the zeroes. This was the case for the functions that came up in my research, and I settled on the following scheme which was much more stable and accurate.

1. Use polyroot to find zeroes of transformed function.
2. Eliminate infeasible roots that fall outside the domain/range of the substitution rule.
3. Transform the remaining ones back to the original domain.
4. Partition domain into intervals each containing one of the transformed roots.
5. Call uniroot on each interval in the above partition.

The main limitation of this approach is that not all functions can be made into polynomials. Some functions, like sums of powers of for example, can be transformed by a trigonometric substitution into a complex polynomial using Euler’s famous formula. But there are still many more functions which cannot be transformed into polynomials by any stretch. And for these we are out of luck. We can only use something like uniroot.all and hope it works.

If you’re aware of other options I haven’t mentioned here, please let me know by commenting or contacting me directly!