optMethods.py

"""Optimization methods used by the path-comparison figures."""

import numpy as np


def _start(start_point):
    x = np.asarray(start_point, dtype=float).reshape(-1).copy()
    return x, [x.copy()]


def _track(points):
    return np.column_stack(points)


def _gradient(obj_fun, x):
    return np.asarray(obj_fun.g_x(x), dtype=float).reshape(-1)


def newton(start_point, obj_fun, modified=0, iterations=5, tolerance=1e-10):
    """Newton's method with the tutorial's optional fixed Hessian shift."""
    x, points = _start(start_point)
    identity = np.eye(x.size)
    for _ in range(iterations):
        gradient = _gradient(obj_fun, x)
        if np.linalg.norm(gradient) <= tolerance:
            break
        hessian = np.asarray(obj_fun.G_x(x), dtype=float) + modified * identity
        x = x + np.linalg.solve(hessian, -gradient)
        points.append(x.copy())
    return _track(points)


def BFGS(start_point, obj_fun, iteration=10, interpolation=0, tolerance=1e-10):
    """Inverse-BFGS with a guarded curvature update."""
    x, points = _start(start_point)
    inverse_hessian = np.eye(x.size)
    gradient = _gradient(obj_fun, x)

    for _ in range(iteration):
        if np.linalg.norm(gradient) <= tolerance:
            break
        direction = -inverse_hessian @ gradient
        if direction @ gradient >= 0:
            inverse_hessian = np.eye(x.size)
            direction = -gradient
        alpha = _line_search(obj_fun, x, direction, interpolation)
        step = alpha * direction
        next_x = x + step
        next_gradient = _gradient(obj_fun, next_x)
        y = next_gradient - gradient
        curvature = float(step @ y)
        if curvature > 1e-12 * np.linalg.norm(step) * np.linalg.norm(y):
            rho = 1.0 / curvature
            identity = np.eye(x.size)
            left = identity - rho * np.outer(step, y)
            inverse_hessian = (
                left @ inverse_hessian @ left.T + rho * np.outer(step, step)
            )
        x, gradient = next_x, next_gradient
        points.append(x.copy())
    return _track(points)


def quasi_newton(
    start_point, obj_fun, iteration=10, interpolation=0, tolerance=1e-10
):
    """Symmetric-rank-one inverse-Hessian quasi-Newton method."""
    x, points = _start(start_point)
    inverse_hessian = np.eye(x.size)
    gradient = _gradient(obj_fun, x)

    for _ in range(iteration):
        if np.linalg.norm(gradient) <= tolerance:
            break
        direction = -inverse_hessian @ gradient
        if direction @ gradient >= 0:
            inverse_hessian = np.eye(x.size)
            direction = -gradient
        alpha = _line_search(obj_fun, x, direction, interpolation)
        step = alpha * direction
        next_x = x + step
        next_gradient = _gradient(obj_fun, next_x)
        y = next_gradient - gradient
        residual = step - inverse_hessian @ y
        denominator = float(residual @ y)
        threshold = 1e-8 * np.linalg.norm(residual) * np.linalg.norm(y)
        if abs(denominator) > threshold:
            inverse_hessian += np.outer(residual, residual) / denominator
        x, gradient = next_x, next_gradient
        # Record the accepted point, not the previous iterate.
        points.append(x.copy())
    return _track(points)


def fletcher_reeves(
    start_point, obj_fun, iteration=10, alpha=0.1, tolerance=1e-10
):
    """Fletcher-Reeves nonlinear conjugate gradient with safe restarts."""
    x, points = _start(start_point)
    gradient = _gradient(obj_fun, x)
    direction = -gradient

    for _ in range(iteration):
        if np.linalg.norm(gradient) <= tolerance:
            break
        if direction @ gradient >= 0:
            direction = -gradient
        step_size = _armijo_line_search(obj_fun, x, direction, alpha0=alpha)
        next_x = x + step_size * direction
        next_gradient = _gradient(obj_fun, next_x)
        denominator = float(gradient @ gradient)
        beta = float(next_gradient @ next_gradient) / denominator if denominator else 0.0
        direction = -next_gradient + beta * direction
        x, gradient = next_x, next_gradient
        points.append(x.copy())
    return _track(points)


def _line_search(obj_fun, x, direction, interpolation):
    if interpolation > 0:
        return _backtracking_line_search(obj_fun, x, direction)
    return _armijo_line_search(obj_fun, x, direction)


def _backtracking_line_search(
    obj_fun, x, direction, c=1e-4, alpha0=1.0, contraction=0.5
):
    """Geometric Armijo backtracking."""
    if not 0 < c < 1 or not 0 < contraction < 1:
        raise ValueError("line-search constants must lie between zero and one")
    alpha = float(alpha0)
    if not np.isfinite(alpha) or alpha <= 0:
        raise ValueError("alpha0 must be positive and finite")
    value = obj_fun.f_x(x)
    slope = float(_gradient(obj_fun, x) @ direction)
    if slope >= 0:
        raise ValueError("line-search direction must be a descent direction")
    for _ in range(60):
        if obj_fun.f_x(x + alpha * direction) <= value + c * alpha * slope:
            # The old implementation contracted once more before returning.
            return alpha
        alpha *= contraction
    raise RuntimeError("backtracking line search failed")


def _armijo_line_search(obj_fun, x, direction, c=1e-4, alpha0=1.0):
    """Safeguarded quadratic-interpolation Armijo search."""
    if not 0 < c < 1:
        raise ValueError("c must lie between zero and one")
    alpha = float(alpha0)
    if not np.isfinite(alpha) or alpha <= 0:
        raise ValueError("alpha0 must be positive and finite")
    value = obj_fun.f_x(x)
    slope = float(_gradient(obj_fun, x) @ direction)
    if slope >= 0:
        raise ValueError("line-search direction must be a descent direction")

    for _ in range(60):
        trial_value = obj_fun.f_x(x + alpha * direction)
        if trial_value <= value + c * alpha * slope:
            return alpha
        denominator = 2.0 * (trial_value - value - alpha * slope)
        candidate = -slope * alpha**2 / denominator if denominator > 0 else alpha / 2
        # Safeguard interpolation so every rejected trial makes real progress.
        alpha = float(np.clip(candidate, 0.1 * alpha, 0.5 * alpha))
    raise RuntimeError("Armijo interpolation failed")