LSQRR Example — Least Squares Polynomial Fitting

Demonstrates overdetermined polynomial least-squares fitting using linearsystems.lsqrr().

Example Code

"""IMSL LSQRR example: polynomial least-squares fitting.

Generates an overdetermined system (20 equations, 3 unknowns) for fitting
a quadratic polynomial to noisy data, solves with lsqrr, and plots the fit.

True model: b = 1 + 2*t + 3*t^2 + noise

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

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from linearsystems import lsqrr


def run_demo_imsl_lsqrr() -> Dict[str, object]:
    """Run IMSL LSQRR example: overdetermined least-squares polynomial fit.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``x`` (np.ndarray — fitted
            polynomial coefficients), ``residual_norm`` (float), ``rank`` (int),
            and ``plot_path`` (str).
    """
    rng = np.random.default_rng(0)
    m = 20
    n = 3
    t = np.linspace(0.0, 1.0, m)
    x_true = np.array([1.0, 2.0, 3.0])
    noise = 0.05 * rng.standard_normal(m)
    b = 1.0 + 2.0 * t + 3.0 * t ** 2 + noise

    # Vandermonde matrix
    a = np.column_stack([np.ones(m), t, t ** 2])

    x, residuals, rank, sv = lsqrr(a, b)
    residual_norm = float(np.linalg.norm(a @ x - b))

    print("\nIMSL LSQRR Example: Overdetermined polynomial least-squares fit")
    print("-" * 55)
    print(f"  Equations m = {m},  unknowns n = {n}")
    print(f"  True coefficients  : {x_true}")
    print(f"  Fitted coefficients: {np.round(x, 6)}")
    print("-" * 55)
    print(f"{'‖Ax-b‖':<30} {residual_norm:.4e}")
    print(f"{'Effective rank':<30} {rank}")
    print(f"{'Singular values':<30} {np.round(sv, 4)}")
    print("-" * 55)

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

    t_fine = np.linspace(0.0, 1.0, 200)
    b_true_fine = x_true[0] + x_true[1] * t_fine + x_true[2] * t_fine ** 2
    b_fit_fine = x[0] + x[1] * t_fine + x[2] * t_fine ** 2

    fig, axes = plt.subplots(1, 2, figsize=(13, 5))

    # Left: data + fit
    ax = axes[0]
    ax.scatter(t, b, s=40, color="#1f77b4", zorder=5, label="Noisy data")
    ax.plot(t_fine, b_true_fine, "g--", linewidth=1.5, label="True model")
    ax.plot(t_fine, b_fit_fine, "r-", linewidth=2.0, label="Fitted curve")
    ax.set_xlabel("t")
    ax.set_ylabel("b")
    ax.set_title("LSQRR: Quadratic Polynomial Fit (m=20, n=3)")
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Right: coefficient comparison
    ax2 = axes[1]
    idx = np.arange(n)
    ax2.bar(idx - 0.2, x_true, width=0.35, label="True", color="#2ca02c", alpha=0.85)
    ax2.bar(idx + 0.2, x, width=0.35, label="Fitted", color="#d62728", alpha=0.85)
    ax2.set_xticks(idx)
    ax2.set_xticklabels(["c₀", "c₁", "c₂"])
    ax2.set_title("Coefficient Comparison")
    ax2.set_ylabel("Value")
    ax2.legend()
    ax2.grid(True, axis="y", alpha=0.3)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved to: {plot_path}")

    return {"x": x, "residual_norm": residual_norm, "rank": rank, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_lsqrr()

Plot Output

Fitted polynomial curve through noisy data and recovered coefficient comparison.

Console Output

IMSL LSQRR Example: Overdetermined polynomial least-squares fit
-------------------------------------------------------
  Equations m = 20,  unknowns n = 3
  True coefficients  : [1. 2. 3.]
  Fitted coefficients: [1.035443 1.774627 3.199016]
-------------------------------------------------------
‖Ax-b‖                         1.7200e-01
Effective rank                 3
Singular values                [5.3304 1.6376 0.2552]
-------------------------------------------------------

Plot saved to: test_output\example_imsl_lsqrr.svg