Newton's method for systems in Python
Solve F(x) = 0 in n unknowns by linearizing: solve J(x) s = -F(x) and set x to x + s. Quadratic convergence near a root with nonsingular Jacobian, no guarantee far from it.
- File
newton_method_for_systems.py, 129 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/newton-method-for-systems/. The authors are listed in prov/EB6BTR.jsonl.
"""Newton's method for F(x) = 0 with F from R^n to R^n: solve J(x) s = -F(x), x <- x + s."""
import math
from dataclasses import dataclass, field
from typing import Any, Callable, List, Optional, Sequence
@dataclass
class NewtonSystemResult:
"""The returned point, why the loop ended, and the work done."""
x: List[Any]
reason: str # "root", "converged", "singular", "failure" or "kmax reached"
iterations: int
f_evals: int
j_evals: int
steps: List[Any] = field(default_factory=list) # infinity norms of the steps
@property
def converged(self) -> bool:
return self.reason in ("root", "converged")
def _finite(v: Any) -> bool:
return math.isfinite(float(v))
def norm_inf(v: Sequence[Any]) -> Any:
return max((abs(t) for t in v), default=0)
def solve(a: Sequence[Sequence[Any]], b: Sequence[Any]) -> Optional[List[Any]]:
"""Solve a s = b by Gaussian elimination with partial pivoting.
Returns None when a pivot is exactly zero, that is when the matrix is
singular (exactly, for Fraction entries) or numerically singular to the
point of a zero pivot. Cost: about 2n^3/3 multiplications.
"""
n = len(b)
m = [list(row) + [b[i]] for i, row in enumerate(a)]
for c in range(n):
p = max(range(c, n), key=lambda r: abs(m[r][c]))
if m[p][c] == 0:
return None
m[c], m[p] = m[p], m[c]
piv = m[c][c]
for r in range(c + 1, n):
t = m[r][c] / piv
if t != 0:
row_r, row_c = m[r], m[c]
for j in range(c, n + 1):
row_r[j] -= t * row_c[j]
s = [0] * n
for i in range(n - 1, -1, -1):
acc = m[i][n]
for j in range(i + 1, n):
acc -= m[i][j] * s[j]
s[i] = acc / m[i][i]
return s
def newton_system(F: Callable[[List[Any]], Sequence[Any]],
J: Callable[[List[Any]], Sequence[Sequence[Any]]],
x0: Sequence[Any], tol_abs: float = 0.0,
tol_rel: float = 2.0 ** -50,
max_iter: int = 100) -> NewtonSystemResult:
"""Run Newton's method on F from x0, with the Jacobian J.
Works for float and fractions.Fraction entries. Precondition: tol_abs >= 0,
tol_rel >= 0, not both zero, max_iter >= 1. Convergence is guaranteed
only near a root where J is nonsingular and Lipschitz. Cost per pass: one
call of F, one of J and one linear solve of O(n^3) operations.
"""
x = list(x0) # @line L1
fe = je = 0
steps: List[Any] = []
for k in range(1, max_iter + 1): # @line L2
fx = list(F(x)) # @line L3
fe += 1
if all(t == 0 for t in fx): # @line L4
return NewtonSystemResult(x, "root", k, fe, je, steps)
if not all(_finite(t) for t in fx): # @line L5
return NewtonSystemResult(x, "failure", k, fe, je, steps)
a = J(x) # @line L6
je += 1
if not all(_finite(t) for row in a for t in row):
return NewtonSystemResult(x, "failure", k, fe, je, steps)
s = solve(a, [-t for t in fx]) # @line L7
if s is None or not all(_finite(t) for t in s):
return NewtonSystemResult(x, "singular", k, fe, je, steps)
x_new = [xi + si for xi, si in zip(x, s)] # @line L8
if not all(_finite(t) for t in x_new): # @line L9
return NewtonSystemResult(x, "failure", k, fe, je, steps)
ns = norm_inf(s)
steps.append(ns)
if ns <= tol_abs + tol_rel * norm_inf(x_new): # @line L10
return NewtonSystemResult(x_new, "converged", k, fe, je, steps)
x = x_new # @line L11
return NewtonSystemResult(x, "kmax reached", max_iter, fe, je, steps) # @line L12
def forward_difference_jacobian(F: Callable[[List[float]], Sequence[float]],
h_rel: float = 2.0 ** -26):
"""Return a function x -> forward-difference Jacobian of F at x.
Each call evaluates F n + 1 times; column j uses the step
h_j = h_rel * max(|x_j|, 1), the square root of the unit roundoff by
default. The counter attribute .calls counts the evaluations of F.
"""
def jac(x: List[float]) -> List[List[float]]:
n = len(x)
f0 = list(F(x))
jac.calls += 1
cols = []
for j in range(n):
h = h_rel * max(abs(x[j]), 1.0)
xp = list(x)
xp[j] = x[j] + h
h = xp[j] - x[j]
fj = list(F(xp))
jac.calls += 1
cols.append([(fj[i] - f0[i]) / h for i in range(len(f0))])
return [[cols[j][i] for j in range(n)] for i in range(len(f0))]
jac.calls = 0
return jac
test_newton_method_for_systems.py 157 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/newton-method-for-systems/. The authors are listed in prov/EB6BTR.jsonl.
import math
import random
from decimal import Decimal, getcontext
from fractions import Fraction
from newton_method_for_systems import (forward_difference_jacobian, newton_system,
norm_inf, solve)
def F(v):
return [v[0] ** 2 + v[1] ** 2 - 4, v[0] * v[1] - 1]
def J(v):
return [[2 * v[0], 2 * v[1]], [v[1], v[0]]]
def exact_root():
getcontext().prec = 60
s3 = Decimal(3).sqrt()
return (2 + s3).sqrt(), (2 - s3).sqrt()
def dec(q):
return Decimal(q.numerator) / Decimal(q.denominator)
def test_circle_hyperbola_example():
r = newton_system(F, J, [2.0, 0.0])
assert r.reason == "converged"
assert (r.iterations, r.f_evals, r.j_evals) == (6, 6, 6)
X, Y = exact_root()
assert r.x[0] == float(X) and abs(Decimal(r.x[1]) - Y) <= Decimal(2) ** -53
assert r.steps[0] == 0.5
def test_observed_order_two_in_exact_arithmetic():
X, Y = exact_root()
x = [Fraction(2), Fraction(0)]
errs = []
for _ in range(7):
errs.append(max(abs(dec(x[0]) - X), abs(dec(x[1]) - Y)))
s = solve(J(x), [-t for t in F(x)])
x = [x[0] + s[0], x[1] + s[1]]
ratios = [errs[k + 1] / errs[k] ** 2 for k in range(6)]
assert all(Decimal("0.2") < q < Decimal("0.75") for q in ratios)
assert abs(ratios[5] - Decimal(2).sqrt() / 2) < Decimal("1e-6")
# observed order log(e2/e1)/log(e1/e0) close to 2
q = (errs[6] / errs[5]).ln() / (errs[5] / errs[4]).ln()
assert abs(q - 2) < Decimal("0.01")
def test_exact_root_with_fractions():
# F(x, y) = (x + y - 3, x - y - 1) is linear: one step lands on (2, 1)
G = lambda v: [v[0] + v[1] - 3, v[0] - v[1] - 1]
JG = lambda v: [[1, 1], [1, -1]]
r = newton_system(G, JG, [Fraction(10), Fraction(-7)])
assert r.x == [2, 1] and r.reason == "root"
assert (r.iterations, r.f_evals, r.j_evals) == (2, 2, 1)
def test_start_at_root():
r = newton_system(F, J, [Fraction(2), Fraction(0)], max_iter=1)
assert r.reason == "kmax reached" and r.x == [2, Fraction(1, 2)] and r.f_evals == 1
G = lambda v: [v[0] - 1, v[1] - 2]
r = newton_system(G, lambda v: [[1, 0], [0, 1]], [1, 2])
assert r.reason == "root" and (r.f_evals, r.j_evals) == (1, 0)
def test_singular_jacobian():
# J(0, y) is singular for F(x, y) = (x^2 + y^2 - 4, x y - 1)
r = newton_system(F, J, [0.0, 0.0])
assert r.reason == "singular" and r.f_evals == 1 and r.j_evals == 1
r = newton_system(F, J, [Fraction(1), Fraction(1)])
assert r.reason == "singular"
def test_singular_root_is_linear_and_needs_absolute_tolerance():
G = lambda v: [v[0] ** 2, v[1]]
JG = lambda v: [[2 * v[0], 0.0], [0.0, 1.0]]
r = newton_system(G, JG, [1.0, 1.0], max_iter=100)
assert r.reason == "kmax reached" and r.x[0] == 2.0 ** -100
for a, b in zip(r.steps[1:], r.steps[2:]):
assert b == a / 2
r = newton_system(G, JG, [1.0, 1.0], tol_abs=1e-8, max_iter=100)
assert r.reason == "converged" and r.iterations == 27
def test_no_solution_fails_or_stops():
# x^2 + y^2 + 1 = 0, x - y = 0 has no real solution
G = lambda v: [v[0] ** 2 + v[1] ** 2 + 1, v[0] - v[1]]
JG = lambda v: [[2 * v[0], 2 * v[1]], [1.0, -1.0]]
r = newton_system(G, JG, [1.0, 1.0], max_iter=50)
assert not r.converged
assert r.f_evals <= 50
def test_nan_and_inf():
G = lambda v: [math.log(v[0]) if v[0] > 0 else float("nan"), v[1]]
JG = lambda v: [[1 / v[0], 0.0], [0.0, 1.0]]
r = newton_system(G, JG, [10.0, 1.0])
assert r.reason == "failure"
H = lambda v: [math.exp(v[0]) - 1, v[1]]
JH = lambda v: [[math.exp(v[0]), 0.0], [0.0, 1.0]]
r = newton_system(H, JH, [-800.0, 0.0])
assert r.reason in ("singular", "failure")
def test_freudenstein_roth_from_standard_start():
G = lambda v: [-13 + v[0] + ((5 - v[1]) * v[1] - 2) * v[1],
-29 + v[0] + ((v[1] + 1) * v[1] - 14) * v[1]]
JG = lambda v: [[1.0, 10 * v[1] - 3 * v[1] ** 2 - 2],
[1.0, 3 * v[1] ** 2 + 2 * v[1] - 14]]
r = newton_system(G, JG, [0.5, -2.0])
assert r.converged and abs(r.x[0] - 5) < 1e-12 and abs(r.x[1] - 4) < 1e-12
assert r.iterations == 44
def test_finite_difference_jacobian():
jac = forward_difference_jacobian(F)
r = newton_system(F, jac, [2.0, 0.0])
assert r.converged and norm_inf([r.x[0] - 1.9318516525781366,
r.x[1] - 0.5176380902050416]) <= 2.3e-16
assert jac.calls == 3 * r.j_evals
def test_solve_pivoting_exact():
a = [[Fraction(0), Fraction(1)], [Fraction(1), Fraction(1)]]
assert solve(a, [Fraction(2), Fraction(3)]) == [1, 2]
assert solve([[1, 2], [2, 4]], [1, 2]) is None
def test_random_well_conditioned_systems():
rng = random.Random(20261108)
for _ in range(500):
n = rng.randint(1, 4)
z = [rng.uniform(-2, 2) for _ in range(n)]
c = [rng.uniform(-0.3, 0.3) for _ in range(n)]
# F_i(x) = x_i - z_i + c_i (x_i - z_i)^2 + 0.1 * sum_j (x_j - z_j)^3 / n
def G(v, z=z, c=c, n=n):
d = [v[i] - z[i] for i in range(n)]
t = sum(dj ** 3 for dj in d) / (10 * n)
return [d[i] + c[i] * d[i] ** 2 + t for i in range(n)]
def JG(v, z=z, c=c, n=n):
d = [v[i] - z[i] for i in range(n)]
return [[(1 + 2 * c[i] * d[i] if i == j else 0) + 3 * d[j] ** 2 / (10 * n)
for j in range(n)] for i in range(n)]
x0 = [z[i] + rng.uniform(-0.5, 0.5) for i in range(n)]
r = newton_system(G, JG, x0)
assert r.converged
assert norm_inf([r.x[i] - z[i] for i in range(n)]) <= 1e-14
assert r.f_evals <= 10 and r.f_evals == r.iterations
What the next level needs
- the tests did not run in this build