# Multivariate function interpolation

I have a (nonlinear) function which takes as input 4 parameters and produces a real number as output. It is quite complex to compute the function value given a set of parameters (as it requires a very big summation).

I'd like to answer queries on this function efficiently so I was thinking of trying to use some interpolation methods. I have used Chebyshev polynomials to interpolate univariate functions, but I haven't been able to find (or understand) anything on interpolating multivariate functions. I'm not set on using Chebyshev polynomials, I have just had some exposure to them and know they tend to be efficient (in terms of # of necessary coefficients and interpolation error).

I was wondering if anyone could give me (an engineer) any pointers for how to go about interpolating a multi-variate function? Simple examples or sample code would be awesome, but I'll take any attempts to explain how interpolation would work in higher dimensions, including (readable) references.

• Have you seen this? One possibility you can try is to take suitably scaled Chebyshev nodes in each variable (essentially, a tensor product) and then perform a method for multidimensional polynomial interpolation, like in the paper I linked to. – J. M. is a poor mathematician Dec 16 '11 at 4:01
• @J.M. Thanks for the reference, I'll check it out. – Nick Dec 16 '11 at 17:09

Since nobody's answered yet, I throw this out there (I'm no expert in multivariate interpolation, hopefully someone with expertise will eventually weigh in).

I imagine you've already looked at the Wikipedia Article on Multivariate Interpolation. There's a lot of stuff out there. However, if you just want a "quick fix", I'd go with some sort of linear regression.

(1) Pick your favorite template function, say, $f(x,y,z,t) = Ax^2t+B\sin(y+z)+Cxyz$

(2) Plug-in some "known" values, so something like:

• $w_1 = f(x_1,y_1,z_1,t_1) = Ax_1^2t_1+B\sin(y_1+z_1)+Cx_1y_1z_1$
• $w_2 = f(x_2,y_2,z_2,t_2) = Ax_2^2t_2+B\sin(y_2+z_2)+Cx_2y_2z_2$
• $w_3 = f(x_3,y_3,z_3,t_3) = Ax_3^2t_3+B\sin(y_3+z_3)+Cx_3y_3z_3$
• $w_4 = f(x_4,y_4,z_4,t_4) = Ax_4^2t_4+B\sin(y_4+z_4)+Cx_4y_4z_4$

(3) Translate to a linear system:

$${\bf w} = \begin{bmatrix} w_1 \\ w_2 \\ w_3 \\ w_4 \end{bmatrix} = \begin{bmatrix} x_1^2t_1 & \sin(y_1+z_1) & x_1y_1z_1 \\ x_2^2t_2 & \sin(y_2+z_2) & x_2y_2z_2 \\ x_3^2t_3 & \sin(y_3+z_3) & x_3y_3z_3 \\ x_4^2t_4 & \sin(y_4+z_4) & x_4y_4z_4 \end{bmatrix} \begin{bmatrix} A \\ B \\ C \end{bmatrix} = M {\bf x}$$

(4) Solve the Normal Equations: $M^T M {\bf x} = M^T {\bf w}$ to find the solution which best fits your data. Of course, if you pick a general enough function $f$ to begin with, the system from (3) will be consistent and you'll be able to solve ${\bf w} = M{\bf x}$ directly.

I imagine finding error bounds is quite daunting. But if you used a general multivariate polynomial of fairly high degree, is ought to do the trick. Even something like: $f(x,y,z,t) = Ax^3+By^3+Cz^3+Dt^3+Ex^2y+Fxy^2+Gxyz+\cdots+Hx+Iy+Jz+K$ ought to do a decent job.

• This is interesting, thank you for the thoughtful response. I'm wondering how you came up with your 'template' function? and in general how does one identify a 'good' template? Please excuse my naivety, but shouldn't the templates be a sum of orthogonal functions to guarantee we can approximate the input function? – Nick Dec 16 '11 at 17:07
• My template function is known to experts as "Bill's Function" :) Ok. Just kidding, I made it up as a random example. For general functions, I would use a multivariate polynomial (they are versatile and simple to compute with). On the other hand, if your function is periodic, then sines and cosines might be better. – Bill Cook Dec 16 '11 at 22:46

I would suggest taking a look at Moving Least Squares (MLS), which can generate quite smooth interpolators using low-order functions. One benefit of MLS is that it can easily work with scattered data, which means you can sample your expensive function at relevant locations only. If you want to perform uniform sampling of your function you can look at the simplicial lattice described in [1], which is a nice alternative to the standard regular grid.

[1] Fast High-Dimensional Filtering Using the Permutohedral Lattice. Adams et al. Eurographics 2010.

• Thanks for the reference. Is MLS an approach for function interpolation or regression? I'm wondering what the difference between these methods would be in terms of the approximation error? – Nick Dec 16 '11 at 17:14
• MLS uses a weight function to define the influence of each data point (based on distance to the query point). If you chose a weight function which has a singularity at zero you will get a perfect interpolation of the data points. If you chose a weight function with no such singularity, you will get some degree of outlier-robustness (ie, more like a regression). If the weight function is constant for all points MLS becomes standard Least Squares. – Eduardo Dec 20 '11 at 2:20
• Sorry but I cannot help you with any error analysis on MLS. You could probably find some references somewhere. Also your best bet is gathering some empirical error data for your specific application. – Eduardo Dec 20 '11 at 2:24