This may be a trivial question, and I apologize if so. Consider the following simple problem: We have a 2D, regular grid of points (say $X = [0,5000] \times [0,5000]$) spaced uniformly by units of 1 (e.g., $(0,0), (0,1), (1,0) \in X$ and there exists a function $f(x,y)$ that we aim to approximate on $X$. Let's assume it's nice and smooth, at least $C^2$. Assume that we have data samples on a uniform subgrid (say every 100 points in a direction we have a sample value $f(x_i,y_i)$). We all know there exists many different methods to approximate such a function. What I'm curious about is suggestions on which methods provide me the fastest evaluation of a new given point not on the subgrid (e.g., at a point $(50,25)$).
Option 1: Multivariate polynomial interpolation.
We could choose a multivariate polynomial basis of some fixed degree (let's assume 3rd order) so , $p(x,y) = c_0 + c_1 x + c_2y + c_3 x^2 + c_4 xy + c_5 y^2 + c_6 x^3 + c_7 x^2 y + c_8 xy^2 + c_n y^3$.
There are a number of advantages to this method: storage is small (we just need to store the coefficients, which is a very small number relative to the number of points in the grid) and evaluation on a new point is relatively fast (we have several multiplications, plus a few additions).
However, we also have to fit this polynomial. This requires us to solve a least squares problem, where we will have a matrix of $N \times d$ where d is the number of coefficients and $N$ is the number of sample points. I believe my best option is an SVD solve (using maybe Armadillo or Eigen). This will definitely be the slowest part of the whole ordeal, so there is significant up front cost.
The next disadvantage is issues such as Runge's Phenomenon. I am usually quite cautious about fitting with multivariate polynomials due to concerns such as poor accuracy in certain regions. Furthermore, the choice of basis significantly affects the solvability of the system. If this were 1D, I'd do Chebyschev polynomials if possible. Is there a better basis for such a problem in 2D?
Option 2: Piecewise Polynomials on the square grids
I'm a fan of this option. I can use linear or quadratic elements on the actual subregions (of say 100 by 100 units). The disadvantage is that I'm required to store the $N$ known function samples. So this requires more storage, so copying and moving around this information is more significant than moving around for example the polynomial coefficients (which were a fixed, tiny number). Evaluation should be quite quick. I need to know the values at my vertices of the cell, and given a new point (x,y), I need to compute which square cell it's in. I assume I can do this by dividing $x$ and $y$ by the stride (say 100) and flooring it to find out the bottom left vertex of the cell the new point is in. The huge advantage i see here is that I'm not solving large systems to get coefficients for an approximating function.
The only disadvantage I really see is the case where my grid spacing varies (say 100 by 100 cells for a bit, then 50 by 50 cells). Given a new point to evaluate a function value at, (x,y), I need to efficiently and quickly compute which cell (x,y) is in, then do my piecewise polynomial evaluation in that cell quickly. Maybe I'm being dumb, but what is a good efficient method to quickly compute which cell a point lies in when the cells are not all the same size?
3: Radial Basis Function methods
This requires solving large systems and lots of coefficients, so I'm not seeing an advantage here. If I had irregularly sampled points, I would consider this.
Summary: Speed is of the essence for me; I will have lots of incoming points and evaluations of the function value at these points desired and I need to also as quickly as possible set up the approximating function (e.g. solving for polynomial or RBF coefficients). Suggestions would be great. Thanks!
EDIT: I was unclear in my question. I am aiming for real-time (or as near as possible) fitting of an incoming dataset (as described above) with some method (polynomials, piecewise polys, whatever), which will then subsequently, in real-time receive many requests for evaluations. I will receive more than just one dataset, possibly many per second (tens? hundreds?). Hence, the upfront cost of building the interpolation/fitting scheme and subsequent evaluations both needs to be fast. I will definitely need to parallelize as much as possible.