Broyden Quasi-Newton Examples — NEQBF and NEQBJ

These examples solve the Rosenbrock system

\[ \begin{align}\begin{aligned}f_1(x, y) = 1 - x = 0\\f_2(x, y) = 10\,(y - x^2) = 0\end{aligned}\end{align} \]

whose unique solution is \((x, y) = (1, 1)\), starting from \((-0.5,\; 0.5)\). Two Broyden quasi-Newton variants are compared:

NEQBF — nonlinearequations.solve_system_broyden_fd()

No Jacobian required. Uses scipy.optimize.root with method='broyden1': the approximate Jacobian is updated by rank-one Broyden corrections without any explicit derivative information.

NEQBJ — nonlinearequations.solve_system_broyden()

Analytic Jacobian supplied. Uses scipy.optimize.root with method='hybr' (Powell’s hybrid method) which exploits the exact Jacobian to achieve rapid convergence with far fewer function evaluations.

The analytic Jacobian of the Rosenbrock system is:

\[\begin{split}J = \begin{pmatrix} -1 & 0 \\ -20x & 10 \end{pmatrix}\end{split}\]

Residual norms are recorded at every function evaluation by wrapping the system function, then plotted to compare convergence rates.

Example Code

"""IMSL Broyden quasi-Newton examples: NEQBF and NEQBJ.

Solves the Rosenbrock system:

    f1(x, y) = 1 - x        → 0
    f2(x, y) = 10*(y - x²)  → 0

Exact solution: (x, y) = (1, 1).

Two Broyden variants are compared:

* :func:`nonlinearequations.solve_system_broyden_fd` (IMSL NEQBF) — no
  Jacobian required; scipy uses an internal quasi-Newton update.
* :func:`nonlinearequations.solve_system_broyden` (IMSL NEQBJ) — supplied
  analytic Jacobian, backed by Powell's hybrid method.

Residual norms are captured at every function evaluation by wrapping the
system function, then plotted as convergence curves.

Outputs:
- Solutions and iteration counts printed to stdout.
- SVG convergence-curve figure saved to ``test_output/example_imsl_broyden.svg``.
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict, List

import matplotlib.pyplot as plt
import numpy as np

from nonlinearequations import solve_system_broyden_fd, solve_system_broyden


# ---------------------------------------------------------------------------
# Rosenbrock system definition
# ---------------------------------------------------------------------------

def rosenbrock_system(v: np.ndarray) -> np.ndarray:
    """Evaluate the Rosenbrock-derived nonlinear system at *v*.

    Args:
        v (np.ndarray): Current point [x, y].

    Returns:
        np.ndarray: Residuals [1 - x, 10*(y - x²)].
    """
    x, y = v[0], v[1]
    return np.array([1.0 - x, 10.0 * (y - x**2)])


def rosenbrock_jacobian(v: np.ndarray) -> np.ndarray:
    """Analytic Jacobian of the Rosenbrock system.

    Args:
        v (np.ndarray): Current point [x, y].

    Returns:
        np.ndarray: 2×2 Jacobian df_i/dx_j.
    """
    x = v[0]
    return np.array([
        [-1.0,    0.0],   # d(f1)/dx, d(f1)/dy
        [-20.0 * x, 10.0],  # d(f2)/dx, d(f2)/dy
    ])


# ---------------------------------------------------------------------------
# Residual-tracking wrappers
# ---------------------------------------------------------------------------

def _make_tracked(f) -> tuple:
    """Wrap *f* to record residual norms at each call.

    Args:
        f (callable): System function to wrap.

    Returns:
        tuple: (wrapped_f, residuals_list) where residuals_list is filled
            in-place as the solver calls wrapped_f.
    """
    history: List[float] = []

    def wrapped(v: np.ndarray) -> np.ndarray:
        """Tracked system function.

        Args:
            v (np.ndarray): Current point.

        Returns:
            np.ndarray: Function value from the wrapped callable.
        """
        val = f(v)
        history.append(float(np.linalg.norm(val)))
        return val

    return wrapped, history


def run_demo_imsl_broyden() -> Dict[str, object]:
    """Compare NEQBF (Broyden FD) and NEQBJ (Broyden analytic Jacobian).

    Solves the Rosenbrock system f = [1-x, 10*(y-x²)] = 0 from starting
    point (-0.5, 0.5) using both Broyden variants and plots convergence.

    Args:
        None

    Returns:
        Dict[str, object]: Keys ``sol_fd`` (RootResult), ``sol_jac``
            (RootResult), ``history_fd`` (list of float residual norms),
            ``history_jac`` (list of float residual norms), ``plot_path``
            (str).
    """
    x0 = np.array([-0.5, 0.5])

    # ---- NEQBF: Broyden with finite-difference Jacobian ----------------
    f_fd, hist_fd = _make_tracked(rosenbrock_system)
    result_fd = solve_system_broyden_fd(f_fd, x0.copy(), x_tol=1e-8, max_iter=200)

    # ---- NEQBJ: Broyden with analytic Jacobian (hybr backend) ----------
    f_jac, hist_jac = _make_tracked(rosenbrock_system)
    result_jac = solve_system_broyden(
        f_jac, rosenbrock_jacobian, x0.copy(), x_tol=1e-8, max_iter=200
    )

    # ---- Console summary -----------------------------------------------
    print("\nIMSL Broyden Examples: Rosenbrock System")
    print("  f1(x, y) = 1 - x")
    print("  f2(x, y) = 10*(y - x²)")
    print("  Exact solution: (1, 1),  starting point: (-0.5, 0.5)")
    print("-" * 65)
    print(f"{'Method':<28} {'x':>10} {'y':>10} {'Residual':>10} {'Feval':>6}")
    print("-" * 65)
    for label, res in [("NEQBF (Broyden FD)", result_fd), ("NEQBJ (Broyden Analytic Jac)", result_jac)]:
        print(
            f"{label:<28} {res.x[0]:>10.7f} {res.x[1]:>10.7f}"
            f" {res.fval:>10.2e} {len(hist_fd if 'FD' in label else hist_jac):>6}"
        )
    print("-" * 65)
    print(f"NEQBF converged: {result_fd.success}  ({result_fd.message})")
    print(f"NEQBJ converged: {result_jac.success}  ({result_jac.message})")

    # ---- Plot: convergence curves ---------------------------------------
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_broyden.svg"

    fig, ax = plt.subplots(figsize=(8, 5))

    iters_fd = np.arange(1, len(hist_fd) + 1)
    iters_jac = np.arange(1, len(hist_jac) + 1)

    ax.semilogy(iters_fd, hist_fd, color="#0e7490", linewidth=2.0,
                marker="o", markersize=4, label="NEQBF — Broyden (FD Jacobian)")
    ax.semilogy(iters_jac, hist_jac, color="#d62728", linewidth=2.0,
                linestyle="--", marker="s", markersize=4,
                label="NEQBJ — Broyden (Analytic Jacobian / hybr)")

    ax.axhline(1e-8, color="#aaaaaa", linewidth=1.0, linestyle=":", label="Tolerance 1e-8")
    ax.set_xlabel("Function Evaluation")
    ax.set_ylabel("Residual Norm  ‖f(x)‖₂")
    ax.set_title("Broyden Quasi-Newton Convergence — Rosenbrock System")
    ax.legend(loc="upper right")
    ax.grid(True, which="both", alpha=0.25)
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {
        "sol_fd": result_fd,
        "sol_jac": result_jac,
        "history_fd": hist_fd,
        "history_jac": hist_jac,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_broyden()

Plot Output

Residual norm \(\|f(x)\|_2\) versus function evaluation count for both Broyden variants. NEQBJ (analytic Jacobian, hybr backend) converges in roughly an order of magnitude fewer evaluations than NEQBF (quasi-Newton FD update).

Console Output

IMSL Broyden Examples: Rosenbrock System
  f1(x, y) = 1 - x
  f2(x, y) = 10*(y - x┬▓)
  Exact solution: (1, 1),  starting point: (-0.5, 0.5)
-----------------------------------------------------------------
Method                                x          y   Residual  Feval
-----------------------------------------------------------------
NEQBF (Broyden FD)            1.0000000  1.0000000   3.17e-09    254
NEQBJ (Broyden Analytic Jac)  1.0000000  1.0000000   0.00e+00     10
-----------------------------------------------------------------
NEQBF converged: True  (Converged)
NEQBJ converged: True  (Converged)