You will not find a better root finder in the known universe. Jack knows what he's talking about; he's tried them all.
In recent months, I worked on ways to find the minimum of some function:
I lamented more than once that life would be much simpler if we only needed to find the root of a similar function. By “root,” I mean the place where:
I also hinted that I had a killer solution to the root-finding problem, and promised to show it to you some day. That day has arrived.
I first wrote about this problem in 1994.[1,2,3] Those of you (and there are some) who clip out and save all my columns will enjoy comparing what I do here to the previous pass.
Finding a root is easier than finding a minimum. If you begin with two values of x that straddle the root-that is, one of them has a positive value and the other a negative one-you just know that somewhere in between them, the function has to cross zero. It will always be so for any function that's continuous and has a continuous first derivative. The algorithm I'm going to show you is absolutely guaranteed to find that minimum within some epsilon error, as long as you don't make epsilon so small as to create numerical accuracy problems.
You don't find many iterative methods that guarantee a solution, so this one is pure gold.
In the case of a function minimum, the curve typically looks like a parabola near the minimum. As we crank up the gain-that is, narrow down the range of the x-values, that curve tends to flatten out and show less variation in the y-value. In the limit, it gets so flat that we can't tell the value of one f(x) from another, and we can lose the minimum entirely.
This doesn't happen with a root finder. For a root finder, the curve becomes more and more like a straight line crossing the x axis, for which the solution is trivial. That's the big difference between the two problems, and the reason I found myself wishing the one were more like the other.
Let's state the problem more rigorously. Let there be some function f(x), and two values, x0 and x1 , such that the function changes sign somewhere between them:
Then, seek the value x for which f(x) = 0.
When I last looked at this problem, I took three months to cover it. That's because I looked at a number of possible solutions, including bisection, the secant method, reguli falsi, and so on. We don't have to do that this time, because I only want to show you the ultimate solution.
Figure 1 Parabolic interpolation
A first-order method is based on the idea of fitting a straight line between two points lying on the function f(x). Not surprisingly, a second-order method is based on fitting a quadratic through three points. Consider the curve of Figure 1. We can fit a quadratic through the three points shown, to get the parabola described by the dashed curve. The general form of such a quadratic is:
(4)but by playing around with the coefficients, we can eliminate the first-order term to get:
where xv and yv are the values at the vertex of the parabola.
This equation describes a parabola opening up (for positive A) or down (for negative A). By looking at the values of the function at x = x0 , x1 , and x2 , we could conceivably determine the constants xv, yv, and A. The root is then given by setting y = 0, to get:
This is Mueller's method. Now, looking at Equation 6, tell me, what's wrong with this picture? That's right, we have two roots (potentially complex), not just one. The function is double-valued when solved for x as a function of y.
While it's true that, for most cases, the roots are liable to be far enough apart to tell us which one we should use, we can't guarantee it. So how do we decide which of the two roots is the right one? And what if the roots are complex? In that case, we get no useful solution. While we can conceivably test the results to be sure Equation 6 leads to the proper root, wouldn't it be nice if our algorithm could assure us one real root in all cases?
This is actually an “inverse quadratic interpolation.” Instead of a parabola opening along the y-axis, let's use one that opens along the x-axis; that is, it will be a quadratic in y rather than x. A function like this will always be single-valued, since each value of y gives us only one value of x. Since y can take on any value, we needn't worry about complex or imaginary roots. An inverse quadratic fit is shown in Figure 2.
Figure 2 Inverse parabolic interpolation
I don't think there's much to recommend one or the other method, in terms of speed of convergence. I'm sure there are functions for which an upward-opening parabola fits better. There will be other functions for which the sideways-opening parabola fits best. In the long run, any quadratic function is bound to be better than the straight line of a linear interpolation. My guess is that, over many different functions, the convergence of the two methods is going to be about the same. The difference is in the math used to do the curve fit.
Fitting the curve
When you're trying to fit a curve through known points, the ease of the process depends on the form of the function. A general form like Equation 4 looks simple, but when you substitute in the values at the known points, you're in for some tedious and error-prone math. Here's a trick that simplifies things quite a bit. Write the equation in a factored form:
You can see the advantage to this form when you evaluate it at the known points. At point p0 , we have:
so the equation simply gives us:0 = A
which, of course, forces A to be zero.
The next coefficient is equally easy, since at point p1 the equation becomes:
(This is the inverse of the slope between points p0 and p1 .)
As you can see, we've already determined two of the three coefficients, with very little work. We get the last coefficient by evaluating at p2 :
Solving for C, we get:
Using Equation 9, this becomes:
We haven't yet placed any restrictions on the values of x0 , x1 , and x2 . Let's now assume that x1 is halfway between x0 and x2 , so that:
Then Equation 9 becomes:
And Equation 10 becomes:
Continuing with the solution for C gives:
We now have all three coefficients in terms of the fitted points. Remember this “factored polynomial” trick, which I also have been using in the minimization series. It's very useful, and works for polynomials of any order. As you can see, the computation of the coefficients gets progressively harder, but it will still be easier than with any other form.
Once we have the coefficients, we can find the desired root by setting y = 0 in Equation 7:
You might think we're home free at this point, since we need only solve Equation 14 for x. Unfortunately, a lot of snakes still lurk in the equations. Note that if any pair of y-values are equal, one of the terms in the denominators of Equations 12 or 15 will go to zero, and our equations blow up.
By requiring that points p0 and p2 straddle the root, we can be sure that:
(15) so at least one term can't vanish. But we can't control the value of y1 , so we can't promise one of the other two terms won't vanish. What do we do? We simply refuse to do the interpolation. In finding y1 , we have, in effect, taken one step of a bisection method. If we find that we can't do the interpolation, all is not lost. We simply re-order the data points so that we still straddle the root, take another bisection step, and repeat the process. Eventually, we get a useable interpolant.
Figure 3 A bad interpolation
Figure 4 A legal, but bad interpolation
If the terms in Equations 12 and 13 don't blow up, we might still find that our interpolation formula doesn't help. Consider the curve of Figure 4. While our formula works, it predicts a root well outside the range x0 ..x2 . Since we know that the desired root is in this range, we can't accept the interpolated value.
We need a mathematical test to tell us when we can trust the root estimated by the inverse interpolation. The obvious approach is to look at the predicted value of x. If it's not inside the range of our search, we refuse to accept the step. Unfortunately, we can get into trouble even then. In Figure 5, we get a mathematically correct fit to the data, and the predicted root is in the range of interest. But the accuracy of the fit is terrible; the predicted root, near p2 , is far from the actual one near p0 . We'd have been better off just accepting a bisection step. What's more, if we carry the method to the next step, we'll find that the next predicted root will be out of range.
Figure 5 Acceptance criterion
Deriving a suitable acceptance criterion is tricky, and I spent a long time finding one. But the final criterion is quite simple, and easily evaluated in software. The derivation is also fascinating in its own right, so bear with me as we delve into some tricky math.
We start by assuming that the root lies between points p0 and p2 , and neither point is the root (if it is, we're finished). This assumption means that y0 and y2 have different signs, and so, as we said in Equation 3:
After bisection to get y1 , we can tell from its sign whether the root lies to its left or its right. This is the basis for the method of bisection. We don't know which of the two cases we'll get, but I wish we did. Looking at Equation 12, we can see that the relationship between the three points is not symmetric. Points p0 and p2 enter into the equation much differently, thanks to our factored form. To avoid having to deal with multiple cases, we'd like much to always know that the root lies between p0 and p1 . We can assure this by the expedient of relabeling the points, if necessary.
If the root lies between p0 and p1 , then y0 and y1 must also have opposite signs, and y1 and y2 must have similar ones. Thus we can be assured that:
Now we'll assert our acceptance criterion, that the predicted root must also lie between x0 and x1 . This is equivalent to the assertion: (18) From Equation 14:
(19) (Because of our relabeling, h may no longer be positive.)
Substituting for B from Equation 9:
The inequality in Equation 13 becomes:
Now, h2 and the denominator are perfect squares, so their signs won't affect the outcome. We know from Equation 16 that the sign of y0 , y1 must be negative. We're left with the condition:
We substitute the value for C from Equation 13. After some horrible algebra, we find:
To simplify, let:
Then we have:
When we multiply these equations together, we get perfect squares in the denominator, which won't affect the sign of the result. So our inequality becomes:
We can identify four cases, depending on the signs of u and v:
Case 1: u > 0, v > 0 Okay
Case 2: u > 0, v < 0="">
Case 3: u < 0,="" v=""> 0 Forbidden
Case 4: u < 0,="" v="">< 0="">
Pushing the envelope
We must somehow relate these cases to the corresponding values of the function at points p0 , p1 , and p2 . We can explore the boundaries between the cases by setting u and v, separately, equal to zero. For example, when u = 0, we have: (28)
This equation is quadratic in y0 and y2 , but only linear in y1 , so we can solve for that parameter. We get:
Here, let's introduce two new variables:
Due to our restrictions on the ranges (see Equations 15 and16), both b and g must be negative.
Using these definitions, we rewrite Equation 29 in the form:
Rearranging the terms gives:
This is a quadratic in g, which we can solve using the quadratic formula. Let:
Then the quadratic formula gives:
This curve is plotted as the lighter curve in Figure 5. As you can see, it has two lobes, enclosing the two areas marked Case 4. This curve is a hyperbola centered at (-2, -1), with a tilt of slope 1/2. The two asymptotes have slopes of 0 and 1, respectively. A zero slope means a horizontal line; you can see that two legs of the curve trend to that direction.I can mark the area bounded by this curve as Case 4 is because it is, or at least appears, wholly enclosed by the second hyperbola, coming up.We can get the other boundary by setting v = 0. From Equation 25, this gives:
This time, it's y0 that we can solve for:
Using Equations 30 as before, this becomes:
A bit of rearranging gives us:
The hyperbola has limiting values of b given by:
I'm not showing the precise form of the equation for the hyperbola because we don't really need it. We can simply solve Equation 37 for g, just as we did for Equation 32. The equation has the form of a quadratic, with:
The quadratic formula gives us:
This curve is plotted as the heavier line in Figure 5. Between the two curves, defined by u = 0 and v = 0, we've defined the boundaries of all four cases. You might be wondering what happened to Case 3, since it's not shown on the curve. As it turns out, this case is represented by a tiny sliver where the two curves just overlap near the point b = 0.85, g = 0.5. The region is too small to show up on this graph. Its precise location and shape turns out to be immaterial anyway, because, while we've shown the graph for both positive and negative values of b and g, the restrictions given by Equations 16 and 17 demand that both b and g be negative. Thus, we're really only concerned with the lower left-hand quadrant of Figure 5.
Only Cases 1 and 4 satisfy the inequality in Equation 27. Within the lower left quadrant of Figure 5, we're allowed to operate in the triangular-shaped part of Case 1, and the part of Case 4 curving in from the left. I'm going to assert now that only the Case 1 portion will do us any good. Consider what happens as we start getting closer to convergence. Presumably, our boundary points p0 through p2 will get closer together, and the segment of the function they define will become more nearly a straight line. When the curve is a straight line, the quadratic of Equation 7 degenerates as the quadratic coefficient C approaches zero. Let's see what must happen to b and g as this degeneration takes place. Setting C to zero in Equation 13 gives:
This, of course, is a straight line, which I've plotted as the dashed line in Figure 5. As we are iterating and getting closer to the true root, we can expect the values for each iteration to migrate towards this line; the actual point will depend on the slope of the function curve. In the lower-left quadrant, this line lies wholly in the Case 1 region. Thus, our point can indeed migrate to the straight line. The Case 4 region, on the other hand, is cut off from the solution line by the Case 2 region surrounding it. Because of this, starting a quadratic step from within Case 4 is a fruitless exercise; we'll end up having to switch back to bisection anyway. Case 4, in fact, corresponds to the case of an acceptable, but poor, fit similar to that of Figure 4.
What we've gained out of these equations and graphs is the conviction that only the small triangular region in the lower right corner of the lower left quadrant will be suitable for inverse quadratic interpolation. In all other regions, we should stick with bisection.
Armed with this conviction, we're finally prepared to define a criterion for accepting the interpolated root. Look again at the four cases; the only difference between Case 1 and Case 3 is the sign of u. If we had to worry about falling into Case 3, we'd have to test the signs of both u and v. But since Case 3 only lies in that tiny sliver in the upper right quadrant, we needn't worry about it. We can never fall into that region. This means that we can ignore the sign of u, and accept as our criterion the condition:
(42) The value of v is given by Equation 25. This leads us to a simple test for acceptable interpolation:
Having struggled with all those equations and graphs, you can now forget them. They served only to establish the condition in Equation 43.
Now we can give the complete algorithm for this method:
- Given two points x0 and x2 straddling a root, bisect the region and evaluate the function to get y1 .
- Test the sign of y1 y2 . If it is negative, exchange the points p0 and p2 to make it positive.
- Test the sign of v = y2 (y2 — y0 ) — 2y1 (y1 — y0 ). If the sign is negative, replace p2 with p1 and repeat. If the sign is positive, compute the coefficients:
- Compute the estimated root from:
and evaluate the new function value, y = f(x).
- Test the sign of y to determine which points to use for the next cycle. If yy0 < 0,="" replace="">2 with (x,y). Otherwise replace p0 with (x,y) and p1 with p1 .
- Repeat until f(x)
I had hoped to give you code this month, but from the length of this document already, I can see that this isn't going to happen. So I'll have to hold off until next month to show you the finished code, but trust me, it'll be worth the wait. You won't find a better root-finder in the known universe.
Jack W. Crenshaw is a senior software engineer at Spectrum-Astro. He is the author of Math Toolkit for Real-Time Programming, from CMP Books. He holds a PhD in physics from Auburn University. Jack enjoys contact and can be e-mailed at .
1. Crenshaw, Jack. “Roots of a Function,”ESP, August 1994, p. 13.
2. Crenshaw, Jack. “None Other than e,” ESP, September 1994, p. 13.
3. Crenshaw, Jack. “Roots by Iteration,” ESP, October 1994, p. 13.