Nonlinear Least Squares Examples

Example Code

Representative script:

"""IMSL UNLSF example: fit y = a*exp(b*x) to data using nonlinear least squares.

Demonstrates nonlinear_least_squares_fd on a curve-fitting problem:
- Model: y = a * exp(b * x)
- Parameters: a, b
- 5 synthetic data points with known true values a=2.0, b=0.5

Outputs:
- Table printed to stdout
- SVG plot saved to test_output/demo_imsl_unlsf.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from optimization import nonlinear_least_squares_fd


# True parameters for synthetic data
_TRUE_A = 2.0
_TRUE_B = 0.5
_X_DATA = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
_Y_DATA = _TRUE_A * np.exp(_TRUE_B * _X_DATA)


def run_demo_imsl_unlsf() -> Dict[str, object]:
    """Run IMSL UNLSF example: fit exponential model y=a*exp(b*x).

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``a``, ``b``, ``fval``,
            ``n_fev``, ``residuals``, and ``plot_path``.
    """
    def residuals(params: np.ndarray) -> np.ndarray:
        a, b = params[0], params[1]
        return a * np.exp(b * _X_DATA) - _Y_DATA

    x0 = np.array([1.0, 1.0])
    result = nonlinear_least_squares_fd(
        residuals, m=len(_X_DATA), x_guess=x0, max_iter=200, max_fev=1000
    )

    a_fit = float(result.x[0])
    b_fit = float(result.x[1])
    res_vals = residuals(result.x)

    print("\nIMSL UNLSF Example: Nonlinear Least Squares (FD Jacobian)")
    print("Model: y = a * exp(b * x)")
    print("-" * 55)
    print(f"{'Parameter':<20} {'True':>10} {'Fitted':>15}")
    print("-" * 55)
    print(f"{'a':<20} {_TRUE_A:>10.4f} {a_fit:>15.6f}")
    print(f"{'b':<20} {_TRUE_B:>10.4f} {b_fit:>15.6f}")
    print("-" * 55)
    print(f"\n{'Residuals at solution:'}")
    for xi, yi, ri in zip(_X_DATA, _Y_DATA, res_vals):
        print(f"  x={xi:.1f}  y_data={yi:.4f}  residual={ri:.2e}")
    print(f"\nObjective (sum sq / 2): {result.fval:.6e}")
    print(f"Function evaluations:  {result.n_fev}")
    print(f"Converged:             {result.success}")

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "demo_imsl_unlsf.svg"

    x_fine = np.linspace(_X_DATA[0] - 0.2, _X_DATA[-1] + 0.2, 200)
    y_true = _TRUE_A * np.exp(_TRUE_B * x_fine)
    y_fit = a_fit * np.exp(b_fit * x_fine)

    fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=True)
    axes[0].plot(x_fine, y_true, "k--", linewidth=1.5, label="True model")
    axes[0].plot(x_fine, y_fit, color="#1f77b4", linewidth=2.0, label=f"Fit: a={a_fit:.4f}, b={b_fit:.4f}")
    axes[0].scatter(_X_DATA, _Y_DATA, color="#d62728", s=50, zorder=5, label="Data")
    axes[0].set_ylabel("y")
    axes[0].set_title("IMSL UNLSF: Nonlinear Least Squares Curve Fit")
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)

    axes[1].stem(_X_DATA, res_vals, linefmt="#ff7f0e", markerfmt="o", basefmt="k-")
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("Residual")
    axes[1].set_title("Residuals at Solution")
    axes[1].grid(True, alpha=0.3)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {
        "a": a_fit,
        "b": b_fit,
        "fval": result.fval,
        "n_fev": result.n_fev,
        "residuals": res_vals,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_unlsf()

Input (Console)

Run the nonlinear least squares script from the package root:

python examples/example_imsl_unlsf.py

Plot Output

Generated SVG plot:

UNLSF: exponential curve fitting with nonlinear least squares

Output Console

Summary console output:

IMSL UNLSF Example: Nonlinear Least Squares (FD Jacobian)
Model: y = a * exp(b * x)
-------------------------------------------------------
Parameter                  True          Fitted
-------------------------------------------------------
a                        2.0000        2.000000
b                        0.5000        0.500000
-------------------------------------------------------

Residuals at solution:
  x=0.0  y_data=2.0000  residual=0.00e+00
  x=0.5  y_data=2.5681  residual=0.00e+00
  x=1.0  y_data=3.2974  residual=0.00e+00
  x=1.5  y_data=4.2340  residual=0.00e+00
  x=2.0  y_data=5.4366  residual=0.00e+00

Objective (sum sq / 2): 0.000000e+00
Function evaluations:  6
Converged:             True