ENH: interpolate: Python bivariate spline fitter `regrid_python`

Contributor & Development Discussion scipy
1

Overview

We’re excited to present a new addition to scipy.interpolate: a Python + C++ based implementation of the classic REGRID workflow for fitting 2-D B-spline surfaces to gridded data.

The new regrid_python function provides a readable, maintainable, and Fortran-free alternative to the existing implementation, while maintaining full mathematical and API compatibility with RectBivariateSpline. This work has been developing across iterations and community feedback, see the detailed discussion in the prototyping PR in my fork.

Final PR at SciPy Repository: gh-23962

What’s New?

Core Functionality

The PR introduces two private API functions to scipy.interpolate._regrid_python:

  • regrid_python(x, y, z, *, bbox=None, kx=3, ky=3, s=0.0, maxit=50)
    Fits a tensor-product B-spline surface to gridded (x, y, z) data. Key features:

    • Adaptive knot-growth strategy guided by residual magnitude
    • Automatic smoothing parameter search (when s is specified)
    • Returns an NdBSpline object for efficient evaluation
    • Full FITPACK mathematical compatibility

  • ndbspline_call_like_bivariate(ndbs, x, y, dx=0, dy=0, grid=True)
    A thin wrapper enabling the returned NdBSpline to be evaluated using the familiar RectBivariateSpline.__call__ calling pattern.

Design Highlights

The implementation decomposes the algorithm into explicit, composable steps:

  1. Separable two-pass solves: Fits are built by solving 1-D problems in the x- and y-directions using FITPACK kernels (C++ implementations in __fitpack.cc), rather than constructing a single monolithic system.

  2. Knot-growth strategy: Residuals are projected into row and column spaces; knots are inserted where the residual is largest, respecting storage constraints.

  3. Smoothing parameter search: Uses the same penalty scaling as FITPACK (penalty = 1/p), with special handling for interpolatory fits (p = -1).

  4. Banded linear algebra: Leverages banded QR decomposition and back-substitution for efficiency, converting to CSR format only when needed for residual evaluation.

  5. Clear mathematics: The algorithm explicitly minimizes
    |A \mathbf{c} - \mathbf{z}|^2 + \frac{1}{p}|D \mathbf{c}|^2,
    with full documentation of the penalty convention and knot-insertion heuristics.

Why This Matters

  • Readability: Every step is written in Python or carefully documented, making the algorithm transparent for debugging, teaching, and future optimization.
  • Maintainability: No external Fortran dependency; the code is self-contained and well-commented.
  • Extensibility: The modular structure makes it straightforward to experiment with variants (e.g., knot-insertion strategies).
  • Compatibility: Passes the same test suite as RectBivariateSpline for supported scenarios; input validation and output format are identical.

Next Steps

This is currently a private API. The team plans to eventually replace the existing RectBivariateSpline backend with this implementation. Timeline and public API surface will be coordinated with the maintainers.

To experiment with this work, you can import from the private module:

from scipy.interpolate._regrid_python import regrid_python, ndbspline_call_like_bivariate

Feedback Welcome

We’d like to hear from the community! Please share:

  • API design: Any suggestions for the eventual public interface?
  • Extensions: What features or variants would be most valuable?

Discussion are ongoing in PR #23962.

Technical Details: See the detailed PR description and module docstring for a deep dive into the mathematics, knot-insertion strategy, smoothing parameter search, and design trade-offs.

Thanks to @ev-br @j-bowhay @ilayn for their reviews/feedback on the PRs.

2 Likes
2

Public API Proposals for Exposing regrid_python

Based on the discussion in PR #23962, here are two complementary public API proposals for exposing the Python-based bivariate spline fitting implementation:

Proposal 1: Public Function rect_bivariate_spline()

Description:
Expose regrid_python() as a public function renamed to rect_bivariate_spline(). This provides a functional interface for fitting bivariate splines on rectangular grids.

Function Signature:

def rect_bivariate_spline(x, y, z, *, bbox=None, kx=3, ky=3, s=0.0, maxit=50):
    """
    Fit a bivariate spline surface to gridded data.

    Parameters
    ----------
    x, y : array_like
        Strictly increasing 1-D coordinate vectors.
    z : array_like, shape (len(x), len(y))
        Data grid.
    bbox : sequence of 4 scalars or None, optional
        Bounding box (xb, xe, yb, ye). Use None to skip clipping.
    kx, ky : int, optional
        Spline degrees (default 3 for cubic).
    s : float, optional
        Target smoothing residual (fp target). s=0 means interpolation.
    maxit : int, optional
        Maximum iterations for adaptive knot insertion (default 50).

    Returns
    -------
    NdBSpline
        The fitted 2D spline surface.
    """

Pros:

  • Simple, lightweight functional interface
  • Complements existing OOP classes (like RectBivariateSpline)
  • Clear that it returns modern NdBSpline object
  • Users can directly leverage NdBSpline API for advanced usage

Cons:

  • Doesn’t replace existing RectBivariateSpline class
  • Users must learn new function name if migrating from old API
  • Less discoverability compared to class-based interface
  • May create confusion about which to use (old class vs. new function)

Proposal 2: New Class RectBivariateSplineNdBSpline

Description:
Create a new class that wraps the Python-based implementation, using regrid_python() internally. This provides a class-based interface identical to RectBivariateSpline but returns and stores an NdBSpline object.

Class Signature:

class RectBivariateSplineNdBSpline(BivariateSpline):
    """
    Bivariate spline approximation over a rectangular mesh (NdBSpline-based).

    This is a modern implementation of bivariate spline fitting that returns
    an NdBSpline object instead of traditional (knots, coefficients) tuples.
    Uses Python-based FITPACK composition for improved transparency and control.

    Can be used for both smoothing and interpolating data.

    Parameters
    ----------
    x, y : array_like
        1-D arrays of coordinates in strictly ascending order.
    z : array_like, shape (x.size, y.size)
        2-D grid of data values.
    bbox : array_like, optional
        Sequence of length 4 specifying the boundary of the rectangular
        approximation domain. Default: [min(x), max(x), min(y), max(y)].
    kx, ky : int, optional
        Degrees of the bivariate spline. Default is 3 (cubic).
    s : float, optional
        Positive smoothing parameter. s=0 gives interpolation (default).
    maxit : int, optional
        Maximum iterations for adaptive knot insertion. Default is 50.

    Returns
    -------
    An NdBSpline-based spline object with __call__ method.

    Notes
    -----
    This class uses adaptive knot insertion internally and returns a modern
    NdBSpline object, making it suitable for advanced use cases and future
    API enhancements.

    See Also
    --------
    RectBivariateSpline : Classic FITPACK-based rectangular bivariate spline
    BivariateSpline : Base class for unstructured data
    NdBSpline : Modern n-dimensional B-spline representation
    """

    def __init__(self, x, y, z, bbox=[None] * 4, kx=3, ky=3, s=0, maxit=50):

        # Use Python-based implementation, all the input validation happens inside `regrid_python`
        self._ndbs = _regrid_python.regrid_python(
            x, y, z, bbox=bbox, kx=kx, ky=ky, s=s, maxit=maxit)

        # Store metadata
        self.degrees = kx, ky

    def __call__(self, x, y, dx=0, dy=0, grid=True):
        """
        Evaluate the spline or its derivatives.

        Parameters
        ----------
        x, y : array_like
            Evaluation coordinates.
        dx, dy : int, optional
            Derivative orders (default 0).
        grid : bool, optional
            If True, evaluate on Cartesian product (default).
            If False, treat as paired coordinates.

        Returns
        -------
        ndarray
            Evaluated spline values.
        """
        return _regrid_python.ndbspline_call_like_bivariate(
            self._ndbs, x, y, dx=dx, dy=dy, grid=grid)

Pros:

  • Familiar class-based interface (like RectBivariateSpline)
  • Clear naming: the suffix NdBSpline explicitly indicates what it returns
  • Smooth migration path for users (same signature, different internals)
  • Can coexist with RectBivariateSpline without breaking changes
  • Enables gradual deprecation of old class over time
  • Users get benefits of NdBSpline (modern, composable, extendable)
  • Two-phase migration: new API first, then deprecation if desired

Cons:

  • Adds another class to the public API namespace
  • Creates two similar classes that might confuse users initially
  • Requires deprecation communication/documentation to explain the difference

Keep ndbspline_call_like_bivariate Private (Recommended)

Description:
Keep ndbspline_call_like_bivariate() as a private utility function (prefixed with _). It’s primarily useful as an internal helper for compatibility and is not needed as a public API since users can directly call NdBSpline.__call__().

Rationale:

# What users would do with public ndbspline_call_like_bivariate:
result = ndbspline_call_like_bivariate(ndbs, x, y, dx=1, grid=True)

# What they can do directly with NdBSpline (more flexible):
X, Y = np.meshgrid(x, y, indexing='ij')
xi = np.stack([X, Y], axis=-1)
result = ndbs(xi, nu=(1, 0), extrapolate=ndbs.extrapolate)

Pros:

  • Reduces API surface area
  • NdBSpline.__call__() is more powerful and flexible
  • Avoids API duplication
  • Simplifies documentation burden
  • Keep it internal for now; can be exposed later if demand warrants

Cons:

  • Users migrating from old bivariate classes lose convenience wrapper
  • Slightly higher barrier for users unfamiliar with NdBSpline API
  • Could expose later if needed (no backward compatibility concern)

Recommended Implementation Strategy

Immediate (after ENH: interpolate: Python bivariate spline fitter `regrid_python` by czgdp1807 · Pull Request #23962 · scipy/scipy · GitHub is merged):

  1. Introduce RectBivariateSplineNdBSpline class – users can opt-in to new behavior
  2. Keep ndbspline_call_like_bivariate private with _ prefix

Future:

  1. Consider deprecating RectBivariateSpline in favor of RectBivariateSplineNdBSpline
  2. Possibly expose ndbspline_call_like_bivariate if demand justifies it

This approach:

  • Maintains backward compatibility (no breaking changes)
  • Provides modern API options for new code
  • Enables clean migration path over time
  • Keeps implementation details private