3 – Fitting and ML Correction

From Physics to Machine Learning

Fitting, Residuals, and Hybrid Models

NISP dispersion calibration — optimization and ML correction

1. Fitting the Physical Model

The physical model has free parameters (focal lengths, prism angles, grating period…).
These must be adjusted to match the calibration data.

We have:

observations: measured centroids (y_mosaic, z_mosaic) for each source position and wavelength
predictions: model output (ŷ, ẑ) = model(y_nisp, z_nisp, λ; θ)
parameters: vector θ to optimize

Goal: find θ* that minimizes the discrepancy between predictions and observations.

The Cost Function

The cost function measures how far the model is from the data. A standard choice is the Mean Squared Error (MSE):

\[\mathcal{L}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \left[ (\hat{y}_i - y_i)^2 + (\hat{z}_i - z_i)^2 \right]\]

Why MSE?
— Penalizes large errors more than small ones
— Smooth and differentiable → easy to minimize

Alternatives
— MAE: more robust to outliers
— Weighted MSE: if some measurements are noisier
— Sum over y and z separately

A good cost function reflects what you care about physically.
Here: minimizing the position prediction error.

What is SciPy?

SciPy is a Python library built on top of NumPy, providing algorithms for scientific computing.

Key submodules
— scipy.optimize: minimization, root finding, curve fitting
— scipy.linalg: linear algebra
— scipy.interpolate: interpolation
— scipy.stats: statistical distributions and tests

Why use it here?
— scipy.optimize.minimize provides general-purpose optimizers
— No need to implement gradient descent by hand
— Works seamlessly with NumPy arrays

pip install scipy          # install
from scipy.optimize import minimize  # import

Full documentation: https://docs.scipy.org/doc/scipy/

Choosing an Optimizer

Use scipy.optimize.minimize. The choice of method depends on the problem.

Method	Needs gradient?	Good when
`Nelder-Mead`	No	Few parameters, noisy cost
`L-BFGS-B`	Numerical	Many parameters, smooth cost
`Powell`	No	Medium size, no gradient

Start with Nelder-Mead or Powell: robust for ~5–15 parameters, no derivatives needed.

from scipy.optimize import minimize

result = minimize(cost_function, initial_params, method="Nelder-Mead",
                  options={"maxiter": 5000, "xatol": 1e-6, "fatol": 1e-8})

Docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

Structure of the Optimization

The optimizer calls your cost function many times, varying θ. You need to:

1. Define a parameter vector θ (a 1D NumPy array)
2. Write cost(θ): unpack θ → run model forward → return scalar MSE
3. Choose initial values θ₀ (physically motivated)
4. Call scipy.optimize.minimize(cost, θ₀, ...)
5. Extract result.x (optimal parameters) and result.fun (final cost)

Important: the cost function must be vectorized — process all data points at once with NumPy, not in a Python loop. Otherwise optimization will be too slow.

Exercise 1 — Implement the Optimization

Implement the full optimization pipeline for a single dataset (e.g. rgs000_0.csv).

Steps to follow:

Load and clean the data with your functions from section 2.
Identify which model parameters are free and which are fixed.
Write a cost(theta) function: unpack theta, run the physical model on all data points, return the MSE.
Set initial values theta0 from physical knowledge (nominal focal lengths, etc.).
Run the optimizer. Print result.fun and result.x.
Compare predictions before and after optimization.

Hint: start with only 1–2 free parameters to verify the setup, then add more.

2. Evaluating the Fit

A scalar cost value alone does not tell you where or why the model fails.
You must analyze the residuals to understand the fit quality.

For each data point:

\[r_y^{(i)} = \hat{y}_i - y_i \qquad r_z^{(i)} = \hat{z}_i - z_i\]

Good fit
Residuals small, random, centered on zero, no spatial structure.

Bad fit
Residuals show a pattern — trend with position, wavelength, or detector location.

Evaluation Metrics

Metric	Formula	Unit	Interpretation
MSE	\(\frac{1}{N}\sum r_i^2\)	mm²	Sensitive to outliers.
RMSE	\(\sqrt{\text{MSE}}\)	mm	Same unit as data. Easy to interpret.
MAE	\(\frac{1}{N}\sum \|r_i\|\)	mm	More robust to outliers.
R²	\(1 - \frac{\sum r_i^2}{\sum(y_i-\bar{y})^2}\)	—	Fraction of variance explained. 1 = perfect.
Max error	\(\max \|r_i\|\)	mm	Worst-case prediction error.

For calibration: RMSE is the standard metric. A pixel on NISP is ~0.3 mm — aim for RMSE well below that.

Visualizing Residuals

Plots reveal patterns that scalar metrics cannot.

Residuals vs position
Plot r_y vs y_nisp or y_mosaic.
A trend means a systematic bias that depends on field position.

2D residual map
Scatter of (y_mosaic, z_mosaic) colored by r_y.
Reveals spatial structure (corner effects, detector tilt).

Histogram of residuals
Should be roughly Gaussian, centered on zero.
Heavy tails → outliers still present.
Offset → systematic bias.

Residuals vs wavelength
Reveals dispersion errors: model wrong about how much the grism bends each wavelength.

Exercise 3 — Compute and Visualize Residuals

After fitting, analyze the residuals.

Steps to follow:

Compute residuals r_y = pred_y - obs_y and r_z for all data points.
Compute RMSE, MAE, and max error. Print a summary.
Plot a histogram of r_y and r_z. Are they Gaussian? Centered?
Scatter plot of residuals colored by y_nisp. Is there a trend?
2D map: plot (y_mosaic, z_mosaic) colored by |r_y|.
Repeat for all 6 datasets. Do residuals have the same pattern?

Use plt.hist, plt.scatter(..., c=..., cmap=...), and plt.colorbar().
See: https://matplotlib.org/stable/gallery/index.html

3. When Physics Is Not Enough

After optimization, residuals may still show systematic structure.
Some error is beyond what the physical model can capture.

Possible causes:

Model too simple (paraxial approximation fails at large field angles)
Unknown optical aberrations not in the model
Detector geometry deviations (chip warping, gap offsets)
Parameters that are coupled and cannot be individually separated

Idea: train a ML model to predict the residual from the input coordinates.
The ML model learns the systematic error that the physics model cannot capture.

The Residual Learning Idea

inputs: (y_nisp, z_nisp, λ)
↓
Physical model → (ŷ_phys, ẑ_phys)
↓
residuals: r_y = y_obs − ŷ_phys, r_z = z_obs − ẑ_phys
↓
ML model trained on (y_nisp, z_nisp, λ) → (r_y, r_z)
↓
Hybrid prediction ŷ_hybrid = ŷ_phys + δy_ML

What is Machine Learning?

Machine Learning is a family of algorithms that learn a function from data.

Supervised learning
Learn a mapping from inputs X to outputs y given labeled examples.
Our case: predict residuals (r_y, r_z) from (y_nisp, z_nisp, λ).

Why ML here?
The physical model leaves structured residuals that cannot be captured by adjusting parameters.
ML can learn these patterns directly from data.

Key concepts
— Training set / test set
— Overfitting vs generalization
— Loss function, model capacity

FIDLE CNRS training (recommended): https://fidle.cnrs.fr/w3/sequences/01-Concepts.html

What is scikit-learn?

scikit-learn is the standard Python library for machine learning. It provides a unified API for dozens of algorithms.

Key features
— Regression, classification, clustering
— Model selection and evaluation tools
— Preprocessing pipelines
— Consistent fit / predict / score interface

Minimal usage pattern

from sklearn.ensemble import RandomForestRegressor

model = RandomForestRegressor()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

pip install scikit-learn

User guide: https://scikit-learn.org/stable/user_guide.html

Preparing the ML Dataset

Input matrix X (shape: N × 3):
columns: y_nisp, z_nisp, wavelength

Target vector y (shape: N):
r_y = y_mosaic_obs − y_mosaic_phys (same for r_z)

Important: split into train and test sets before fitting.
Never evaluate performance on data the model was trained on.
Use sklearn.model_selection.train_test_split with 80% train / 20% test.

Docs: https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html

Data Preprocessing and Feature Scaling

Many algorithms are sensitive to the scale of input features (and output variables). Always preprocess before training.

Why scaling matters
— y_nisp ∈ [−90, 90] mm, λ ∈ [1.0, 1.9] µm: very different ranges
— Algorithms that use distances or gradients (SVM, MLP, linear models) will be dominated by the largest-scale feature
— Tree-based models (RF, GBT) are not affected by scale

Common scalers
— StandardScaler: zero mean, unit variance. Best default.
— MinMaxScaler: maps to [0, 1]. Good when bounds are known.
— RobustScaler: uses median and IQR. Better when outliers are present.

Fitting the Scaler

Critical rule: fit the scaler only on the training set, then apply it to the test set.
Never fit on the full dataset — that leaks test information.

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)  # fit + transform on train
X_test_scaled  = scaler.transform(X_test)        # transform only on test

Docs: https://scikit-learn.org/stable/modules/preprocessing.html

Choosing the ML Model

Start simple. Increase complexity only if needed.

Model	Complexity	Training time	When to use
Linear Regression	Low	Instant	Baseline. Does the residual have a linear trend?
Random Forest	Medium	Seconds	Good default. Non-linear, no tuning needed.
Gradient Boosting	Medium	Seconds	Often better than RF. Try `HistGradientBoosting`.
MLP (neural net)	High	Minutes	Powerful but needs careful tuning. Use last.

Recommended order: LinearRegression → RandomForest → comparison.
All run on a laptop without GPU.

Docs: https://scikit-learn.org/stable/supervised_learning.html

Exercise 4 — Build and Evaluate the ML Model

Train a ML model to predict the residuals of your physical model.

Steps to follow:

Compute residuals r_y and r_z from your fitted physical model.
Build and scale the feature matrix X = [y_nisp, z_nisp, wavelength] and target y = r_y.
Split into train/test (80/20). Use a fixed random_state.
Train a LinearRegression. Evaluate RMSE on the test set.
Train a RandomForestRegressor. Evaluate RMSE on the test set.
Compare RMSE of both models. Which performs better?
Plot predicted vs actual residuals for the test set (plt.scatter(r_true, r_pred)).

Overfitting and Cross-Validation

A model that performs very well on training data but poorly on test data is overfitting.
It has memorized training noise instead of learning the underlying pattern.

Signs of overfitting
— Train RMSE ≪ Test RMSE
— Erratic predictions between nearby points
— Poor performance on a different dataset

How to reduce it
— Limit complexity (e.g. max_depth in RF)
— Use cross-validation to select hyperparameters
— Use more training data

from sklearn.model_selection import cross_val_score

# 5-fold cross-validation
scores = cross_val_score(model, X, y, cv=5,
                         scoring="neg_root_mean_squared_error")
print(-scores.mean(), "±", scores.std())

Docs: https://scikit-learn.org/stable/modules/cross_validation.html

4. Building the Hybrid Model

The hybrid model combines the physical prediction and the ML correction: \[\hat{y}_\text{hybrid} = \hat{y}_\text{phys}(\theta^*, \mathbf{x}) + \delta y_\text{ML}(\mathbf{x})\] \[\hat{z}_\text{hybrid} = \hat{z}_\text{phys}(\theta^*, \mathbf{x}) + \delta z_\text{ML}(\mathbf{x})\]

The physical model provides the correct global behavior (dispersion, wavelength dependence)
The ML model corrects local, structured residuals that the physics misses
Outside the training range: the ML correction should be small — verify this
The hybrid is interpretable: each component has a clear role

Exercise 5 — Implement and Compare the Hybrid Model

Build the hybrid predictor and compare it to the physical model alone.

Steps to follow:

Compute ŷ_hybrid = ŷ_phys + δy_ML on the test set.
Compute RMSE for: (a) physical model alone, (b) hybrid model.
Plot residual histograms for both models on the same figure (overlaid).
Plot 2D residual maps for both models side by side.
Is the improvement uniform across the focal plane, or localized?
Test on a different dataset (e.g. rgs180_0 if trained on rgs000_0). Does the ML correction generalize?

If the ML correction does not generalize to other datasets, it may be overfitting to dataset-specific noise.

Model Comparison Table

Fill this table with your results. Use the test set only.

Model	RMSE y (mm)	RMSE z (mm)	Generalizes?
Physical model (before fit)	—	—	—
Physical model (after fit)	—	—	Yes
Physical + Linear correction	—	—	?
Physical + Random Forest	—	—	?
Physical + Gradient Boosting	—	—	?

This table is the key output of this section. It compares physical and hybrid approaches quantitatively.

Summary

Step	Goal	Key tool
Fit physical model	Find optimal parameters θ*	`scipy.optimize.minimize`
Evaluate residuals	Identify systematic errors	RMSE, MAE, residual plots
ML correction	Predict residuals from inputs	`sklearn` regressors
Overfitting check	Verify generalization	Train/test split, cross-val
Hybrid model	Physics + ML correction	Additive combination
Comparison	Quantify improvement	RMSE table, residual maps
Next step	Full calibration pipeline (section 4)	Autonomous work