Far away points can have large influence...
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.
...all points are needed to reconstruct a given value. The calculation for each point involves summation over N points, which may be large, leading to accumulation of small errors.
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.
Additionally, solving the matrix equation for large N may be unstable as well
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.
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?
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.
EDIT in reply to the comments
"...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."
RE: "Interpolating, smooth, local. Pick any two"
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.
RE: "minimizing some roughness functional ... you cannot minimize it locally, even if your basis support is local"
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.