Householder's methods

These are a class of root-finding methods; that is, for the numerical solution of a single nonlinear equation, developed by Alston Scott Householder in 1970. They may be considered a generalisation of the well known Newton-Raphson method (also known more simply as Newton's method) defined by

\[ x\leftarrow x-\frac{f(x)}{f'(x)}. \]

where the equation to be solved is \(f(x)=0\).

From a starting value \(x_0\) a sequence of iterates can be generated by

\[ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}. \]

As is well known, Newton's method exhibits quadratic convergence; that is, if the sequence of iterates converges to a root value \(r\), then the limit

\[ \lim_{n\to\infty}\frac{x_{n+1}-r}{(x_n-r)^2} \]

is finite. This means, in effect, that the number of correct decimal places doubles at each step. Householder's method for a rate of convergence \(d+1\) is defined by

\[ x\leftarrow x-d\frac{(1/f)^{(d-1)}(x)}{(1/f)^{(d)}(x)}. \]

We show how this definition can be rewritten in terms of ratios of derivatives, by using Python and its symbolic toolbox SymPy.

We start by defining some variables and functions.

1from sympy import *
2x = Symbol('x')
3f = Function('f')(x)

Now we can define the first Householder formula, with \(d=1\):

1d = 1
2H1 = x + d*diff(1/f,x,d-1)/diff(1/f,x,d)
3H1

\[ x-\frac{f(x)}{\frac{d}{dx}f(x)} \]

which is Newton's formula. Now for \(d=2\):

1d = 2
2H2 = x + d*diff(1/f,x,d-1)/diff(1/f,x,d)
3H2

\[ x - \frac{2 \frac{d}{d x} f{\left (x \right )}}{- \frac{d^{2}}{d x^{2}} f{\left (x \right )} + \frac{2 \left(\frac{d}{d x} f{\left (x \right )}\right)^{2}}{f{\left (x \right )}}} \]

This is a mighty messy formula, but it can be greatly simplified by using ratios of derivatives defined by

\[ r_k=\frac{f^{(d-1}(x)}{f^{(d)}(x)} \] This means that \[ r_1=\frac{f}{f'},\quad r_2=\frac{f'}{f^{\prime\prime}} \] To make the substitution into the current expression above, we can use the substitutions \[ f^{\prime\prime}=f'/r_2,\quad f'=f/r_1 \] to be done sequentially (first defining the new symbols)

1r_1,r_2,r_3 = symbols('r_1,r_2,r_3')
2H2r = H2s.subs([(Derivative(f,x,2), Derivative(f,x)/r_2), (Derivative(f,x), f/r_1)]).simplify()
3H2r

\[ -\frac{2r_1r_1}{r_1-2r_2} \] Dividing the top and bottom by \(2r_2\) produces the formulation \[ \frac{r_1}{1-\displaystyle{\frac{r_1}{2r_2}}} \] and so Householder's method for \(d=2\) is defined by the recurrence \[ x\leftarrow x-\frac{r_1}{1-\displaystyle{\frac{r_1}{2r_2}}}. \] This is known as Halley's method, after Edmond Halley, also known for his comet. This method has been called the most often rediscovered iteration formula in the literature.

It would exhibit cubic convergence, which means that the number of correct figures roughly triples at each step.

Apply the same sequence of steps for \(d=3\), and including the substitution \[ f^{\prime\prime\prime} = f^{\prime\prime}/r_3 \] produces the fourth order formula \[ x\leftarrow x-\frac{3 r_{1} r_{3} \left(2r_{2} - r_{1}\right)}{r_{1}^{2} - 6 r_{1} r_{3} + 6 r_{2} r_{3}} \]

A test

We'll use the equation \[ x^5+x-1=0 \] which has a root close to \(0.7\). First Newton's method, which is the Householder method of order \(d=1\), and we start by defining the symbol \(x\) and the function \(f\):

1x = Symbol('x')
2f = x**5+x-1

Next define the iteration of Newton's method, which can be turned into a function with the handy tool lambdify:

1nr = lambdify(x, x - f/diff(f,x))

Now, a few iterations, and print them as strings:

 1y = 0.7
 2ys = [y]
 3for i in range(10):
 4    y = N(nr(y),100)
 5    ys += [y]
 6
 7for i in ys:
 8    print(str(i))
 9
100.7
110.7599545557827765973613054484303575009107589721679687500000000000000000000000000000000000000000000000
120.7549197891599746887794253559985793967456078439525201893202319456623650882121929457935763902468565963
130.7548776691557956141971506438033504033307707534709697222674827264390889507161368160254597915269779252
140.7548776662466927739251146002523856449587324643131536407777773148939177229546284200355119465808326870
150.7548776662466927600495088963585290075677963335246916447723036615900830138144428153523526591809355834
160.7548776662466927600495088963585286918946066177727931439892839706462440390043279509776806970677946058
170.7548776662466927600495088963585286918946066177727931439892839706460806551280810907382270928422503037
180.7548776662466927600495088963585286918946066177727931439892839706460806551280810907382270928422503037
190.7548776662466927600495088963585286918946066177727931439892839706460806551280810907382270928422503037
200.7548776662466927600495088963585286918946066177727931439892839706460806551280810907382270928422503037

We can easily compute the number of correct decimal places each time by simply finding the first place in each string where it differs from the previous one:

1for i in range(1,7):
2    d = [ys[i][j] == ys[i+1][j] for j in range(102)]
3    print(d.index(False)-2)

\begin{array}{r} 2\cr 3\cr 8\cr 16\cr 32\cr 66 \end{array}

and we see a remarkable closeness with doubling of the number of correct values each iteration.

Now, the fourth order method, with \(d=3\):

1r1 = lambdify(x,g(x)/diff(g(x),x))
2r2 = lambdify(x,diff(g(x),x)/diff(g(x),x,2))
3r3 = lambdify(x,diff(g(x),x,2)/diff(g(x),x,3))
4h3 = lambdify(x,x-3*r1(x)*r3(x)*(2*r2(x)-r1(x))/(r1(x)**2-6*r1(x)*r3(x)+6*r2(x)*r3(x)))

Now we basically copy down the above commands, except that we'll use 1500 decimal places instead of 100:

1y = 0.7
2ys = [str(x)]
3for i in range(10):
4    y = N(h3(x),1500)
5    ys += [str(y)]
6
7for i in range(1,6):
8    d = [xs[i][j] == xs[i+1][j] for j in range(1502)]
9    print(d.index(False)-2)

\begin{array}{r} 4\\ 19\\ 76\\ 308\\ 1233 \end{array}

and we that the number of correct decimal places at each step is indeed increased by a factor very close to 4.