Halley's Method: Third-Order Root Finding

• Halley’s method is a third-order root-finding algorithm that refines Newton’s approach by incorporating second-derivative information for cubic convergence.
• It extends to Banach spaces and generalized equations, using a predictor-corrector form and majorant functions for semilocal analysis.
• Practical implementations leverage its explicit formula in high-performance computing, automatic differentiation, and nonlinear system solvers.

Halley’s method is a classical third-order root-finding algorithm that augments Newton’s method by incorporating second-derivative information. In its standard scalar form for solving $f(x)=0$, the iteration is

$$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$

and, under the usual simple-root regularity assumptions, it is locally cubically convergent (Ling et al., 2012). Contemporary work treats the method not only in the classical scalar and Banach-space settings, but also for generalized equations, rational dynamical systems, floating-point kernels, large-scale nonlinear solvers, and stochastic spectral-chaos discretizations (Roubal et al., 24 Apr 2025).

1. Classical iteration and equivalent formulations

For a scalar nonlinear equation $f(x)=0$, Halley’s method can be viewed as a higher-order correction to Newton’s method. One useful representation is

$$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$

which makes explicit that the Newton step is multiplied by a curvature-dependent rational factor (Ling et al., 2012). In the notation

$q=\frac{f(x)f''(x)}{[f'(x)]^2},$

the same factor is $C(q)=1/(1-q/2)$; near a simple root, $q\approx 0$, so Halley reduces to Newton plus a higher-order correction (Martin, 2023).

In Banach spaces, Halley’s method admits a predictor-corrector form that is especially useful for extensions beyond single-valued equations. With $f:X\to Y$ twice continuously Fréchet differentiable, the predictor $u_{k+1}$ solves the Newton linearization

$f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$

and the corrector $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$0 solves

$$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$1

When $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$2, this reduces to the familiar scalar Halley formula, but the predictor-corrector form is the one that extends naturally to generalized equations (Roubal et al., 24 Apr 2025).

Halley’s method is also the $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$3 member of the Householder family

$$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$4

for which $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$5 is Newton’s method and $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$6 is Halley’s method (Li et al., 2018). This placement is operationally important because it situates Halley as the lowest-order method beyond Newton that uses second derivatives while still retaining a simple closed form.

2. Local order, exactness properties, and structural limitations

The standard local statement is cubic convergence: if $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$7 is a simple root and $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$8 is sufficiently smooth, then

$$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$9

whereas Newton’s method is only quadratic (Kouba, 2011). In the scalar Householder perspective, this is the $f(x)=0$0 order corresponding to $f(x)=0$1 (Tan et al., 28 Jan 2025). In Banach spaces, the same cubic regime appears when the relevant differentiability and invertibility assumptions hold (Ling et al., 2012).

Several papers emphasize exactness properties that distinguish Halley from Newton. One characterization is that Halley’s method is exact for Möbius functions, whereas Newton is exact for affine functions (Martin, 2023). A related characterization arises from the Schwarzian viewpoint: the Schwarzian-Newton method is exact for functions with constant Schwarzian derivative, and Halley appears as the limiting case when the Schwarzian derivative tends to zero; in that limit, the SNM iteration reduces to Halley’s formula (Segura, 2015). These formulations clarify why Halley is often more accurate than Newton for functions whose local geometry is closer to rational than affine.

The same literature also records a structural limitation. In the correction-factor form $f(x)=0$2, the factor is positive for $f(x)=0$3 and negative for $f(x)=0$4. The paper “Root-finding: from Newton to Halley and beyond” identifies $f(x)=0$5 as the regime where Halley’s step reverses direction relative to the Newton step and describes this as undesirable behavior (Martin, 2023). This does not contradict cubic local convergence; rather, it shows that local order alone does not control global step geometry.

A more specialized exact representation appears in real-line dynamics for $f(x)=0$6. In that case Halley’s iteration becomes

$f(x)=0$7

and with the substitution $f(x)=0$8 one gets

$f(x)=0$9

This makes the map a ternary Bernoulli shift in angle coordinates and exposes periodic, divergent, and chaotic regimes in exact arithmetic (Li et al., 2018).

3. Banach-space and generalized-equation extensions

A substantial modern generalization replaces the equation $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$0 by the generalized equation

$$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$1

where $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$2 is single-valued and twice continuously Fréchet differentiable, while $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$3 is set-valued with closed graph (Roubal et al., 24 Apr 2025). Classical root-finding is recovered by taking $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$4.

The Josephy–Halley method extends Halley’s predictor-corrector structure to this setting. Given $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$5, one first computes a predictor $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$6 from the partially linearized inclusion

$$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$7

and then a corrector $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$8 from

$$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}\cdot\frac{1}{1-\tfrac12 L_f(x_n)}, \qquad L_f(x):=\frac{f(x)f''(x)}{2[f'(x)]^2},$$9

If $q=\frac{f(x)f''(x)}{[f'(x)]^2},$0, these reduce exactly to the classical Halley equations (Roubal et al., 24 Apr 2025).

The central regularity notion is metric regularity of the frozen linearization

$q=\frac{f(x)f''(x)}{[f'(x)]^2},$1

which serves as the generalized-equation analogue of invertibility or a quantitative inverse mapping property. Under metric regularity at $q=\frac{f(x)f''(x)}{[f'(x)]^2},$2 and Hölder continuity of $q=\frac{f(x)f''(x)}{[f'(x)]^2},$3 of order $q=\frac{f(x)f''(x)}{[f'(x)]^2},$4, the method converges locally with order $q=\frac{f(x)f''(x)}{[f'(x)]^2},$5, and therefore cubically when $q=\frac{f(x)f''(x)}{[f'(x)]^2},$6 (Roubal et al., 24 Apr 2025).

The Banach-space operator version used in semilocal theory can also be written as

$q=\frac{f(x)f''(x)}{[f'(x)]^2},$7

for nonlinear operator equations $q=\frac{f(x)f''(x)}{[f'(x)]^2},$8 on Banach spaces (Ling et al., 2012). This form isolates the Halley correction operator and is the standard bridge between classical scalar formulas and abstract convergence theory.

4. Semilocal theory, majorants, and Kantorovich-type analysis

Semilocal analysis for Halley’s method proceeds through majorant functions that dominate the operator derivatives. For Banach-space equations $q=\frac{f(x)f''(x)}{[f'(x)]^2},$9, one introduces a scalar majorant $C(q)=1/(1-q/2)$0 satisfying conditions on $C(q)=1/(1-q/2)$1, $C(q)=1/(1-q/2)$2, $C(q)=1/(1-q/2)$3, and the first zero $C(q)=1/(1-q/2)$4, together with the majorant inequality

$C(q)=1/(1-q/2)$5

Halley’s iteration for $C(q)=1/(1-q/2)$6 is then compared to Halley’s iteration for $C(q)=1/(1-q/2)$7, producing a majorizing scalar sequence $C(q)=1/(1-q/2)$8 (Ling et al., 2012).

A notable feature of this analysis is that it drops the requirement that the majorant function have a second root. Under assumptions (A1)–(A3), the Halley sequence remains in the ball $C(q)=1/(1-q/2)$9, converges to a zero $q\approx 0$0, and does so with Q-cubic rate; the paper also gives an explicit cubic error constant involving $q\approx 0$1, $q\approx 0$2, and the left directional derivative $q\approx 0$3 (Ling et al., 2012).

For generalized equations, semilocal analysis takes a more explicitly Kantorovich-type form. The scalar majorant is

$q\approx 0$4

where $q\approx 0$5 is a metric regularity constant, $q\approx 0$6 bound $q\approx 0$7 and its Lipschitz constant, and $q\approx 0$8 encodes the initial residual (Roubal et al., 24 Apr 2025). Halley-type scalar sequences $q\approx 0$9 and $f:X\to Y$0 are then defined by applying a scalar predictor-corrector Halley iteration to $f:X\to Y$1, and their convergence controls the vector iterates.

Under the semilocal assumptions in that framework, the Josephy–Halley sequence stays in an explicit neighborhood $f:X\to Y$2 and converges $f:X\to Y$3-cubically: $f:X\to Y$4 for some $f:X\to Y$5 and $f:X\to Y$6 (Roubal et al., 24 Apr 2025). In the single-valued case this yields a semilocal theorem for classical Halley under metric regularity of $f:X\to Y$7, which the paper identifies with surjectivity plus a bounded right inverse.

5. Rational dynamics, basin geometry, and global behavior

When applied to a polynomial $f:X\to Y$8, Halley’s method defines a rational map on the Riemann sphere,

$f:X\to Y$9

and the resulting Fatou–Julia decomposition governs global convergence from complex initial data (Canela et al., 30 Jul 2025). In this formulation, simple roots are superattracting fixed points, while additional non-root fixed points can appear. A later 2025 study states more generally that, for a polynomial $u_{k+1}$0, a finite fixed point of $u_{k+1}$1 is either a root of $u_{k+1}$2 or a critical point of $u_{k+1}$3 that is not a root, and that $u_{k+1}$4 is repelling (Liu et al., 1 Aug 2025).

For the symmetric family

$u_{k+1}$5

the Halley map has explicit symmetry and degree $u_{k+1}$6, and the global basin geometry differs sharply from Newton’s method. The main theorem in “Boundedness and simple connectivity of the basins of attraction for some numerical methods” shows that the immediate basins $u_{k+1}$7 of the nonzero roots are unbounded and simply connected, that $u_{k+1}$8 is unbounded for $u_{k+1}$9, and that $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$0 is bounded for $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$1 (Canela et al., 30 Jul 2025). This furnishes explicit examples in which a root basin for Halley’s method is bounded, a phenomenon the paper contrasts with Newton’s method.

The later dynamics paper refines the picture by separating classes. For unicritical polynomials, certain symmetric cubics, and certain quartics with nontrivial symmetry, $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$2 is convergent, $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$3 is connected, the immediate basins of the roots are unbounded, and the symmetry group of $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$4 coincides with that of the polynomial (Liu et al., 1 Aug 2025). The same paper then extends the analysis to broader classes and again shows that the immediate basin corresponding to a root can be bounded. It also proves that Halley’s method is not generally convergent for cubic polynomials by exhibiting non-root superattracting 2-cycles for certain parameter values (Liu et al., 1 Aug 2025).

These results correct a common oversimplification: cubic local convergence does not imply Newton-like global basin geometry. A plausible implication is that the global usability of Halley’s method is strongly problem-dependent even when its local asymptotics are favorable.

6. Implementations, specialized variants, and applications

At the hardware level, Halley’s method has been used as a compensated refinement for reciprocal square root. For $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$5, a correctly rounded Newton step can be followed by a Halley compensation implemented with fused multiply-add: $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$6 where $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$7. The paper reports that Halley compensation requires one additional FMA and one additional multiplication relative to Newton compensation and experimentally always matched the correctly rounded reference in the reported tests (Borges, 2021).

For large-scale nonlinear systems, higher-order automatic differentiation changes the practical cost model. The paper “Scalable higher-order nonlinear solvers via higher-order automatic differentiation” identifies Halley’s method as the most valuable Householder variant, computes only higher-order directional derivatives via Taylor-mode AD, and derives a multivariate update that reuses one Jacobian factorization for two linear solves. In the reported experiments, Halley was about 40% faster than Newton at $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$8 for a dense Chandrasekhar $f(x_k)+f'(x_k)(u_{k+1}-x_k)=0,$9-function problem, about 25% faster at $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$00 ($$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$01) for a sparse Brusselator steady-state problem, and gave an average 5–10% speedup for the entire stiff ODE solve when used inside implicit time-steppers (Tan et al., 28 Jan 2025).

Several domain-specific formulations exploit the structure of the underlying differential equations. For solutions of $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$02, a modified Halley method freezes $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$03 at the initial point,

$$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$04

and the paper proves that this modified scheme retains third-order convergence while reducing per-iteration cost; it then develops algorithms for all nodes and weights of Gauss–Legendre and Gauss–Hermite quadratures (K et al., 18 Mar 2026). In fractional two-point boundary value problems with Robin boundary conditions, Halley is used as the shooting method for the scalar residual $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$05, coupled with a quadratic-interpolation high-order predictor-corrector method; the reported overall convergence rate is $$x_{n+1} = x_n-\frac{2 f(x_n) f'(x_n)}{2[f'(x_n)]^2-f(x_n)f''(x_n)},$$06 (Lee et al., 2020). In stochastic eigenproblems solved by a spectral-chaos expansion, the projected system is quadratic, and the authors argue that Halley’s cubic convergence is maximal because higher-order Householder methods cannot improve the rate when third derivatives vanish identically (Esquivel et al., 24 Jun 2026).

Across these settings, the recurring pattern is consistent: Halley’s method occupies the first genuinely higher-order position beyond Newton that remains algebraically explicit, structurally extensible, and competitive in implementation when second-derivative information is either analytically available or efficiently generated.