I recently faced the problem of efficiently approximating a very large set of data points and, neither having a model of the empiric function, nor of the error distribution, my method of choice would have been a least squares fit with B-Splines.
The resulting system of equations would however have been forbiddingly large and a peculiarity of my problem was, that the value of the approximating function had only to be determined for a small set of clustered arguments, which means, that there were large intervals, where the function would never have to be evaluated.
My idea was then to approximate the sets of data points, that are inside the support of the B-Spline basis function, by an individual polynomial.
The evaluation of empiric function would then be done via a convex combination of the approximating polynomials that are associated to those supports of the B-Spline basis functions, that contain the argument point.
The weights for the convex combination would be obtained from an evaluation of the B-Spline basis functions for the argument.
The advantages of the method would be that there is no need to completely determine the approximating function if only a few function values are needed and, the effect of outliers is limited to a small neighborhood, which could be desirable, if an approximation of shape is more important than error minimization (e.g. in case of a local deformation due to damage).
Has such a kind of approximation already been described and, where can I find information about it?