Root Finding and Fixed-Point Iteration

Given $f : \mathbb{R}^n \to \mathbb{R}^n$, find $\mathbf{x}^\star$ such that:

$$f(\mathbf{x}^\star) = \mathbf{0}$$

For $n=1$, this is finding the zeros of a scalar function. For $n > 1$, this is solving a system of nonlinear equations.

from scipy.optimize import fsolve, root, brentq

# Scalar: brentq (bracketing, guaranteed for continuous f)
x_star = brentq(f, a, b)  # f(a) and f(b) must have opposite signs

# System: fsolve (Levenberg-Marquardt)
x_star = fsolve(func, x0)

# System: root (multiple methods)
result = root(func, x0, method='hybr')

For a continuous $f : \mathbb{R} \to \mathbb{R}$ with $f(a) \cdot f(b) < 0$ (sign change), the bisection method repeatedly halves the interval:

1. Compute $c = (a + b) / 2$
2. If $f(a) \cdot f(c) < 0$, set $b = c$; else set $a = c$
3. Repeat until $|b - a| < \epsilon$

Convergence is linear with rate $1/2$: error halves each iteration. After $k$ iterations, $|x^\star - c_k| \le (b-a)/2^{k+1}$.

from scipy.optimize import brentq

# brentq combines bisection with inverse quadratic interpolation
x = brentq(f, a, b, xtol=1e-12)

brentq is guaranteed to converge for continuous functions and is the most robust scalar root-finder in SciPy.

For solving $f(x) = 0$, Newton's method linearizes $f$ at $x_k$:

$$x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$

For systems $\mathbf{f}(\mathbf{x}) = \mathbf{0}$:

$$\mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{J}(\mathbf{x}_k)^{-1} \mathbf{f}(\mathbf{x}_k)$$

where $\mathbf{J}$ is the Jacobian matrix $J_{ij} = \partial f_i / \partial x_j$.

Convergence is quadratic near the root (digits of accuracy double each iteration), but the method can diverge from poor starting points.

from scipy.optimize import newton

# Scalar with derivative
x = newton(f, x0, fprime=df)

# System: use root with method='hybr' (Powell's hybrid)
res = root(func, x0, jac=jacobian, method='hybr')

A fixed point of $g : \mathbb{R}^n \to \mathbb{R}^n$ is a point $\mathbf{x}^\star$ where $g(\mathbf{x}^\star) = \mathbf{x}^\star$.

The iteration $\mathbf{x}_{k+1} = g(\mathbf{x}_k)$ converges to $\mathbf{x}^\star$ if $g$ is a contraction: $\|g(\mathbf{x}) - g(\mathbf{y})\| \le L\|\mathbf{x} - \mathbf{y}\|$ with $L < 1$.

from scipy.optimize import fixed_point

# Solve x = cos(x)
x_star = fixed_point(np.cos, 0.5)
# x_star ≈ 0.7390851332 (Dottie number)

Root finding $f(\mathbf{x}) = \mathbf{0}$ is equivalent to fixed-point iteration with $g(\mathbf{x}) = \mathbf{x} - \alpha f(\mathbf{x})$ for suitable $\alpha$.

Anderson acceleration speeds up fixed-point iteration by combining the last $m$ iterates:

$$\mathbf{x}_{k+1} = (1-\beta)\sum_{i=0}^{m_k} \theta_i \mathbf{x}_{k-i} + \beta \sum_{i=0}^{m_k} \theta_i g(\mathbf{x}_{k-i})$$

where the weights $\theta_i$ are chosen to minimize the residual norm in a least-squares sense.

from scipy.optimize import root

# Anderson acceleration via root
def residual(x):
    return g(x) - x  # fixed point <=> root of residual

res = root(residual, x0, method='anderson')

Anderson acceleration is particularly effective for self-consistent field (SCF) iterations in electronic structure calculations and for ADMM in optimization.