Least-Squares, Smoothing, and Rational Approximation Examples¶
Example Code¶
Representative script:
"""IMSL LSQ/smoothing/rational family examples with plots.
Covers:
- RLINE, RCURV, FNLSQ
- BSLSQ, BSVLS
- CSSED, CSSMH, CSSCV
- RATCH
"""
from __future__ import annotations
from pathlib import Path
from typing import Any, Dict
import sys
import matplotlib.pyplot as plt
import numpy as np
REPO_ROOT = Path(__file__).resolve().parent.parent
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
from interpolation import (
bspline_evaluate,
bspline_optimal_knots,
cubic_spline_evaluate,
cubic_spline_smooth,
cubic_spline_smooth_cv,
cubic_spline_smooth_error,
least_squares_bspline_fixed_knots,
least_squares_bspline_variable_knots,
least_squares_general,
least_squares_linear,
least_squares_polynomial,
rational_chebyshev,
)
def run_demo_imsl_lsq_smooth_ratch_family() -> Dict[str, Any]:
"""Run website-style demos for LSQ/smoothing/rational families.
Returns:
Dict[str, Any]: Results and SVG artifact paths.
"""
out: Dict[str, Any] = {}
out_dir = Path("test_output")
out_dir.mkdir(parents=True, exist_ok=True)
# RLINE / RCURV / FNLSQ --------------------------------------------------
x_lin = np.linspace(-1.0, 1.0, 41)
y_lin_true = 1.5 + 2.0 * x_lin
y_lin = y_lin_true + 0.03 * np.sin(9.0 * x_lin)
rline = least_squares_linear(x_lin, y_lin)
rcurv = least_squares_polynomial(x_lin, y_lin, degree=3)
fnlsq = least_squares_general(
x_lin,
y_lin,
[
lambda t: np.ones_like(t),
lambda t: t,
lambda t: t**2,
lambda t: t**3,
],
)
p3 = np.asarray(rcurv["coefficients"], dtype=float)
y_rcurv = np.polyval(p3, x_lin)
y_fnlsq = (
fnlsq["coefficients"][0]
+ fnlsq["coefficients"][1] * x_lin
+ fnlsq["coefficients"][2] * x_lin**2
+ fnlsq["coefficients"][3] * x_lin**3
)
out["RLINE"] = rline
out["RCURV"] = {
"coefficients": p3,
"residual_norm": float(rcurv["residual_norm"]),
}
out["FNLSQ"] = {
"coefficients": np.asarray(fnlsq["coefficients"], dtype=float),
"residual_norm": float(fnlsq["residual_norm"]),
}
fig1, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(x_lin, y_lin, "o", markersize=3.5, alpha=0.8, label="Data")
ax1.plot(x_lin, rline["intercept"] + rline["slope"] * x_lin, linewidth=2.0, label="RLINE")
ax1.plot(x_lin, y_rcurv, "--", linewidth=1.8, label="RCURV (deg=3)")
ax1.plot(x_lin, y_fnlsq, ":", linewidth=2.2, label="FNLSQ")
ax1.set_title("IMSL RLINE / RCURV / FNLSQ")
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.grid(True, alpha=0.3)
ax1.legend(loc="best")
fig1.tight_layout()
p1 = out_dir / "demo_imsl_lsq_family.svg"
fig1.savefig(p1, format="svg")
plt.close(fig1)
out["lsq_plot_path"] = str(p1)
# BSLSQ / BSVLS ----------------------------------------------------------
xb = np.linspace(0.0, 1.0, 51)
yb = np.sin(2.0 * np.pi * xb)
k = 4
n_coeff = 12
knots = bspline_optimal_knots(xb, order=k, n_coeff=n_coeff)
rep_bslsq = least_squares_bspline_fixed_knots(xb, yb, knots=knots, order=k)
rep_bsvls = least_squares_bspline_variable_knots(xb, yb, order=k, n_coeff=n_coeff)
x_eval = np.linspace(0.0, 1.0, 101)
y_true = np.sin(2.0 * np.pi * x_eval)
y_bslsq = np.array([bspline_evaluate(rep_bslsq, float(x)) for x in x_eval], dtype=float)
y_bsvls = np.array([bspline_evaluate(rep_bsvls, float(x)) for x in x_eval], dtype=float)
out["BSLSQ"] = {
"rmse": float(np.sqrt(np.mean((y_bslsq - y_true) ** 2))),
"sample_value": float(bspline_evaluate(rep_bslsq, 0.37)),
}
out["BSVLS"] = {
"rmse": float(np.sqrt(np.mean((y_bsvls - y_true) ** 2))),
"sample_value": float(bspline_evaluate(rep_bsvls, 0.37)),
}
fig2, ax2 = plt.subplots(figsize=(8, 5))
ax2.plot(x_eval, y_true, "k--", linewidth=1.4, label="True sin(2pi x)")
ax2.plot(x_eval, y_bslsq, linewidth=2.0, label="BSLSQ")
ax2.plot(x_eval, y_bsvls, linewidth=1.8, label="BSVLS")
ax2.set_title("IMSL BSLSQ / BSVLS")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.grid(True, alpha=0.3)
ax2.legend(loc="best")
fig2.tight_layout()
p2 = out_dir / "demo_imsl_bslsq_bsvls_family.svg"
fig2.savefig(p2, format="svg")
plt.close(fig2)
out["bslsq_plot_path"] = str(p2)
# CSSED / CSSMH / CSSCV --------------------------------------------------
# Use random Gaussian noise (seeded) so the three methods produce
# visually distinguishable results at different smoothing strengths.
rng = np.random.default_rng(42)
xs = np.linspace(0.0, 1.0, 120)
y_clean = np.sin(2.0 * np.pi * xs)
y_noisy = y_clean + 0.22 * rng.standard_normal(len(xs))
# CSSED: very light smoothing (error_scale << actual noise) -- stays close to data
pp_cssed = cubic_spline_smooth_error(xs, y_noisy, error_scale=0.05)
# CSSCV: cross-validation picks the balanced optimum
cv = cubic_spline_smooth_cv(xs, y_noisy)
pp_csscv = cv["spline"]
# CSSMH: heavy regularisation (large lam) -- over-smoothed, polynomial-like
pp_cssmh = cubic_spline_smooth(xs, y_noisy, lam=0.5)
y_cssed = np.array([cubic_spline_evaluate(pp_cssed, float(x)) for x in xs], dtype=float)
y_cssmh = np.array([cubic_spline_evaluate(pp_cssmh, float(x)) for x in xs], dtype=float)
y_csscv = np.array([cubic_spline_evaluate(pp_csscv, float(x)) for x in xs], dtype=float)
out["CSSED"] = {
"rmse_to_clean": float(np.sqrt(np.mean((y_cssed - y_clean) ** 2))),
"sample_value": float(cubic_spline_evaluate(pp_cssed, 0.33)),
}
out["CSSMH"] = {
"rmse_to_clean": float(np.sqrt(np.mean((y_cssmh - y_clean) ** 2))),
"sample_value": float(cubic_spline_evaluate(pp_cssmh, 0.33)),
}
out["CSSCV"] = {
"lambda": float(cv["lambda"]),
"gcv_scores_min": float(np.min(cv["gcv_scores"])),
"rmse_to_clean": float(np.sqrt(np.mean((y_csscv - y_clean) ** 2))),
}
fig3, (ax3, ax3_resid) = plt.subplots(
2,
1,
figsize=(8, 7),
sharex=True,
gridspec_kw={"height_ratios": [3, 1.4]},
)
ax3.scatter(xs, y_noisy, s=4, color="#aaaaaa", zorder=1, label="Noisy data")
ax3.plot(xs, y_clean, "k--", linewidth=1.5, label="Clean reference", zorder=2)
# CSSED drawn last so it sits on top -- it tightly tracks the data
ax3.plot(xs, y_cssmh, color="#d62728", linewidth=2.2, linestyle="-", label="CSSMH (λ=0.5, over-smooth)", zorder=3)
ax3.plot(xs, y_csscv, color="#1f77b4", linewidth=2.0, linestyle="--", label="CSSCV (auto-optimal)", zorder=4)
ax3.plot(xs, y_cssed, color="#ff7f0e", linewidth=1.6, linestyle="-.", label="CSSED (err_scale=0.05, light)", zorder=5)
ax3.set_title("IMSL CSSED / CSSMH / CSSCV — three smoothing levels")
ax3.set_ylabel("y")
ax3.grid(True, alpha=0.3)
ax3.legend(loc="upper right", fontsize=8.5)
ax3_resid.plot(xs, y_cssmh - y_clean, color="#d62728", linewidth=1.6, linestyle="-", label="CSSMH - clean")
ax3_resid.plot(xs, y_csscv - y_clean, color="#1f77b4", linewidth=1.6, linestyle="--", label="CSSCV - clean")
ax3_resid.plot(xs, y_cssed - y_clean, color="#ff7f0e", linewidth=1.4, linestyle="-.", label="CSSED - clean")
ax3_resid.axhline(0.0, color="#333333", linewidth=0.9, alpha=0.6)
ax3_resid.set_xlabel("x")
ax3_resid.set_ylabel("residual")
ax3_resid.grid(True, alpha=0.25)
ax3_resid.legend(loc="upper right", ncol=3, fontsize=8)
fig3.tight_layout()
p3 = out_dir / "demo_imsl_css_family.svg"
fig3.savefig(p3, format="svg")
plt.close(fig3)
out["css_plot_path"] = str(p3)
# RATCH ------------------------------------------------------------------
xr = np.linspace(0.1, 2.5, 80)
yr = np.log1p(xr)
ratch = rational_chebyshev(xr, yr, numerator_degree=3, denominator_degree=2)
pcoef = np.asarray(ratch["numerator"], dtype=float)
qcoef = np.asarray(ratch["denominator"], dtype=float)
num = np.polyval(pcoef[::-1], xr)
den = np.polyval(qcoef[::-1], xr)
yhat = num / den
out["RATCH"] = {
"numerator": pcoef,
"denominator": qcoef,
"max_abs_error": float(ratch["max_abs_error"]),
"rmse": float(np.sqrt(np.mean((yhat - yr) ** 2))),
}
fig4, ax4 = plt.subplots(figsize=(8, 5))
ax4.plot(xr, yr, "k--", linewidth=1.4, label="log(1+x)")
ax4.plot(xr, yhat, linewidth=2.0, label="RATCH")
ax4.set_title("IMSL RATCH")
ax4.set_xlabel("x")
ax4.set_ylabel("y")
ax4.grid(True, alpha=0.3)
ax4.legend(loc="best")
fig4.tight_layout()
p4 = out_dir / "demo_imsl_ratch.svg"
fig4.savefig(p4, format="svg")
plt.close(fig4)
out["ratch_plot_path"] = str(p4)
return out
if __name__ == "__main__":
run_demo_imsl_lsq_smooth_ratch_family()
Input (Console)¶
Run the least-squares and smoothing scripts from the package root:
python examples/example_imsl_lsq_smooth_ratch_family.py
python examples/example_imsl_bsls2.py
python examples/example_imsl_bsls3.py
Plot Output¶
Generated SVG plots:
Output Console¶
Summary console output:
script input setup representative numeric outputs
example_imsl_lsq_smooth_ratch_family.py 41 linear, 51 B-spline, 120 smoothing, and 80 rational samples RLINE slope = 2.00961673409346; intercept = 1.5; BSLSQ RMSE = 0.0001841959679721844; BSVLS RMSE = 0.0001841959679721844; CSSMH RMSE = 0.08471105728072979; CSSCV lambda = 1e-06; RATCH max abs error = 7.517921174385322e-06
example_imsl_bsls2.py 4x6 parity grid and 90x120 plotting grid mean abs error = 0.495449311150341; max abs error = 0.9458701788744891; spline head = [-0.1799418682616353, 0.3792275358671994, 0.3076189352441675]
example_imsl_bsls3.py 24 parity points from structured 3D sample axes max abs error = 0.431595810457456; true head = [0, 0, 2.287355287178842]; spline head = [-0.2852155134461971, -0.4107214322657342, 2.385896890218667]
Detailed integration result tables:
Least-squares and smoothing full metric output
parameter value
RLINE_intercept 1.5
RLINE_slope 2.00961673409346
RLINE_residual_norm 0.1323833248001724
FNLSQ_residual_norm 0.1216957656404049
BSLSQ_rmse 0.0001841959679721844
BSVLS_rmse 0.0001841959679721844
CSSED_rmse_to_clean 0.08471717177146082
CSSMH_rmse_to_clean 0.08471105728072979
CSSCV_lambda 1e-06
CSSCV_gcv_scores_min 2.188039866410482e-09
CSSCV_rmse_to_clean 0.08472385408477985
RATCH_max_abs_error 7.517921174385322e-06
RATCH_rmse 1.511313183359482e-06
BSLS2 parity full output
x y true spline error
0 0 0 -0.1799418682616353 0.1799418682616353
0 1 0.8414709848078965 0.3792275358671994 0.4622434489406971
0 2 0.9092974268256817 0.3076189352441675 0.6016784915815142
0 3 0.1411200080598672 0.2671787090079698 -0.1260587009481026
0 4 -0.7568024953079282 -1.164923584215012 0.4081210889070839
0 5 -0.9589242746631385 -1.28583913180298 0.326914857139842
1 0 2.287355287178842 1.659455254615714 0.6279000325631285
1 1 2.471726672004819 3.002081942436757 -0.5303552704319383
1 2 0.383603953541131 1.218261310226298 -0.8346573566851674
1 3 -2.057202470728003 -1.825630875575844 -0.2315715951521593
1 4 -2.606626430685079 -2.529560676423036 -0.07706575426204365
1 5 -0.7595300713439712 0.1863401075305179 -0.9458701788744891
2 0 6.71884969742825 6.002021306741963 0.7168283906862873
2 1 1.042743656235904 0.2823733236350315 0.7603703326008729
2 2 -5.592056093640982 -6.291170752235975 0.6991146585949934
2 3 -7.085545260112314 -6.67610202450181 -0.409443235610504
2 4 -2.064616791102519 -1.199053528414262 -0.8655632626882572
2 5 4.854510834178772 5.020539626351583 -0.1660287921728116
3 0 2.834471132487004 2.18706468575647 0.6474064467305345
3 1 -15.20078446306795 -15.94778035435259 0.7469958912846337
3 2 -19.26050892528742 -19.39822379337659 0.1377148680891693
3 3 -5.612210305985402 -5.204139232138545 -0.4080710738468571
3 4 13.19592858660572 12.68771262206538 0.5082159645403372
3 5 19.87179159281414 20.34444349982927 -0.4726519070151234
BSLS3 parity full output
idx x y z true spline error
0 0 0 0 0 -0.2852155134461971 0.2852155134461971
1 0 0 1 0 -0.4107214322657342 0.4107214322657342
2 0 1 0 2.287355287178842 2.385896890218667 -0.09854160303982429
3 0 1 1 0.8414709848078965 0.7504540402370715 0.09101694457082499
4 0 2 0 6.71884969742825 6.931600531594422 -0.2127508341661724
5 0 2 1 2.471726672004819 2.718618662888146 -0.2468919908833276
6 1 0 0 0.8414709848078965 0.9422441184010596 -0.1007731335931631
7 1 0 1 0.3095598756531122 0.09369495670793548 0.2158649189451767
8 1 1 0 2.471726672004819 2.174065888679456 0.2976607833253628
9 1 1 1 0.9092974268256817 0.8209165176639955 0.08838090916168617
10 1 2 0 1.042743656235904 1.061582765998652 -0.01883910976274783
11 1 2 1 0.383603953541131 0.3561025203943294 0.02750143314680165
12 2 0 0 0.9092974268256817 1.340893237283138 -0.431595810457456
13 2 0 1 0.3345118292392623 0.422696409118712 -0.08818457987944972
14 2 1 0 0.383603953541131 0.05264012478865676 0.3309638287524743
15 2 1 1 0.1411200080598672 0.5105375757369264 -0.3694175676770591
16 2 2 0 -5.592056093640982 -5.262807453637794 -0.3292486400031871
17 2 2 1 -2.057202470728003 -2.035379504460865 -0.02182296626713764
18 3 0 0 0.1411200080598672 -0.2757028545917273 0.4168228626515945
19 3 0 1 0.05191514970317339 -0.2932769300346445 0.3451920797378179
20 3 1 0 -2.057202470728003 -1.78506340800126 -0.2721390627267435
21 3 1 1 -0.7568024953079282 -0.8597788167492768 0.1029763214413486
22 3 2 0 -7.085545260112314 -7.499399301561888 0.4138540414495742
23 3 2 1 -2.606626430685079 -2.481258610262492 -0.1253678204225879