A cubic method that is Newton's in disguise

Newton’s method squares its error because it uses the first derivative to build a straight line that touches the curve, and the line’s zero is the next guess. A double root halves the error showed what happens when that line is a poor model. The obvious improvement is a better model: something that matches the curve not only in value and slope but in curvature too, whose zero should then be closer still.

Edmond Halley published one in 1694, eight years after Newton’s Principia, in a paper offering “a new, exact and easy method of finding the roots of any equations generally”, and it is still taught. He derived it for the cube root and then for polynomials, by an argument about successive corrections rather than about curves; the interpretation as a fitted hyperbola came much later. Its step is

$$z - \frac{2\,p(z)\,p'(z)}{2\,p'(z)^2 - p(z)\,p''(z)},$$

and near a simple root it cubes the error rather than squaring it — three times as many correct digits after each step instead of twice. What makes it worth an essay is not the speed, which turns out to be a small advantage, but a fact that makes the speed obvious and says a good deal more: Halley’s method is Newton’s method, applied to a different function.

Newton on another function #

Take $g = p/\sqrt{p'}$. It has the same roots as $p$ wherever $p' \ne 0$, since dividing by a non-zero number does not move a zero. Newton’s method for $g$ takes the step $g/g'$, and differentiating,

$$g' = \frac{p'}{\sqrt{p'}} - \frac{p\,p''}{2\,p'\sqrt{p'}} = \frac{2p'^2 - p\,p''}{2\,p'\sqrt{p'}}, \qquad \frac{g}{g'} = \frac{2\,p\,p'}{2p'^2 - p\,p''}.$$

That is Halley’s step, exactly. The figure’s comparison checked it at four points of the plane, twice over — once against the formula above and once against a numerical derivative of $g$ — so the identity is not a quotation.

Why that particular function? Because at a simple root, $g$ is unusually straight. Its second derivative vanishes at the root: dividing by $\sqrt{p'}$ is precisely the correction that cancels the curvature $p$ has there. Newton’s method on a function whose tangent line is also its best quadratic fit gains an extra order, so Newton’s squaring of $g$’s error becomes a cubing of $p$’s. The third-order method is the second-order method, run on a function chosen to have no curvature at its root.

The same trick reappears throughout numerical analysis. Newton’s method on $p/p'$ removed the multiplicity of a root, as the essay on double roots showed; Newton’s method on $p/\sqrt{p'}$ removes its curvature. In each case the method is not changed; the function is.

Cubing, drawn #

The speed map for Halley’s method on $z^3 - 1$ is paler everywhere than Newton’s. The discs of starts that arrive in three steps or fewer are much larger — each one fills most of its basin’s central lobe — and the average number of steps falls from 6.93 to 3.91. On the real line, starting from 3 and aiming at $\sqrt 2$, Halley reaches the limit of double precision in four steps where Newton needs five; its errors run $1.6$, $0.14$, $0.00028$, $2.8 \times 10^{-12}$ — the exponent roughly tripling each time. For $x^2 - 2$ the step simplifies to $x(x^2 + 6)/(3x^2 + 2)$, so from 3 the guesses are $45/29 \approx 1.5517$ and then $1.41449$: a rational function of the previous guess, as Newton’s $(x + 2/x)/2$ is, with every step exact arithmetic on fractions if one wants it.

The basins themselves have changed shape. Newton’s three basins interleave in a fractal along three rays from the origin, with a chain of ever-smaller lobes between every pair; Halley’s boundaries are much closer to three straight rays, with only small lobes along them. But they are not straight, and they are still fractal: every point of the boundary between two basins has points of the third basin arbitrarily close by. The same argument that where Newton’s method goes instead made for Newton’s map applies to any rational map with three attracting fixed points, and Halley’s map is one. A faster method has calmer basins, and the same kind of boundary.

Walking out of Newton’s trap #

An area that never finishes found the classic failure of Newton’s method: for $z^3 - 2z + 2$, the Newton map sends 0 to 1 and 1 back to 0, the two-cycle attracts, and a whole region of starting points around them — with area, not a thin boundary — never converges. In the left panel it is the warn-coloured patches near the real axis, 1.2% of the square.

The right panel has none. Halley’s method on the same polynomial sends 0 to 1, but 1 to 1.5, not back to 0 — a single computation, since at 1 the polynomial is 1, its derivative 1 and its second derivative 6, and the Halley step is $1 - 2 \cdot 1 \cdot 1/(2 - 6) = 1.5$. The cycle that trapped Newton does not exist for Halley, and every sampled start converges.

The escape is not a general theorem. Halley’s map is a different rational map, with its own free critical points, and for other polynomials it can have attracting cycles of its own — the escape here is the fact that this polynomial’s trap was built for Newton’s map and does not transfer. What does transfer is the lesson of the trap: whether a root-finder can get stuck depends on where its map sends its own critical points, and changing the map changes where they go.

What faster buys in practice #

A root-finder in real use is rarely started blind. The practical recipe is to locate each root roughly first and then polish it: count and isolate the real roots with Sturm’s chain, or start from a covering set of points guaranteed to reach every root, and only then iterate. In the polishing phase every start is already inside a pale disc, and the difference between the methods is the difference between squaring and cubing a small error — five steps against four to full precision, which the efficiency calculation below shows is nearly a wash.

Where Halley’s method earns its place is before the polishing, in the wild first steps from a poor start. The comparison panels show that region: Newton’s filaments of slow and scattered starts are thicker and more numerous than Halley’s, and on the polynomial with a trap, some of Newton’s starts never arrive. A method whose first steps are tamer needs a less careful initial search. The speed-up that matters is not the order; it is the shape of the basins at a distance. The natural hybrid takes Halley’s steps while a start is far out, where their tameness pays, and Newton’s once it is inside a pale disc, where the cheaper step loses almost nothing — and the two speed maps side by side show where along the way the switch would be worth making.

The double root, a third at a time #

At a double root Halley’s method also loses its high order, because the function $p/\sqrt{p'}$ still has a root of multiplicity greater than one there: $p$ vanishes to order 2 and $\sqrt{p'}$ to order $\tfrac12$, so $g$ vanishes to order $\tfrac32$. Newton’s method on a root of order $\tfrac32$ multiplies the error by $1 - \tfrac23 = \tfrac13$ each step — by the same calculation that gave $\tfrac12$ for Newton on $p$. In general Halley’s method at a root of multiplicity $m$ converges linearly with factor $(m - 1)/(m + 1)$, against Newton’s $(m - 1)/m$.

So Halley is faster at a double root too, but only by a constant factor: reducing an error from 1 to $10^{-7}$ takes about 15 steps at factor $\tfrac13$ and 23 at factor $\tfrac12$. The speed map shows it as a paler version of the same picture — the double root’s basin is still uniformly slow, the average down from 18.80 to 11.61 steps. No method that evaluates $p$ alone can reach a double root quickly without knowing it is double, and none can locate it to more than about eight digits, because the rounding floor belongs to the polynomial.

A family with one member at each order #

Halley’s method is the second member of a family that Alston Householder organised in 1970, with Newton’s method as the first. The method of order $d + 1$ uses the derivatives of $1/p$ up to order $d$: Newton’s step is $1/p$ divided by its first derivative, rescaled; Halley’s uses the ratio of the first and second derivatives of $1/p$; the next member the second and third. Each step lands, near a simple root, $d + 1$ times as many digits deep as the last.

Why the reciprocal? Near a simple root, $1/p$ has a simple pole, and the derivatives of a function with a pole grow in a very regular way — the ratio of successive derivatives tends to a constant that encodes the distance to the pole. That is König’s theorem of 1884, and it is the same phenomenon how fast the staircase arrives meets in the ratio of successive errors: a sequence whose consecutive ratios settle reveals the size of the thing controlling it. Householder’s methods read off the distance to the nearest root from the ratios of derivatives of $1/p$, and a higher order means reading it from further along the sequence.

Every member of the family can also be written as Newton’s method on some modified function, as Halley’s was. The order of a method is then a statement about how flat that function is made at its root, and the cost of a method is the number of derivatives the flattening needs.

What the extra order costs #

Halley’s step needs $p''$ as well as $p$ and $p'$: three evaluations where Newton needs two. A fair comparison divides the gain by the cost. Newton multiplies the number of correct digits by 2 per two evaluations, which is a factor of $\sqrt2 \approx 1.414$ per evaluation; Halley multiplies it by 3 per three evaluations, $\sqrt[3]{3} \approx 1.442$ per evaluation. Per unit of work, Halley’s method is about two per cent faster. For a polynomial, whose derivatives are cheap by Horner’s scheme, it can be worth having; for a function whose second derivative is expensive, it usually is not.

That arithmetic is the reason higher-order methods are not simply better. Householder’s methods continue the sequence — order 4 with the third derivative, and so on — and each costs one more evaluation per step for one more order; the efficiency $\sqrt[d]{d}$ peaks at $d = 3$ and falls after. The practical gains of Halley’s method are the ones the pictures show rather than the ones its order suggests: larger regions of fast arrival, fewer wild first steps, and, for some polynomials, an escape from Newton’s traps.

Why the basins are calmer #

The larger pale discs have a geometric reason. Newton’s method models the curve by its tangent line, and the tangent’s zero can land far away wherever the slope is small — near any critical point of $p$, the step $p/p'$ is enormous and the next guess is flung across the plane. That flinging is what builds the fractal filaments: a region near a critical point is scattered over every basin. And the critical points are never far from the action — the roots of the slope stay inside the smallest convex region containing the roots — so every polynomial with several roots has places between them where the tangent line is nearly flat and Newton’s step is wild.

Halley’s method models the curve by a hyperbola — a function of the form $(z - a)/(bz + c)$ — chosen to match $p$ in value, slope and curvature, which is why it is also called the method of tangent hyperbolas. A hyperbola can bend back towards the axis where a line would run off to infinity, so near a critical point Halley’s step stays moderate: the denominator $2p'^2 - p\,p''$ does not vanish where $p'$ does, unless $p$ or $p''$ also vanishes there. Fewer starts are flung far, and the regions of fast arrival grow. Its map has critical points too, and near them the scattering returns — which is why the boundaries in the right-hand panels still carry small lobes.

What the side-by-side panels cannot show #

The panels are sampled grids, 130 starting points on a side, and a basin’s fractal boundary is finer than any grid; “misses on 0.0% of the square” means no sampled point failed, not that no point fails. For $z^3 - 1$ under either method the boundary itself, a set of no area, consists of points that never converge to any root, and the samples miss it with probability one.

The comparison is of steps, not of time. Halley’s 3.91 steps each cost more than Newton’s 6.93 steps; the efficiency calculation above is the fair comparison, and the maps are not.

The identity is checked at four points, not proved by the figure. The algebra above proves it; the four sample points confirm that the implementation of Halley’s step and the formula for $g/g'$ agree, which guards against a transcription error in the step rather than establishing the identity.

And one trap escaped is one polynomial. That Halley converges from everywhere sampled for $z^3 - 2z + 2$ says nothing about polynomials designed to trap Halley’s map, which exist by the same construction that produced this one for Newton’s.

Still open: a universal starting set for Halley #

For Newton’s method there is a finite set of starting points, depending only on the degree, from which every root of every polynomial of that degree is found; its size grows roughly like $d \log^2 d$. The proof uses specific properties of Newton maps — that every root’s basin reaches out to infinity through channels of controlled width — which are known because Newton maps of polynomials have been classified.

No such result is known for Halley’s method or for the other methods of Householder’s family. Whether their basins also reach infinity in controlled channels, and so whether a universal starting set exists and how large it must be, is open; the covering essay’s closing section raised the question, and it has not been settled for any method of order higher than two.

The difficulty is structural. Newton’s map for a polynomial of degree $d$ is a rational map of degree $d$ whose only poles are at the critical points of $p$, and its dynamics have been classified in detail; Halley’s map for the same polynomial has degree up to $2d - 1$, extra free critical points, and a parameter space — the analogue of the picture of the quadratic family’s parameters — that has not been mapped. Each extra free critical point is another place an attracting cycle might hide.

The method was never the point #

Halley’s method is faster than Newton’s, and the reason is not a cleverer iteration but a better-chosen function. Divide $p$ by $\sqrt{p'}$ and Newton’s method, unchanged, cubes the error; divide $p$ by $p'$ and it recovers from a multiple root. The basins change with the function, sometimes dramatically — a trap disappears — and the character of the boundary does not, because both methods iterate a rational map with several attracting fixed points, and such maps always have fractal boundaries between their basins.

The pictures agree with the algebra at every point where they can be compared: the calmer basins are Newton’s basins for a straighter function, the escape is Newton’s map for that function having no cycle where $p$’s map had one, and the double root’s factor of a third is Newton’s factor for a root of order one and a half.

When a method is improved, ask whether the method changed or the problem did. Here the problem did: Halley’s method is Newton’s method asked a slightly different question, and everything that distinguishes the two — the cubing, the calmer basins, the escape from the cycle, the factor of a third at a double root — follows from the question.