Trust region methods - Executable Engineering: Interactive Numerical Methods in Python

Trust region methods define a region (a ball with radius $\Delta_k$ ) around the current guess in which the actual function is trusted to behave like a local model function (typically a second-order Taylor expansion). They are designed to improve robustness for complex surfaces, including indefinite and ill-conditioned regions.

They are geometrically dual to line search methods:

Line search: Choose a search direction $\mathbf{p}_k$ first, then search along that 1D ray for step length $\alpha_k$ .
Trust region: Choose a maximum step radius $\Delta_k$ first, then solve for the optimal direction and step vector $\mathbf{s}_k$ satisfying $\|\mathbf{s}_k\| \le \Delta_k$ .

If the actual reduction $f(\mathbf{x}_k) - f(\mathbf{x}_k + \mathbf{s}_k)$ agrees well with the model prediction, the step is accepted and the trust radius $\Delta_{k+1}$ is expanded; if agreement is poor, the step is rejected and $\Delta_k$ is shrunk.

import numpy as np
import scipy.optimize as optimize
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

from executable_engineering import plot_optimization_path_3d

Comparison of optimization methods¶

The primary unconstrained optimization families available in numerical scientific packages like scipy.optimize are:

Category	Optimization Technique	SciPy Method Identifier
Line Search (Quasi-Newton)	Broyden-Fletcher-Goldfarb-Shanno (BFGS)	`'BFGS'`
	Limited-memory BFGS (bound constrained)	`'L-BFGS-B'`
	Newton-Conjugate-Gradient	`'Newton-CG'`
Trust Region	Trust-Region Nearly Exact	`'trust-exact'`
	Trust-Region Newton-CG	`'trust-ncg'`
	Trust-Region Truncated Generalized Lanczos	`'trust-krylov'`
Direct / Derivative-Free	Nelder-Mead Downhill Simplex	`'Nelder-Mead'`
	Powell’s Conjugate Direction Method	`'Powell'`
First-Order Gradient	Non-linear Conjugate Gradient	`'CG'`

Example: Comparing optimizers on Himmelblau’s function¶

Himmelblau’s function has four identical global minima where $f(x_1, x_2) = 0$ . Let’s test the 'trust-exact' trust-region solver starting from $\mathbf{x}_0 = (1.5, 2.0)$ :

(Try swapping 'trust-exact' for 'BFGS', 'Newton-CG', or 'Nelder-Mead' to observe how different algorithm families traverse the multiminimum landscape!)

def f(x):
    return (x[0]**2 + x[1] - 11)**2 + (x[0] + x[1]**2 - 7)**2

def grad_f(x):
    x1, x2 = x
    df_dx1 = 4 * (x1**2 + x2 - 11) * x1 + 2 * (x1 + x2**2 - 7)
    df_dx2 = 2 * (x1**2 + x2 - 11) + 4 * (x1 + x2**2 - 7) * x2
    return np.array([df_dx1, df_dx2])

def hess_f(x):
    x1, x2 = x
    d2f_dx1dx1 = 12 * x1**2 + 4 * x2 - 42
    d2f_dx1dx2 = 4 * x1 + 4 * x2
    d2f_dx2dx2 = 2 + 4 * x1 + 12 * x2**2 - 28
    return np.array([[d2f_dx1dx1, d2f_dx1dx2], [d2f_dx1dx2, d2f_dx2dx2]])

x0 = np.array([1.5, 2.0])
guesses = [x0.copy()]

# Options to explore: 'trust-exact', 'trust-ncg', 'BFGS', 'Newton-CG', 'Nelder-Mead', 'CG'
method = 'trust-exact'
result = optimize.minimize(
    f, x0,
    method=method,
    jac=grad_f,
    hess=hess_f,
    callback=lambda xk: guesses.append(xk.copy())
)

print(f"Method {method}: converged in {result.nit} iterations to {np.round(result.x, 5)}")
print(f"Minimum function value: {result.fun:.4e}")

plot_optimization_path_3d(f, guesses, x_range=(-5, 5), y_range=(-5, 5), title=f"Himmelblau Optimization Path ({method})")