Householder's methods in Python
The iteration x <- x + d (1/f)^(d-1)/(1/f)^(d), which uses f, f', ..., f^(d) and has order d+1 at a simple root. d = 1 is Newton's method and d = 2 is Halley's.
- File
householder_methods.py, 124 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/householder-methods/. The authors are listed in prov/M9RH7W.jsonl.
"""Householder's method of order d + 1: x <- x + d (1/f)^(d-1)(x) / (1/f)^(d)(x).
The derivatives of 1/f are never formed directly. With the Taylor
coefficients a_j = f^(j)(x) / j! the scaled coefficients
beta_n = f(x)^(n+1) (1/f)^(n)(x) / n! satisfy
beta_0 = 1, beta_n = -sum_{j=1..n} a_j a_0^(j-1) beta_(n-j),
which has no division, and the step is a_0 beta_(d-1) / beta_d.
d = 1 is Newton's method and d = 2 is Halley's method.
The code works with float, and also with decimal.Decimal or
fractions.Fraction when derivs returns such values.
"""
import math
from dataclasses import dataclass, field
from typing import Callable, List, Sequence
@dataclass
class HouseholderResult:
root: object
iterations: int
evaluations: int # calls of derivs; each call gives f, f', ..., f^(d)
status: str # "converged", "exact", "stalled", "zero-denominator", "max-iter"
iterates: List[object] = field(default_factory=list)
@property
def converged(self) -> bool:
return self.status in ("converged", "exact")
def _finite(v) -> bool:
return not isinstance(v, float) or math.isfinite(v)
def householder_step(values: Sequence, d: int):
"""Return (a_0, beta_(d-1), beta_d) from f, f', ..., f^(d) at one point."""
a = []
fact = 1
for j in range(d + 1):
if j > 0:
fact *= j
a.append(values[j] / fact)
beta = [1]
for n in range(1, d + 1):
acc = 0
p = 1 # a_0^(j-1)
for j in range(1, n + 1):
acc += a[j] * p * beta[n - j]
p *= a[0]
beta.append(-acc)
return a[0], beta[d - 1], beta[d]
def householder(derivs: Callable[[object], Sequence], x0, d: int = 2,
tol=1e-14, max_iter: int = 100) -> HouseholderResult:
"""Run Householder's method of order d + 1 from x0.
derivs(x) must return at least d + 1 values f(x), f'(x), ..., f^(d)(x).
Stops with "converged" when the step is at most tol (1 + |x|), with
"exact" when f(x) = 0, with "stalled" when the step is 0 although
f(x) != 0, with "zero-denominator" when beta_d = 0, and with
"max-iter" after max_iter steps. Raises ArithmeticError when a value
or the step is not finite.
"""
if d < 1:
raise ValueError("d must be at least 1")
x = x0
iterates = [x]
evals = 0
for k in range(1, max_iter + 1): # @line L1
values = derivs(x) # @line L2
evals += 1
if len(values) < d + 1:
raise ValueError("derivs must return f, f', ..., f^(d)")
if not all(_finite(v) for v in values[:d + 1]):
raise ArithmeticError("non-finite value at x = %r" % (x,))
if values[0] == 0: # @line L3
return HouseholderResult(x, k - 1, evals, "exact", iterates)
a0, b_prev, b_d = householder_step(values, d) # @line L4
if b_d == 0: # @line L5
return HouseholderResult(x, k - 1, evals, "zero-denominator", iterates)
s = a0 * b_prev / b_d # @line L6
if not _finite(s):
raise ArithmeticError("non-finite step at x = %r" % (x,))
if s == 0:
return HouseholderResult(x, k - 1, evals, "stalled", iterates)
x = x + s # @line L7
iterates.append(x)
if abs(s) <= tol * (1 + abs(x)): # @line L8
return HouseholderResult(x, k, evals, "converged", iterates)
return HouseholderResult(x, max_iter, evals, "max-iter", iterates) # @line L9
def poly_derivs(coeffs: Sequence, d: int) -> Callable[[object], List]:
"""Return derivs(x) = [p(x), p'(x), ..., p^(d)(x)] for p with coefficients
coeffs, highest degree first, by repeated Horner division."""
cs = list(coeffs)
def derivs(x):
b = list(cs)
out = []
fact = 1
for j in range(d + 1):
if j > 0:
fact *= j
if len(b) == 0:
out.append(0 * x)
continue
acc = 0 * x
nb = []
for c in b:
acc = acc * x + c
nb.append(acc)
out.append(fact * nb[-1]) # j! times the j-th Taylor coefficient
b = nb[:-1]
return out
return derivs
test_householder_methods.py 153 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/householder-methods/. The authors are listed in prov/M9RH7W.jsonl.
import math
import random
import unittest
from decimal import Decimal, localcontext
from fractions import Fraction
from householder_methods import householder, householder_step, poly_derivs
CUBIC = [1, 0, -2, -5] # x^3 - 2x - 5, root 2.0945514815423265...
ROOT = 2.0945514815423265
class TestHouseholder(unittest.TestCase):
def test_converges_for_each_order(self):
for d in range(1, 7):
r = householder(poly_derivs(CUBIC, d), 2.0, d)
self.assertTrue(r.converged, (d, r.status))
self.assertAlmostEqual(r.root, ROOT, delta=1e-15)
self.assertEqual(r.evaluations, len(r.iterates) - 1 + (r.status == "exact"))
def test_step_counts_fall_with_order(self):
steps = [householder(poly_derivs(CUBIC, d), 2.0, d).iterations for d in range(1, 7)]
self.assertEqual(steps, [5, 4, 3, 3, 3, 3])
def test_d1_is_newton_and_d2_is_halley(self):
f = lambda x: x ** 3 - 2 * x - 5
df = lambda x: 3 * x * x - 2
d2f = lambda x: 6 * x
x = Fraction(2)
n = householder(poly_derivs([Fraction(c) for c in CUBIC], 1), x, 1, tol=0, max_iter=2)
h = householder(poly_derivs([Fraction(c) for c in CUBIC], 2), x, 2, tol=0, max_iter=2)
xn, xh = x, x
for k in range(2):
xn = xn - f(xn) / df(xn)
xh = xh - 2 * f(xh) * df(xh) / (2 * df(xh) ** 2 - f(xh) * d2f(xh))
self.assertEqual(n.iterates[k + 1], xn)
self.assertEqual(h.iterates[k + 1], xh)
def test_observed_order_is_d_plus_1(self):
with localcontext() as ctx:
ctx.prec = 400
cs = [Decimal(c) for c in CUBIC]
nd = poly_derivs(cs, 1)
r = Decimal(2)
for _ in range(12):
v = nd(r)
r -= v[0] / v[1]
for d in range(1, 6):
res = householder(poly_derivs(cs, d), Decimal(2), d, tol=Decimal(10) ** -380, max_iter=3)
e = [abs(x - r) for x in res.iterates[:4]]
q = (e[3].ln() - e[2].ln()) / (e[2].ln() - e[1].ln())
self.assertAlmostEqual(float(q), d + 1, delta=0.01)
def test_error_constant(self):
# e_{k+1} / e_k^(d+1) -> (-1)^d f'(x*) psi_d with psi = 1/phi, phi = f/(x - x*)
with localcontext() as ctx:
ctx.prec = 300
cs = [Decimal(c) for c in CUBIC]
nd = poly_derivs(cs, 1)
r = Decimal(2)
for _ in range(12):
v = nd(r)
r -= v[0] / v[1]
p0, p1 = 3 * r * r - 2, 3 * r
psi = [1 / p0]
for n in range(1, 4):
psi.append(-(p1 * psi[n - 1] + (psi[n - 2] if n >= 2 else 0)) / p0)
for d in (1, 2, 3):
res = householder(poly_derivs(cs, d), Decimal(2), d, tol=0, max_iter=6)
e = [x - r for x in res.iterates]
k = max(i for i in range(len(e) - 1) if abs(e[i + 1]) > Decimal(10) ** -250)
ratio = e[k + 1] / e[k] ** (d + 1)
self.assertAlmostEqual(float(ratio), float((-1) ** d * p0 * psi[d]), places=5)
def test_multiple_root_rate(self):
# (x - 1)^m (x + 2): rate (m - 1) / (m + d - 1)
for m in (2, 3):
cs = [1]
for _ in range(m):
cs = [a - b for a, b in zip(cs + [0], [0] + cs)]
cs = [a + 2 * b for a, b in zip(cs + [0], [0] + cs)]
with localcontext() as ctx:
ctx.prec = 60
for d in (1, 2, 3):
res = householder(poly_derivs([Decimal(c) for c in cs], d), Decimal(2), d, tol=0, max_iter=40)
e = [x - 1 for x in res.iterates]
self.assertAlmostEqual(float(e[-1] / e[-2]), (m - 1) / (m + d - 1), delta=1e-6)
def test_pure_power_exact_ratio(self):
for m in (2, 3, 4):
for d in (1, 2, 3):
res = householder(lambda x: [math.perm(m, j) * x ** (m - j) for j in range(d + 1)],
Fraction(1), d, tol=0, max_iter=4)
for a, b in zip(res.iterates, res.iterates[1:]):
self.assertEqual(b, a * Fraction(m - 1, m + d - 1))
def test_root_at_start(self):
r = householder(poly_derivs([1, 0, -4], 3), 2.0, 3)
self.assertEqual((r.status, r.root, r.iterations, r.evaluations), ("exact", 2.0, 0, 1))
def test_zero_derivative(self):
derivs = poly_derivs([1, 0, -2], 3)
self.assertEqual(householder(derivs, 0.0, 1).status, "zero-denominator") # Newton divides by f'
self.assertEqual(householder(derivs, 0.0, 2).status, "stalled") # Halley's step is 0
def test_no_real_root_hits_max_iter(self):
r = householder(poly_derivs([1, 0, 1], 2), 0.5, 2, max_iter=50)
self.assertEqual((r.status, r.iterations, r.evaluations), ("max-iter", 50, 50))
def test_non_finite_raises(self):
with self.assertRaises(ArithmeticError):
householder(lambda x: [float("nan"), 1.0, 0.0], 1.0, 2)
with self.assertRaises(ArithmeticError):
householder(lambda x: [math.inf, 1.0, 0.0], 1.0, 2)
with self.assertRaises(ValueError):
householder(poly_derivs(CUBIC, 1), 2.0, 2)
with self.assertRaises(ValueError):
householder(poly_derivs(CUBIC, 1), 2.0, 0)
def test_transcendental(self):
# f = exp(x) - 2: all derivatives equal exp(x)
for d in range(1, 6):
r = householder(lambda x: [math.exp(x) - 2] + [math.exp(x)] * d, 0.0, d)
self.assertTrue(r.converged)
self.assertAlmostEqual(r.root, math.log(2), delta=1e-15)
def test_step_helper(self):
a0, b1, b2 = householder_step([2.0, 3.0, 4.0], 2)
self.assertEqual((a0, b1, b2), (2.0, -3.0, 9.0 - 2.0 * 2.0))
def test_random_cubics(self):
rng = random.Random(20261107)
for _ in range(600):
roots = sorted(rng.uniform(-10, 10) for _ in range(3))
if min(roots[1] - roots[0], roots[2] - roots[1]) < 0.5:
continue
a, b, c = roots
cs = [1, -(a + b + c), a * b + b * c + a * c, -a * b * c]
i = rng.randrange(3)
gap = min(abs(roots[i] - roots[j]) for j in range(3) if j != i)
x0 = roots[i] + rng.uniform(-0.1, 0.1) * gap
d = rng.randint(1, 5)
r = householder(poly_derivs(cs, d), x0, d, max_iter=60)
self.assertTrue(r.converged, (roots, x0, d, r.status))
self.assertAlmostEqual(r.root, roots[i], delta=1e-9 * (1 + abs(roots[i])))
self.assertLessEqual(r.evaluations, 12)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build