Improvements of the Modified Anderson-Björck (modAB) Root-Finding Algorithm

1. Introduction

Iterative algorithms are generally fast but can easily diverge, if improper starting point is selected. Bracketing algorithms are more reliable in this respect. If the function is continuous in the enclosing interval [a, b] and has opposite signs at both ends, according to the Bolzano’s theorem, at least one root exists inside [a, b] and the bracketing method is guaranteed to find it with a certain precision. Bisection is the simplest method, and also worst-case optimal. The required number of iterations is constant for a certain precision and does not depend on the function properties. It can be precomputed by the following equation:

$$n={\log }_2{\displaystyle \frac{\left|x_2-x_1\right|}{2ϵ}}$$

(1)

However, for well-behaved functions, other algorithms like false position, Illinois [6] and Ridder’s [10] can provide significantly better performance, but for poor-behaved ones they can downgrade much below bisection. So, stability and performance are contradictory goals and achieving both is a challenging task. Some authors have noticed that this can be done by combining bisection with a faster algorithm, like false-position [13], secant or inverse-quadratic interpolation. Such are the Brent’s [11] and ITP [9] methods. Also, initial bisection steps can boost the performance of false-position and similar methods by bringing the starting point to more favorable conditions [1,5]. Although obtained by combining only linear algorithms, the resulting algorithm shows nonlinear (bilinear) behavior and better convergence than using both separately. In all cases, the selection of proper switching criteria is of key importance for applying each method at the most appropriate place and thus maximizing the overall performance.

2. The Original Modified Anderson-Björck (modAB) Algorithm

On that basis, a new hybrid algorithm was proposed in 2022 [1] that combined bisection with Anderson-Björck’s (AB) method. The switching criteria is based on measuring the linearity of the function. In this way, linear interpolation is applied when the deviation of the curve from a straight line becomes small enough according to the equation:

$$\left|y_3-y_m\right|<k\left(\left|y_3\right|+\left|y_m\right|\right),$$

(2)

where y₃ is the function value at midpoint, y_m is the chord ordinate at the same point and k is the limit deviation factor. After switching, the iterations continue as ordinary AB steps as follows:

• An intermediate point is calculated by linear interpolation:

$$x_3={\displaystyle \frac{x_{1}f\left(x_2\right)-x_{2}f\left(x_1\right)}{f\left(x_2\right)-f\left(x_1\right)}}$$

(3)

2.

The ordinate of the end that remained fixed at the previous iteration is reduced by a factor m. For the left end, it is calculated by the expression: y₁ = f (x₁) · m, where m = 1 − f (x₃) / f (x₂) if m > 0, otherwise: m = 0.5. For the right end, the equation is obtained by swapping indices 1 and 2, as follows: y₂ = f (x₂) · m, where m =1 − f (x₃) / f (x₁) if m > 0, otherwise: m = 0.5

An option to fall back to bisection is provided, if AB fails to converge within a certain number of iterations. The algorithm was included in Calcpad [2] and Root-Fortan [3] libraries. Although the proposed algorithm showed superior performance compared to others in the test suite and later in practice, we identified some drawbacks and opportunities for further improvements.

3. Drawbacks of the Original modAB Method

The previously proposed algorithm has several shortcomings that can be addressed:

• The deviation factor k in eq. (2) has a fixed value of k = 0.25 as a good average estimate that is experimentally tuned on a larger set of functions. However, this is not optimal because it does not account for the properties of each function in the selected interval.
• Although the possibility of redundant switch to AB was discussed in [1], nothing much has been done on the subject. Only a simple bisection fallback is included as a safety mechanism, if the algorithm fails to converge within the expected number of iterations. It is proposed empirically to be half of the required bisection steps. This is a simple and conservative approach, but in some cases, it may result in a “waste” of AB steps before it “realizes” that it is not working as required.
• The original algorithm uses the termination criteria of the false-position method:|x_n − x_n−1| < ε. Further tests have shown that it fails for certain functions that are flat around the root, but with large ordinates at endpoints ( > 1/ε). This causes the interpolation point to crawl very slowly along the abscissa, so the step can become lower than ε even if we are still too far from the root.
• The benchmarks in the original paper [1] were based on the number of iterations with a factor of two for Ridder’s method. However, the number of function calls is more closely related to the actual execution speed when the algorithm is implemented in software.
• In the original paper, we used 21 functions to benchmark the method performance. This set is very limited, which can bring some distortions to the results. So, one may consider it insufficient to derive reliable conclusions. It may occur that some algorithms are particularly suitable for a larger part of the functions in the set showing a false advantage. That is why additional tests were performed in [5] with a set of 50 functions.

4.1. Adaptive Value of Deviation Factor k

For highly curved functions, like the one on Figure 1, symmetry plays a great role in the behavior of linear interpolation algorithms. If endpoint ordinates are highly asymmetric like on the left plot, there is a possibility for the root to occur at a greater distance. This will slow down performance since there is not only a longer path to travel, but the step size is also smaller. If the endpoint values are nearly symmetric, the behavior improves significantly, although both cases have the same deviation from the straight line.

$$11·x^{11}-1=0@x\in \left[0.1;1\right]\hspace{1em}\hspace{1em}\hspace{1em}11·x^{11}-5=0@x\in \left[0.1;1\right]$$

That is why it is reasonable to include symmetry into the method switching criterion and the factor k seems to be the right place. For unsymmetric endpoints, we can keep k lower, forcing the algorithm to make more bisection steps. For symmetric endpoints we can keep k larger and allow the algorithm to switch faster to AB since we expect it to be more favorable. However, this is highly dependent on the shape of the function. So, there is no sense in tuning it too much for specific types of functions, because it may not work for other types. Instead, we can just map the value of k against the endpoint ordinates by using a simple relationship that is fast to evaluate and does not bring significant overhead on each iteration.

For that purpose, we will introduce a new factor of symmetry r by using the equation:

$$r=1-\left|{\displaystyle \frac{y_m}{y_2-y_1}}\right|,$$

(4)

where y₁ = f (x₁), y₂ = f (x₂) and y_m = 0.5 (y₁ + y₂) is the average value. Thus, we have r = 1 for symmetric and r = 0.5 for unsymmetric endpoint ordinates. Symmetry also ensures that the distance from the interpolation point to the root cannot exceed 0.5 (x₂ − x₁). Otherwise, it can approach the whole interval length of x₂ − x₁. But how to account for it in the deviation factor k? If we look closely at the above functions, we can notice that the relationship between symmetry and convergence is not linear. Higher symmetry leads to both shorter path and longer step, which is rather quadratic than linear. The simplest quadratic expression we can use is k = r². Thus, we map k to the interval [0.25, 1], where 0.25 is our old value for k, but with allowance to grow up to 1 for fully symmetric endpoints.

4.2. Adaptive Bisection Fallback Criterion

In this paper, we propose a new version of the modAB algorithm aiming to be more robust, reliable and efficient than the first one. This is achieved by fixing the above issues and making it more adaptive to the properties of the function. Instead of using a fixed number of expected iterations, we can apply more precise criterion, by measuring the actual performance of AB algorithm after switching. One very simple and direct approach is to check whether AB has shrunk the interval more than bisection would do. Assume that the algorithm switches from bisection to AB at the interval [a₀, b₀] with initial width:

$$W_0=b_0-a_0$$

(5)

After performing t more steps, a pure bisection method would reduce the width to $W_0/2^t$. So, we will define a fallback threshold as follows:

$$T\left(t\right) = C{\displaystyle \frac{W_0}{2^t}}$$

(6)

where C is a safety factor. Let $W\left(t\right)=b_t-a_t$ denotes the actual interval width after t AB steps. Then, the adaptive bisection fallback criterion becomes:

$$W\left(t\right)>T\left(t\right)$$

(7)

Taking C = 2 ⁿ will allow for a delay of n AB steps before switching back to bisection. Numerical experiments have shown that without a safety factor for certain functions the algorithm starts to oscillate rapidly between both methods, causing more “damage” than improvement. The question is how to select a suitable value for C. One or two iterations proved to be insufficient to prevent oscillations. A reasonable value of n appeared to be about 3 to 4. Larger values caused the algorithm to waste too many iterations on AB. So, we can safely pick n = 4 and thence C = 2 ⁴ = 16. For efficiency, we will not apply equation (6) directly but instead will initialize the threshold at the moment of switching as $T_0 = C{W}_0$. Then we will update it on each AB iteration t as follows: $T_t=T_{t-1}/2$.

To demonstrate how the new method switching criteria work, we will take a nearly symmetrical function f (x) = x³ – 0.001 and solve the equation f (x) = 0 in the interval [[-10, 10]. The plot is displayed on Figure 2.

Both original and proposed algorithms switch to AB at the first iteration, finding itself in a position to deal with a highly nonlinear and unsymmetric problem. But after that, they show very different behavior as summarized in the following table:

This demonstrates that using fixed switching criteria is not optimal and can possibly slow down the algorithm.

Table 1. Comparison between the original and the proposed method switching criteria.

Algorithm	Iteration number when
	switch to AB	fallback to bisection	switch to AB	converge and quit
The original modAB 2022 fixed criteria	1	24	34	40
The proposed new adaptive criteria	1	7	15	24

Table 1. Comparison between the original and the proposed method switching criteria.

Algorithm	Iteration number when
	switch to AB	fallback to bisection	switch to AB	converge and quit
The original modAB 2022 fixed criteria	1	24	34	40
The proposed new adaptive criteria	1	7	15	24

4.3. New Termination Check

Instead of the false-position convergence check|x_n − x_n−1| < ε, we will use the more robust and reliable bisection check b_n − a_n < ε_abs + ε_rеl · |x_n|, where a_n and b_n are the endpoints of the interval at iteration n, and ε_abs and ε_rеl are the absolute and relative tolerances, respectively. This criterion is not applicable for the pure false-position method because of the fixed endpoint issue. But having AB corrections, the interval will eventually shrink from both sides. However, for some functions this may require 1 or 2 extra steps after the solution has actually converged, just to shrink the interval from the other side. But it is a good trade to provide robustness and reliability of the method.

4.4. The Improved Algorithm

After implementing the above improvements and refactoring the structure, we arrive at the following improved version of the modAB algorithm:

Algorithm 1: The improved modAB algorithm

Calculate points p1: y1 = f (x1) and p2: y2 = f (x2) initialize “method” to “bisection” and “side” to “none” for iter = 1 to maxIterations  if “method” is “bisection” then:   calculate p3 as the midpoint such that x3 = (x1 + x2) / 2 and y3 = f (x3)   calculate the average ordinate by the equation: ym = (y1 + y2) / 2   calculate r to equation (4) and k = r2           if the function is close enough to linear, satisfying that           |ym - y3| < k (|y3| + |ym|) then:    change the method to “Anderson-Björck”   end  else if method is “Anderson-Björck” then:   calculate p3 as the secant point such that   x3 = (x1*y2 − y1*x2) / (y2 − y1) and y3 = f (x3)   clamp p3 between p1 and p2 in case floating point  rounding errors shoots it outside the interval  end  if the convergence check:  y3 = 0 or x2 – x1 ≤ ε_abs + ε_rel *|x3| is satisfied then:  stop and return x3 as a result  if y1 and y3 have equal signs then:   if “method” is “Anderson-Björck” then: (apply AB correction)    if “side” is “right” then     calculate m = 1 − y3 / y1     if m ≤ 0 then: divide y2 by 2 else: multiply y2 by m    else:     change “side” to “right”    end   end   for both methods, shrink the interval from left by setting p1 = p3  else   if “method” is “Anderson-Björck” then:    if “side” is “left” then:     calculate m = 1 − y3 / y2     if m ≤ 0 then: divide y1 by 2 else: multiply y1 by m    else:     change “side” to “left”    end   end   for both methods, shrink the interval from right by setting p2 = p3  end  if AB fails to shrink the interval enough, then:   reset the method back to “bisection” end return error, if failed to converge for the expected number of iterations

The source code of the algorithm is listed in Appendix A. It was also translated into Julia, C++ and Python/Cython languages and included in Calcpad [2], SciML/NonlinearSolve.jl [4], JuliaMath/Roots.jl [21], and ROOT.CERN [22] libraries for scientific computations.