Fix interp performance regression for vectorized interpolation

[headtr1ck]

Fixes the performance regression from #9881 when the new coordinates vary along a dimension of the interpolated object, e.g. da[t, r, z].interp(z=new_z[t]).

• New apply_ufunc(..., vectorize="shared"): only loops over the broadcast dimensions along which more than one input varies, all others are passed to func in bulk (which must then broadcast its inputs following numpy rules).
• interp uses it, so it only loops over t instead of t and r like np.vectorize. Works for numpy and dask alike, without changing chunks or the task graph.
• _localize and the datetime conversion run once per call before apply_ufunc instead of once per loop iteration.
• Fix interp with empty new coordinates for method="linear" / "nearest".
• New asv benchmarks InterpolationVectorized / InterpolationVectorizedDask.
• Docs: apply_ufunc docstring examples from basic to advanced, and a vectorize="shared" section in the apply_ufunc example notebook, which is executed in the docs build again.

vectorize="shared" could also enable multidimensional weights in polyfit (#9846): loop over the dimensions the weights vary along and solve all other series in one lstsq call.

This is an alternative to #11496, which declares the free dimensions as core dimensions (numpy only).

benchmark	before #9881	main	#11496	this PR
numpy, linear	1.7 ms	1173 ms	25 ms	9.0 ms
numpy, cubic	110 ms	605 ms	17 ms	18.1 ms
dask, linear	11.8 ms	1271 ms	1265 ms	30.7 ms
dask, cubic	102 ms	789 ms	807 ms	52.5 ms

The remaining gap for linear compared to before #9881 comes from the Python loop over t: each iteration builds its own scipy interpolator. Before #9881, the shared dimensions were added as extra interpolation dimensions, which needs merged dask chunks, only works for n-D capable methods with numeric coordinates, and leaks NaNs between slices.

• Closes interp() performance regression in 2025.1.0+ #10683
• 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]

It's somehow much uglier than just vectorize=True but it seems to work. Any ideas how to make this better understandable?

[dcherian]

I can take a look later. My conclusion a year ago was that there are no real nice solutions without reinventing parts of apply_ufunc. Perhaps there's a vectorization helper for apply-ufunc that we can use here?

But this is why i accepted the slowdown for #9881; indeed it was much worse than expected though.

[dcherian]

: each iteration builds its own scipy interpolator.

we can totally construct an interpolator outside the loop no?

[headtr1ck]

: each iteration builds its own scipy interpolator.

we can totally construct an interpolator outside the loop no?

I let Claude benchmark several variations of that and neither showed a clear performance improvement. There are always tradeoffs for small or large arrays while just adding more code boilerplate.

[headtr1ck]

I can take a look later. My conclusion a year ago was that there are no real nice solutions without reinventing parts of apply_ufunc. Perhaps there's a vectorization helper for apply-ufunc that we can use here?

But this is why i accepted the slowdown for #9881; indeed it was much worse than expected though.

Thanks!

I couln't find any such helper, but thinking about it we could totally turn the new interp code into such a helper. Potentially even add it as a new vectorized="shared" option.

[headtr1ck]

@dcherian I turned the loop into a helper and added it to apply_ufunc as vectorize="shared": it only loops over the broadcast dimensions along which more than one input varies and passes all other broadcast dimensions to func in bulk. interp now just uses apply_ufunc(..., vectorize="shared"), so _interpnd is almost back to its version on main. Performance is the same or slightly better than before, see the updated description.

[This is Claude Code on behalf of Michael Niklas]

[github-actions]

⚡ Zoom zoom! Some benchmarks improved.

Change	Before [4929aa9]	After [4c634f3]	Ratio	Benchmark (Parameter)
🚀	562±4ms	20.5±0.2ms	0.04	interp.InterpolationVectorized.time_interpolation_vectorized('cubic')
🚀	1.22±0s	8.53±0.04ms	0.01	interp.InterpolationVectorized.time_interpolation_vectorized('linear')
🚀	611±2ms	65.7±1ms	0.11	interp.InterpolationVectorizedDask.time_interpolation_vectorized('cubic')
🚀	1.50±0.02s	38.2±1ms	0.03	interp.InterpolationVectorizedDask.time_interpolation_vectorized('linear')

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

Comparing 4c634f3b (this PR merged into the base branch) with 4929aa90 (base branch).
Only changes by more than a factor of 1.5 are shown, see the workflow run for all results.

[headtr1ck]

@dcherian do you like the new vectorize option? In xarray I didn't find any other usecase for it (most applications are differing only by none or a single dim and there np.vectorize works).
But I guess some very advanced users might be able to utilize it.