specular.optimization.minimize

Copy

minimize(objective_function, initial_point, *, direction='speg', step_size='constant', gradient=None, line_search_gradient=None, h=None, max_iter=1000, tol=1e-06, record_history=True, step_options=None)

Combine a direction and a step size rule in the iteration \(x_{n+1}=x_n+\gamma_n d_n\). Also available as specular.minimize.

Parameters

Parameter	Description
objective_function	Callable returning a real scalar. Its value at the initial point must be finite.
initial_point	Finite real scalar or nonempty one-dimensional vector.
direction="speg"	"speg", its alias "specular_gradient", or a callable (n, x) returning a finite direction with the same shape as x.
step_size="constant"	Named step rule, positive finite scalar, or callable (n, x, d) returning a positive finite scalar.
gradient=None	Raw gradient callback gradient(x) for SPEG and the stopping test. If omitted, numerical specular differentiation uses the selected backend.
line_search_gradient=None	Raw gradient callback for a named line search. Overrides its default classical or specular derivative.
h=None	Positive finite scalar interval for numerical specular differentiation; omitted intervals are chosen automatically.
max_iter=1000	Nonnegative integer limiting accepted updates. Zero returns the initial point without a gradient check.
tol=1e-6	Finite, nonnegative threshold for the raw gradient norm.
record_history=True	Whether to record trajectory arrays. Must be a Boolean.
step_options=None	Mapping of named or numeric step rule options; see make_step_size.

Returns

An OptimizationResult. With history enabled, point and objective histories have iteration + 1 entries, and the step history has iteration entries.

Notes

Indices passed to rule callables start at one. A custom direction is not normalized automatically. A schedule accepting only n must first be adapted with make_step_size. For callable step rules, configure any gradient on the callable itself; step_options and line_search_gradient cannot be supplied to this solver for that callable.

The stopping test uses the raw gradient, including when the direction is custom. A custom direction's size alone is therefore not a convergence criterion. The solver also checks the gradient at the last permitted iterate when max_iter is positive.

The solver runs an eager Python loop with NumPy state. Scalar objectives receive scalar points; vector objectives receive NumPy arrays during iteration and backend-specific arrays during numerical differentiation. Selecting Numba or JAX changes the calculation backend without compiling the complete solver. See the backend guide.

Ordinary Armijo and Wolfe searches use classical centered differences by default. The specular_ variants instead use the raw specular gradient. line_search_gradient overrides either choice. A line search's own finite-difference interval can be configured with step_options={"h": ...} when it uses ordinary centered differences. Default specular searches share the solver's h; a supplied line_search_gradient controls its own differentiation. Specular substitutions do not automatically inherit classical Wolfe convergence guarantees.

A LineSearchError is recorded as a failed search in the result. A non-finite proposed point or objective, a zero direction, or a step lost to working precision also stops without replacing the last accepted point. Invalid arguments and callback shape errors raise exceptions.

Examples

Copy

import numpy as np

from specular.optimization import minimize


def quadratic(x):
    return np.dot(x, x)


def quadratic_gradient(x):
    return 2.0 * x


result = minimize(
    quadratic,
    initial_point=[1.0, -2.0],
    direction="speg",
    step_size="strong_Wolfe",
    gradient=quadratic_gradient,
    line_search_gradient=quadratic_gradient,
    max_iter=100,
)
print(result.solution, result.stop_reason)

Source code in specular/optimization/solver.py

Copy

def minimize(
    objective_function: Callable,
    initial_point: Any,
    *,
    direction: str | Callable = "speg",
    step_size: str | float | Callable = "constant",
    gradient: Callable | None = None,
    line_search_gradient: Callable | None = None,
    h: float | None = None,
    max_iter: int = 1000,
    tol: float = 1e-6,
    record_history: bool = True,
    step_options: Mapping[str, Any] | None = None,
) -> OptimizationResult:
    """Iterate ``x_next = x + step_size(n, x, d) * d`` with ``d=direction(n,x)``.

    ``direction`` selects SPEG or a callable accepting ``(n, x)``. A callable
    ``step_size`` accepts ``(n, x, d)``; wrap a schedule accepting only ``n``
    with :func:`make_step_size`. Iteration indices start at one.

    ``gradient`` supplies the raw specular gradient for SPEG and termination.
    When omitted, the selected calculation backend approximates it with mesh
    ``h`` (automatic by default). Scalar and nonempty 1D real states are
    supported; the optimizer itself is an eager NumPy/Python loop.

    Ordinary line searches use ordinary centered differences by default.
    The ``specular_`` rules use the same specular gradient as the direction.
    ``line_search_gradient`` overrides either choice. The classical Wolfe
    theory does not automatically apply to a specular substitution.

    A failed search or invalid proposed update leaves the last accepted point
    intact. Invalid arguments and callback shape errors raise exceptions.
    """
    if not callable(objective_function):
        raise TypeError("objective_function must be callable")
    if gradient is not None and not callable(gradient):
        raise TypeError("gradient must be callable")
    if line_search_gradient is not None and not callable(line_search_gradient):
        raise TypeError("line_search_gradient must be callable")
    if isinstance(max_iter, (bool, np.bool_)) or not isinstance(max_iter, (int, np.integer)):
        raise TypeError("max_iter must be a nonnegative integer")
    if max_iter < 0:
        raise ValueError("max_iter must be a nonnegative integer")
    tol = _scalar(tol, name="tol")
    if tol < 0:
        raise ValueError("tol must be nonnegative")
    if h is not None:
        h = _scalar(h, name="h")
        if h <= 0:
            raise ValueError("h must be positive")
    if not isinstance(record_history, bool):
        raise TypeError("record_history must be a bool")
    x = _as_point(initial_point, name="initial_point")
    shape = np.shape(x)
    options = dict(step_options) if step_options is not None else {}

    def value_at(point):
        # Callbacks receive their own point so they cannot mutate accepted data.
        value = objective_function(_as_point(point))
        return _scalar(value, name="objective value", finite=False)

    def gradient_at(point):
        point = _as_point(point)
        if gradient is not None:
            value = gradient(point)
        elif np.ndim(point) == 0:
            value = derivative(objective_function, point, h=h)
        else:
            # Pass the original callable to preserve JAX tracing / Numba compilation.
            value = specular_gradient_value(objective_function, point, h=h)
        result = _as_point(value, name="gradient")
        if np.shape(result) != shape:
            raise ValueError("gradient shape must match initial_point")
        return result

    direction_fn = make_direction(direction, gradient=gradient_at, h=h)
    builtin_direction = isinstance(direction, str)
    if callable(step_size):
        if options:
            raise ValueError("step_options belong to a named or numeric step_size")
        if line_search_gradient is not None:
            raise ValueError("configure the gradient on the callable step_size itself")
        step_fn = step_size
        builtin_step = False
    else:
        # Never pass SPEG's raw gradient to an ordinary Wolfe rule implicitly.
        search_gradient = line_search_gradient
        is_specular_search = (
            isinstance(step_size, str) and step_size.strip().lower().startswith("specular_")
        )
        if search_gradient is None and is_specular_search:
            search_gradient = gradient_at
        if any(key in options for key in ("f", "gradient")):
            raise ValueError("use objective_function and line_search_gradient explicitly")
        step_fn = make_step_size(step_size, f=value_at, gradient=search_gradient, **options)
        builtin_step = True

    start = time.perf_counter()
    fx = value_at(x)
    if not math.isfinite(fx):
        raise ValueError("objective value at initial_point must be finite")
    points = [np.asarray(x).copy()] if record_history else []
    values = [fx] if record_history else []
    steps: list[float] = []
    iteration = 0
    success = False
    reason = "max_iter reached"

    # The endpoint check also recognizes convergence on the last allowed update.
    for n in range(1, max_iter + 2):
        if max_iter == 0:
            break
        q = gradient_at(x)
        if _norm(q) <= tol:
            success, reason = True, "gradient norm below tolerance"
            break
        if iteration == max_iter:
            break
        if builtin_direction:
            d = direction_fn(n, x, gradient_value=q)
        else:
            d = direction_fn(n, x)
        if _norm(d) == 0.0:
            reason = "zero direction with gradient above tolerance"
            break
        try:
            if builtin_step:
                gamma = step_fn(n, x, d, f_value=fx)
            else:
                gamma = step_fn(n, _as_point(x), _as_point(d))
            gamma = _scalar(gamma, name="step size")
            if gamma <= 0:
                raise ValueError("step size must be positive")
        except LineSearchError as exc:
            reason = f"line search failed: {exc}"
            break
        trial = _advance(x, d, gamma)
        if not np.all(np.isfinite(trial)):
            reason = "non-finite proposed point"
            break
        if np.array_equal(trial, x):
            reason = "step does not change the point at working precision"
            break
        trial = _as_point(trial)
        trial_value = value_at(trial)
        if not math.isfinite(trial_value):
            reason = "non-finite objective at proposed point"
            break
        x, fx = trial, trial_value
        iteration += 1
        if record_history:
            points.append(np.asarray(x).copy())
            values.append(fx)
            steps.append(gamma)

    history = {
        "variables": np.asarray(points, dtype=float).reshape((-1,) + shape),
        "values": np.asarray(values, dtype=float),
        "step_sizes": np.asarray(steps, dtype=float),
    }
    return OptimizationResult(
        solution=_as_point(x),
        func_val=fx,
        iteration=iteration,
        runtime=time.perf_counter() - start,
        history=history,
        stop_reason=reason,
        success=success,
        method=direction if isinstance(direction, str) else "custom",
    )