Aberth method

The Aberth method, or Aberth–Ehrlich method or Ehrlich–Aberth method, named after Oliver Aberth^[1] and Louis W. Ehrlich,^[2] is a root-finding algorithm developed in 1967 for simultaneous approximation of all the roots of a univariate polynomial.

This method converges cubically, an improvement over the Durand–Kerner method, another algorithm for approximating all roots at once, which converges quadratically.^[1][2] (However, both algorithms converge linearly at multiple zeros.^[3])

This method is used in MPSolve, which is the reference software for approximating all roots of a polynomial to an arbitrary precision.

Let p ( x ) = p n x n + p n − 1 x n − 1 + ⋯ + p 1 x + p 0 {\displaystyle p(x)=p_{n}x^{n}+p_{n-1}x^{n-1}+\cdots +p_{1}x+p_{0}} be a univariate polynomial of degree n {\displaystyle n} with real or complex coefficients. Then there exist complex numbers z 1 ∗ , z 2 ∗ , … , z n ∗ {\displaystyle z_{1}^{*},\,z_{2}^{*},\dots ,z_{n}^{*}} , the roots of p ( x ) {\displaystyle p(x)} , that give the factorization:

p ( x ) = p n ⋅ ( x − z 1 ∗ ) ⋅ ( x − z 2 ∗ ) ⋯ ( x − z n ∗ ) . {\displaystyle p(x)=p_{n}\cdot (x-z_{1}^{*})\cdot (x-z_{2}^{*})\cdots (x-z_{n}^{*}).}

Although those numbers are unknown, upper and lower bounds for their absolute values are computable from the coefficients of the polynomial. Now one can pick n {\displaystyle n} distinct numbers in the complex plane—randomly or evenly distributed—such that their absolute values are within the same bounds. (Also, if the zeros are symmetrical, the starting points must not be exactly symmetrical along the same axis, as this can prevent convergence.)^[1] A set of such numbers is called an initial approximation of the set of roots of p ( x ) {\displaystyle p(x)} . This approximation can be iteratively improved using the following procedure.

Let z 1 , … , z n ∈ C {\displaystyle z_{1},\dots ,z_{n}\in \mathbb {C} } be the current approximations of the zeros of p ( x ) {\displaystyle p(x)} . Then offset numbers w 1 , … , w n ∈ C {\displaystyle w_{1},\dots ,w_{n}\in \mathbb {C} } are computed as

w k = p ( z k ) p ′ ( z k ) 1 − p ( z k ) p ′ ( z k ) ⋅ ∑ j ≠ k 1 z k − z j , {\displaystyle w_{k}={\frac {\frac {p(z_{k})}{p'(z_{k})}}{1-{\frac {p(z_{k})}{p'(z_{k})}}\cdot \sum _{j\neq k}{\frac {1}{z_{k}-z_{j}}}}},}

where p ′ ( z k ) {\displaystyle p'(z_{k})} is the polynomial derivative of p {\displaystyle p} evaluated in the point z k {\displaystyle z_{k}} .

The next set of approximations of roots of p ( x ) {\displaystyle p(x)} is then z 1 − w 1 , … , z n − w n {\displaystyle z_{1}-w_{1},\dots ,z_{n}-w_{n}} . One can measure the quality of the current approximation by the values of the polynomial or by the size of the offsets.

Conceptually, this method uses an electrostatic analogy, modeling the approximated zeros as movable negative point charges, which converge toward the true zeros, represented by fixed positive point charges.^[1] A direct application of Newton's method to each approximated zero will often cause multiple starting points to incorrectly converge to the same root. The Aberth method avoids this by also modeling the repulsive effect the movable charges have on each other. In this way, when a movable charge has converged on a zero, their charges will cancel out, so that other movable charges are no longer attracted to that location, encouraging them to converge to other "unoccupied" zeros. (Stieltjes also modeled the positions of zeros of polynomials as solutions to electrostatic problems.)

Inside the formula of the Aberth method one can find elements of Newton's method and the Durand–Kerner method. Details for an efficient implementation, esp. on the choice of good initial approximations, can be found in Bini (1996).^[3]

The updates of the roots may be executed as a simultaneous Jacobi-like iteration where first all new approximations are computed from the old approximations or as a sequential Gauss–Seidel-like iteration that uses each new approximation from the time it is computed.

A very similar method is the Newton-Maehly method. It computes the zeros one after another, but instead of an explicit deflation it divides by the already acquired linear factors on the fly. The Aberth method is like the Newton-Maehly method for computing the last root while pretending you have already found the other ones.^[4]

The iteration formula is the univariate Newton iteration for the function

F ( x ) = p ( x ) ∏ j = 1 ; j ≠ k n ( x − z j ) {\displaystyle F(x)={\frac {p(x)}{\prod _{j=1;\,j\neq k}^{n}(x-z_{j})}}}

If the values z 1 , … , z n {\displaystyle z_{1},\dots ,z_{n}} are already close to the roots of p ( x ) {\displaystyle p(x)} , then the rational function F ( x ) {\displaystyle F(x)} is almost linear with a dominant root close to z k {\displaystyle z_{k}} and poles at z 1 , … , z k − 1 , z k + 1 , … , z n {\displaystyle z_{1},\dots ,z_{k-1},z_{k+1},\dots ,z_{n}} that direct the Newton iteration away from the roots of p(x) that are close to them. That is, the corresponding basins of attraction get rather small, while the root close to z k {\displaystyle z_{k}} has a wide region of attraction.

The Newton step F ( x ) F ′ ( x ) {\displaystyle {\tfrac {F(x)}{F'(x)}}} in the univariate case is the reciprocal value to the logarithmic derivative

F ′ ( x ) F ( x ) = d d x ln ⁡ | F ( x ) | = d d x ( ln ⁡ | p ( x ) | − ∑ j = 1 ; j ≠ k n ln ⁡ | x − z j | ) = p ′ ( x ) p ( x ) − ∑ j = 1 ; j ≠ k n 1 x − z j {\displaystyle {\begin{aligned}{\frac {F'(x)}{F(x)}}&={\frac {d}{dx}}\ln |F(x)|\\&={\frac {d}{dx}}{\big (}\ln |p(x)|-\sum _{j=1;\,j\neq k}^{n}\ln |x-z_{j}|{\big )}\\&={\frac {p'(x)}{p(x)}}-\sum _{j=1;\,j\neq k}^{n}{\frac {1}{x-z_{j}}}\end{aligned}}}

Thus, the new approximation is computed as

z k ′ = z k − F ( z k ) F ′ ( z k ) = z k − 1 p ′ ( z k ) p ( z k ) − ∑ j = 1 ; j ≠ k n 1 z k − z j , {\displaystyle z_{k}'=z_{k}-{\frac {F(z_{k})}{F'(z_{k})}}=z_{k}-{\frac {1}{{\frac {p'(z_{k})}{p(z_{k})}}-\sum _{j=1;\,j\neq k}^{n}{\frac {1}{z_{k}-z_{j}}}}}\,,}

which is the update formula of the Aberth–Ehrlich method.

• MPSolve A package for numerical computation of polynomial roots. Free usage for scientific purpose.