[SciPy-Dev] ENH: improve RBF interpolation

[Stefan van der Walt]

Hi Trever,

On Fri, Jan 29, 2021, at 05:13, Trever Hines wrote:
> I would like to contribute code to scipy to address some issues regarding RBF interpolation. The code can be found on my branch _here_ <https://github.com/treverhines/scipy/blob/improve-rbf-interpolation/scipy/interpolate/rbfinterp.py>. My contribution consists of two new classes for scattered N-D interpolation:


This is fantastic; thank you so much for sharing your expertise on RBFs!

>  1. `RBFInterpolator`: This is intended to be a replacement for `Rbf` that addresses issues mentioned in 9904 <https://github.com/scipy/scipy/issues/9904> and  4790 <https://github.com/scipy/scipy/issues/4790>. Namely, the major differences with `Rbf` are 1) the usage is similar to `NearestNDInterpolator` and `LinearNDInterpolator` making it easier to swap out different interpolation methods, 2) the sign of the smoothing parameter is correct (see page 10 of these lecture notes <http://pages.stat.wisc.edu/~wahba/stat860public/lect/lect8/lect8.pdf>), and 3) the interpolant includes polynomial terms.  
> For some RBF choices (values of “linear”, “thin_plate”, “cubic”, “quintic”, or “multiquadratic” for `function` in `Rbf`), the additional polynomial terms are needed to ensure that the interpolation problem is well-posed (see theorem 3.2.7 in this document <http://amadeus.math.iit.edu/~fass/603_ch3.pdf>). Without the additional polynomial terms for these RBFs, I have noticed that some values for the smoothing parameter (with the corrected sign) result in an obviously erroneous interpolant. Even when the chosen RBF does not require additional polynomial terms, they still can improve the quality of the interpolant. In particular, the polynomial terms are able to accommodate shifts or linear trends in data, which the RBFs tend to struggle with by themselves.

Is there any advantage to keeping the old interface, or should this eventually replace Rbf entirely?

>  1. `KNearestRBFInterpolator`: This class performs RBF interpolation using only the k neaest data points to each interpolation point (which was suggested in 5180 <https://github.com/scipy/scipy/issues/5180>). This class is useful when there are too many observations for `RBFInterpolator` (on the order of tens of thousands) and you want an interpolant that* looks** *smoother than what you get with `NearestNDInterpolator` or `LinearNDInterpolator`.  
> My concern with interpolation using the k nearest neighbors is that it is a bit of an ad hoc strategy to work around computational limitations. That being said, I have seen a similar strategy used in the Kriging <https://pro.arcgis.com/en/pro-app/latest/help/analysis/geostatistical-analyst/search-neighborhoods.htm> world (Kriging is a form of RBF interpolation). 

Superb!  I've been trying to do RBF interpolation with N>1000 and with the N^2 memory requirement it gets you pretty quickly.  This makes the algorithm much more pragmatic to apply to, e.g., images.  We can always add other picking strategies later on.

> I wou;d appreciate your feedback on whether you think these would be valuable contributions to scipy.  If so, I will make the pull request after adding benchmarks, unit tests, and more docs.


I'd say 100% yes; thank you again.

Best regards,
Stéfan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/scipy-dev/attachments/20210204/f4928e0d/attachment.html>