[headtr1ck]

Multidimensional interp with method="cubic" / "quintic" used scipy.interpolate.interpn, which builds the collocation matrix of the whole grid and solves it with an iterative sparse solver (gcrotmk, atol=1e-6). That is slow, needs a lot of memory and is only accurate to ~1e-6 (~1e-5 for quintic).

• Solve for the tensor-product spline coefficients one dimension after the other with make_interp_spline and evaluate them with NdBSpline. This is exact and keeps the interpn semantics (bounds_error, fill_value, NaN destinations, any destination layout).

• Single point in a (288, 71, 73) array (the example from xarray interp with method="cubic" fails for time/lon/lat interpolation #11342): 6.5 s / 1.8 GB → 0.45 s / 67 MB.

• Tests compare against scipy with a direct solver; test_interpolate_nd uses that reference for spline methods too, because quintic interpn exceeds the default rtol.

• New asv benchmarks InterpolationSplineND / InterpolationSplineNDDask (time for a point, a grid and pointwise destinations, plus peak memory).

• Compatible with Fix interp performance regression for vectorized interpolation #11681 (merges cleanly, only the interpolator passed to apply_ufunc changes).

• Closes xarray interp with method="cubic" fails for time/lon/lat interpolation #11342

• Tests added

• User visible changes (including notable bug fixes) are documented in whats-new.rst

AI Disclosure

• This PR contains AI-generated content.
  • I have tested any AI-generated content in my PR.
  • I take responsibility for any AI-generated content in my PR.
  • Tools: Claude Code with Opus 5.5

🤖 Generated with Claude Code

[headtr1ck]

Some background on why this works and why scipy doesn't do it already:

Why solving one dimension at a time is exact. On a regular grid, the values are interpolated at every combination of grid points. So the collocation matrix of the n-D spline is the Kronecker product of the 1-D collocation matrices, A = A_t ⊗ A_y ⊗ A_x, and solving A c = v is the same as solving with each 1-D matrix along its own axis. That is what calling make_interp_spline once per dimension does. Each 1-D matrix is small and banded, so this is fast, needs only a few copies of the data, and is exact. The result matches scipy with a direct solver (RegularGridInterpolator(..., solver=spsolve)) to ~1e-15.

What scipy does instead. RegularGridInterpolator / interpn build the spline with make_ndbspl (since scipy 1.13). It lists all grid points with itertools.product, builds the full sparse collocation matrix (64 nonzeros per grid point for 3-D cubic, which is where the ~2 GB in #11342 comes from), and solves it with gcrotmk at atol=1e-6. That is also why its results are only accurate to ~1e-6 (~1e-5 for quintic), and why test_interpolate_nd now needs an exact reference. make_ndbspl just treats this as a general sparse problem and doesn't take advantage of the Kronecker structure.

Downsides. None that I can see for interp:

• It needs the data on a full grid, which is always the case for interpn and interp.
• The boundary conditions are the same (not-a-knot), so the spline is identical.
• NaNs in the data still raise, like before.
• bounds_error, fill_value and NaN destinations behave like interpn.
• The only thing given up is passing a custom solver to scipy, which interpn doesn't accept anyway, so it was never possible through interp.

This could be fixed in scipy too, which would help other RegularGridInterpolator users.

[github-actions]

🏎️ New personal best! Some benchmarks got better.

Change	Before [78199ec]	After [454712d]	Ratio	Benchmark (Parameter)
🚀	7.33±0.01ms	2.66±0.02ms	0.36	dataset_io.IOReadSingleFile.time_read_dataset('scipy', None)
🚀	262±3ms	12.2±0.08ms	0.05	interp.InterpolationSplineND.time_interpolation_spline_nd('cubic', 'grid')
🚀	245±20ms	7.54±0.02ms	0.03	interp.InterpolationSplineND.time_interpolation_spline_nd('cubic', 'point')
🚀	251±10ms	7.69±0.1ms	0.03	interp.InterpolationSplineND.time_interpolation_spline_nd('cubic', 'pointwise')
🚀	2.92±0.02s	18.6±0.03ms	0.01	interp.InterpolationSplineND.time_interpolation_spline_nd('quintic', 'grid')
🚀	2.91±0s	7.75±0.02ms	0	interp.InterpolationSplineND.time_interpolation_spline_nd('quintic', 'point')
🚀	2.86±0.01s	8.71±0.8ms	0	interp.InterpolationSplineND.time_interpolation_spline_nd('quintic', 'pointwise')
🚀	244±20ms	19.0±1ms	0.08	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('cubic', 'grid')
🚀	240±20ms	13.9±0.9ms	0.06	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('cubic', 'point')
🚀	242±20ms	14.3±0.9ms	0.06	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('cubic', 'pointwise')
🚀	2.87±0.03s	25.4±0.8ms	0.01	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('quintic', 'grid')
🚀	2.85±0.03s	14.3±1ms	0.01	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('quintic', 'point')
🚀	2.85±0.04s	14.6±1ms	0.01	interp.InterpolationSplineNDDask.time_interpolation_spline_nd('quintic', 'pointwise')

🚀 better · 🐌 worse · 💥 now failing · 🩹 fixed · 🤷 not comparable

Comparing 454712dd (this PR merged into the base branch) with 78199ecb (base branch).
Only changes by more than a factor of 1.5 are shown, see the workflow run for all results.

[headtr1ck]

I opened an upstream issue for that. Maybe it's better to fix this at the source instead of ugly custom builds here. Thoughts?

[dcherian]

On the whole, I think this is worth it. We maintain vendored functions for efficiency and backcompat already. It would definitely be good to upstream. In the meantime, we could stick these functions in scipy_compat.py for example, following current code conventions.

I'll leave it to your judgement.

If you'd like to merge, please add some property tests comparing the output to scipy for random inputs. Another property could be the interp of a 1D vector shoul dmatch the interp of that vector broadcasted out to nD.

[headtr1ck]

It seems that scipy is very interested in getting this patch. I might work on that, so maybe this PR is not necessary. I don't think we really need a compatiblity layer in that case. If you want it faster, update your scipy.

[dcherian]

Great. As an aside, I loved that thread. Great example of how LLMs can help out OSS!