Thin Plate Splines are a great method to find a smooth interpolating surface given scattered data. Essentially, the method involves calculating weights for a radial basis function centred around each of the data points:
$$f({\bf{x}}) = a + bx + cy + \sum_{i=1}^N{w_i \phi(|{\bf{x}} – {\bf{x}}_i|)}$$
$$\phi(r) = r^2 \log r$$
However, there are two problems with this method from a practical perspective:
- Non locality – Far away points can have large influence on the functions value, and in addition, all points are needed to reconstruct a given value.
- Numerical instability – The calculation for each point involves summation over $N$ points, which may be large, leading to accumulation of small errors. Additionally, solving the matrix equation for large $N$ may be unstable as well.
Are there other methods or extensions to the TPS method that smoothly solve the 2D interpolation problem from scattered data that don't have the above issues?
Best Answer
Far away points in their mass, not individually, may have effect on filling in the "holes" of corresponding scale. In the presence of enough neighboring points the influence of far trends will be negligible. As a demonstration, write a univariate linear or cubic spline as a N-sum of RBF formally including "far-away" terms.
The evaluation of large RBF could be done with the method of fast multipoles (FMM) to a given precision. The idea is to pre-compute all the influences of hierarchical groups of far-away points. This is similar to multipole decomposition in electrostatics: far away charges are approximated as one net charge or dipole, a finer structure of remote cloud doesn't matter here.
For moderate RBF one could indeed sum over N, perhaps using Kahan summation, for example, to preserve floating-point precision as much as possible.
Fitting a large globally supported RBF interpolant is indeed not an easy problem. One current approach is to use various matrix preconditioning methods. I for example am developing with my mentor an alternative approach, a greedy algorithm, which should allow fitting to very large oversampled datasets.
If your data is large but have more or less regular sparsity, or you can estimate the scales of holes, I would recommend compactly supported RBF or a series of them with different radius, that is how it's done I think in the ALGLIB library, for example.
"...if the data are much further apart than the support size of the compactly supported radial basis functions, we will get a useless approximation to the data, although it always interpolates... Conversely, if the data are far closer to each other on average than the support size of the ra- dial basis function, we lose the benefits of compact support. Thus, the scaling should be related to the local spacing of the centres." -- Buhman, Radial Basis Functions. See also here.
If you don't have an estimate on the diameter of the largest possible gap and can't use a uniform radius, you may want to look for variable scaling schemes, for example: "This algorithm builds multiscale (hierarchical) RBF models... with decreasing radii... Values predicted by the first layer of the RBF model are subtracted from the function values at nodes ... residuals are passed to the next iteration."
That is too short and categorical to live by. A cubic spline minimizes a curvature-like functional while can be written in a locally supported B-spline basis.
Even though CSRBF may look visually unpleasant when the radius is too small, they still minimize a norm, i.e., are the best approximator in a so-called native functional space which for Wendland's CSRBF of minimal degree is a classical Sobolev space. This is Theorem 10.35 in Wendland, Scattered Data Approximation.