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.
 A: 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.
A: 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. 
