It may help to first look at the easiest non-trivial case, i. e. n==2. The function you are looking at maps 2-space to a real number, and if you combine the result y with the argument vector (x[0], x[1]), you get a point in 3-space, (x[0], x[1], y). What you have, now, is a surface. For n>2 you'd get an n-dimensional hypersurface in (n+1)-space instead.
If you want to approximate a surface (i. e. for n==2), there are plenty of spline interpolation and approximation algorithms around, but if the points you know (the yi values) are arbitrarily spread over the 2-space of the xi, then it's a rather complex matter!
Depending on the number of points you know, you can come up with various methods that make some sense: e.g. if you start with exactly n+1 known points, these points will define a hyperplane (for n==2, you need three points to define a plane), and you can simply project any argument vector x onto that plane to get the corresponding y value. If you have only n or less known points, you'd similarly define a linear subspace and project your X onto that subspace.
If you have more than n+1 known points however, and they are not co-(hyper)planar, then not only do the points describe a nonlinear subspace, there may not even be a unique 'best' fit solution to the manifold they describe. In that case, I am not aware of any generic algorithm. For n==1 you might simply order the xi by size, determine the interval [xk,xk+1] that x lies in, and then interpolate y between yk and yk+1. For n=2 the corresponding algorithm would be more complex, since you'd need to figure out a triangulation of the plane based on the xi, check which triangle the argument vector x lies in, and then interpolate y using that triangle. The tricky part here is how to define the triangulation. For higher values of n, the same idea may work, but the task of subdividing the hyperplane of Xs into n-dimensional simplexes gets even more complex an ambiguous!
In any case, this approach would only deliver a linear approximation of a non-linear surface. A more generic approach would be to calculate a weighted sum of the yi, where the weights should diminish based on the distance of the argument vector X from the correponding Xis. The tricky part is to somehow balance the weights so that they always add up to a total of one. At the moment I have no idea how to achieve that for n>2.