sparse_interp_nd, an Octave code which constructs a sparse interpolant to a function f(x) of a multidimensional argument x.
The problem is defined as follows: for any point x in some M-dimensional product space X, such as the unit hypercube, we are able to provide the value f(x) for some function f(). We wish to determine a new function g(), (the interpolant) which will match the value of f() on some specified set of points {xi}. Presumably, we choose the number and location of the points {xi} to ensure that the interpolating function g() is a good fit to f().
There are standard techniques for producing a family of interpolants g(j)(), such that, as we increase the index j, the interpolating process will exactly match any function f() which happens to be a polynomial of limited degree. The usual techniques for doing this are derived from tensor product interpolants, whose expense grows exponentially with the spatial dimension.
For simplicity, we'll assume that the argument space X is the unit hypercube, and that we have a family of 1D interpolation rules g(j)(), where the index j is an indication of the precision of the interpolation rule. A typical tensor product interpolant G(J)() can then be thought of as constructed by the product
G(J)(X) = g(j1)(x1) * g(j2)(x2) * ... * g(jm)(xm).
where, of course, J = (j1,j2,...jm) and X = (x1,x2,...,xm).
A procedure due to Smolyak, which is more typically applied to problems in quadrature, can also be used for the interpolation problem. The Smolyak interpolation rule of level L is defined by
A(L,M) = sum ( L-M+1 <= |J| <= L ) C(|J|) * g(j1)(x1) * g(j2)(x2) * ... * g(jm)(xm).
Here |J| = j1+j2+...+jm, and, for each |J|, the sum must be taken over all
possible vectors J with nonnegative integer entries that sum to |J|.
Thus, a naive implementation of a sparse interpolant would, for a given spatial dimension M, pick a level L, determine the (at most L+1) coefficients C(), construct every tensor product interpolant of total degree L or less, evaluate each interpolant at the point X, and combine these values with the appropriate weights C() to arrive at the sparse grid interpolant value at X.
Some improvements to this approach can be suggested. First, many of the coefficients C() may be zero, because the coefficient vector C for an M-dimensional sparse interpolant of level L will have at most M nonzero coefficients. Secondly, if the 1D interpolation family is chosen so that the interpolant points of successive members are nested, then it is possible to simplify the evaluation process greatly.
The code requires access to the R8LIB library.
The computer code and data files described and made available on this web page are distributed under the MIT license
sparse_interp_nd is available in a C version and a C++ version and a Fortran90 version and a MATLAB version and an Octave version.
lagrange_interp_nd, an Octave code which defines and evaluates the lagrange polynomial p(x) which interpolates a set of data depending on a multidimensional argument x that was evaluated on a product grid, so that p(x(i)) = z(i).
r8lib, an Octave code which contains many utility routines using double precision real (r8) arithmetic.
rbf_interp_nd, an Octave code which defines and evaluates radial basis function (RBF) interpolants to multidimensional data.
shepard_interp_nd, an Octave code which defines and evaluates Shepard interpolants to multidimensional data, based on inverse distance weighting.
spinterp, an Octave code which carries out piecewise multilinear hierarchical sparse grid interpolation; an earlier version of this software is ACM TOMS algorithm 847, by Andreas Klimke;
test_interp_nd, an Octave code which defines test problems for interpolation of data z(x), depending on an m-dimensional argument.