4
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ It's important to note that the size of the system you must solve for the least squares fitting depends on the number of coeficients, not the number of data points. In your example for option 1 this would give a 10x10 matrix which should be quite fast to solve. (Fast as in tens of thousands of times per second.) $\endgroup$
    – Doug Lipinski
    Mar 2, 2016 at 18:50
  • 1
    $\begingroup$ You say "interpolation", but it seems you actually want curve fitting then function evaluation of the fitted curve. Is that correct? $\endgroup$
    – Doug Lipinski
    Mar 2, 2016 at 18:50
  • $\begingroup$ I meant curve fitting via least squares, sorry. So, with regards to the 10x10 vs. N by 10, I think I was making a huge dumb mistake. I set up my system as a Vandermonde matrix where each row represented $1, x_i, y_i, x_i^2, x_i*y_i, ... y_i^3 $ and had one row per data point $(x_i,y_i)$. I then let the linear solver find the least squares solution. Instead, I could analytically write out the equations that minimize the $\ell^2$ norm of $z_i - p(x_i,y_i)$, which leads to a 10 by 10 (for 10 poly coefficients). That should be a huge help $\endgroup$
    – user35959
    Mar 3, 2016 at 0:40
  • $\begingroup$ The easiest thing to do is to form your matrix $M_{\text{N}\times 10}$ for the overdetermined system and then solve the system $(M^\intercal M)x = M^\intercal b$ which gives the least squares solution. See en.wikipedia.org/wiki/… I'm not sure what solver you're using, but a well written code that's capable of solving linear least squares problems will already do this for you. For example, MATLAB's \ aka mldivide should do this very efficiently. $\endgroup$
    – Doug Lipinski
    Mar 3, 2016 at 13:39

1 Answer 1

Reset to default
2
$\begingroup$

First of all, you can turn the fitting problem for 1 and 3 into an iterative optimization problem if you don't want to deal directly with solving a large system of equations using a matrix approach.

For 2, you can use a hash table and design a spatial hash function to map a coordinate to the appropriate subregion, which would be very fast.

For 3, given you use RBFs that are locally impacting the fit to the data, you could use a KD Tree with an m-Nearest Neighbors algorithm to find the RBFs that impact a given input point and only use those to compute the value. I have found this as very efficient in my implementations. The tough point is figuring out how many RBFs to use and how large to set their Region of Influence to be.

Other Idea

If you don't have storage requirements, you could also consider doing Locally Weighted Regression using a local linear or quadratic fit. I have built codes to do this on large data sets with 9-D data points and it was quite fast.

The benefit to this algorithm is you don't have to fit any models beforehand. You end up using a set of data close to your input data point that you want to find the function value at, and do a local fit using a simpler model (linear, quadratic), so it's fast. If you use a KD Tree to find the M Nearest Neighbors, it should be quite quick.

Final Thoughts

I am suspecting option 2 might be your best balance between speed and reduced data storage, especially if you use a hash table/hash function pair to map the input points to the appropriate subregion.

Edit

I almost forgot, but for method 1, you could also incorporate regularization on the coefficients of your polynomial model. This will help you avoid Runge's Phenomenon and should give you the fastest model to evaluate lots of points relative to the other approaches.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for the input! I appreciate it. I like the idea of the hash table for Option 2 as well as a regularization for Option 1. I'll test some of these ideas out later. One other issue I didn't make clear in my original post is that I'm not just fitting one dataset (i.e., doing one polynomial model). I will have constantly new incoming datasets with demands for fitting and then subsequent evaluations. Hence my concern over the startup time (e.g., solving the least squares polynomial system for Option 1 ). $\endgroup$
    – user35959
    Feb 27, 2016 at 18:49
  • 1
    $\begingroup$ @user35959 I think you should better express your requirements for this task in your question. What is the number of function evaluations you expect to make for a given fit of the data? Is this functionality supposed to run in real time or is this more a task that can take weeks to run and you want to shrink that as much as possible? Is the function evaluations expected to be the biggest bottleneck? $\endgroup$
    – spektr
    Feb 27, 2016 at 19:09
  • $\begingroup$ Sorry for the lack of clarity; I edited the question. This goal is an as near-real-time thing as possible, with multiple datasets needing fitting and many, probably hundreds at least, of function evaluations per fit dataset. $\endgroup$
    – user35959
    Feb 27, 2016 at 20:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.