I have been working on adding weighted least-squares fitting of N-dimensional tensor-product B-splines to `scipy.interpolate`. The implementation was originally proposed in scipy/scipy#25472 and has since been split into smaller PRs. The internal scattered-data fitter was added in scipy/scipy#25493, and an internal helper for rectilinear-grid data was added in scipy/scipy#25501. The public API is currently being discussed in scipy/scipy#26213.
My original use case comes from spectroscopy. We needed to fit a point-spread function (PSF) as a smooth function of variables such as radial separation and wavelength. The observations are weighted and can be scattered rather than arranged on a complete tensor-product grid. We also want the same method to extend beyond two dimensions.
For scattered samples \(x_i \in \mathbb{R}^d\), observed values \(y_i\), residual weights \(w_i\), tensor-product B-spline basis matrix \(A\), and coefficients \(c\), the implementation solves
\[
\min_c \left\| W(Ac-y) \right\|_2,
\]
where \(W = \operatorname{diag}(w)\).
The design matrix is sparse because each B-spline basis function has local support. The fitted coefficients are stored in the existing `NdBSpline` class, which is then used for evaluation. This separates fitting from evaluation and avoids introducing another spline representation.
The current proposal requires the user to provide explicit full knot vectors. Automatic knot generation and smoothing penalties are intentionally outside the scope of this API proposal.
My preferred public interface is one function with a `grid` keyword. For scattered data, `grid=False` would be the default:
~~~python
x = np.column_stack((radius, wavelength))
spl = make_lsq_ndbspline(
x,
psf,
(radius_knots, wavelength_knots),
k=(3, 3),
w=weights,
grid=False,
)
fitted_psf = spl(x)
~~~
Here, `x` has shape `(n_samples, n_dimensions)`, and each row contains the coordinates of one observation.
For values sampled on a complete rectilinear grid, the same function would be called with `grid=True`:
~~~python
spl = make_lsq_ndbspline(
(radius_axis, wavelength_axis),
psf_values,
(radius_knots, wavelength_knots),
k=(3, 3),
w=weights,
grid=True,
)
~~~
In this case, the first argument contains the one-dimensional coordinate arrays defining the grid, and the leading dimensions of `psf_values` correspond to those grid axes.
Internally, these two cases can still use separate helpers because their input validation and coordinate construction are different. The `grid` keyword would only provide a single public entry point and select the appropriate internal implementation.
Testing the implementation on the spectroscopy use case highlighted several points:
- Weighted fitting is important because the measurement uncertainties are not uniform.
- Scattered input is essential because the observations do not always form a complete grid.
- Returning `NdBSpline` works well because evaluation is already implemented and tested there.
- Sparse construction and solving are necessary because dense tensor-product design matrices quickly become too large.
- Invalid knot placement or insufficient sample coverage can leave some basis functions unsupported, so the implementation detects this and raises an error.
- Pairwise evaluation is already possible by passing an array with shape `(n_points, n_dimensions)` to the returned `NdBSpline`.
The current PR initially proposed two public functions, one for scattered data and one for gridded data. After considering the alternatives, I prefer one public `make_lsq_ndbspline` function with `grid=False` by default. I think this keeps the public API smaller while still making both input forms directly available.
I would appreciate feedback on:
- Whether a single function with `grid=True` or `grid=False` is the preferred interface.
- Whether `grid=False` should be the default.
- Whether the meaning of `grid` is sufficiently clear.
- The proposed function name.
- Whether the input and weight conventions are consistent with the rest of `scipy.interpolate`.
- Any additional functionality that should be required before making the constructor public.