A Krylov subspace algorithm for multiquadric interpolation in many dimensions

We consider the version of multiquadric interpolation where the interpolation conditions are the equations s((x) under bar (i)) f(i) = 1, 2, ..., n, and where the interpolant has the form s((x) under bar) = Sigma(j=1)(n) lambda(j)(parallel to(x) under bar - (x) under bar (j)parallel to(2) + c(2))(1/2) + alpha, (x) under bar is an element of R(d), subject to the constraint Sigma(j=1)(n) lambda(j) = 0. The points (x) under bar (i) is an element of R(d), the right-hand sides f(i), i = 1, 2, ..., n, and the constant c are data. The equations and the constraint define the parameters lambda(j), j = 1, 2, ..., n, and alpha. The resultant approximation s approximate to f is useful in many applications, but the calculation of the parameters by direct methods requires O(n(3)) operations, and n may be large. Therefore iterative procedures for this calculation have been studied at Cambridge since 1993, the main task of each iteration being the computation of s ((x) under bar (i)), i = 1, 2,..., n, for trial values of the required parameters. These procedures are based on approximations to Lagrange functions, and often they perform very well. For example, ten iterations usually provide enough accuracy in the case d = 2 and c = 0, for general positions of the data points, but the efficiency deteriorates if d and c are increased. Convergence can be guaranteed by the inclusion of a Krylov subspace technique that employs the native semi-norm of multiquadric functions. An algorithm of this kind is specified, its convergence is proved, and careful attention is given to the choice of the operator that defines the Krylov subspace, which is analogous to pre-conditioning in the conjugate gradient method. Finally, some numerical results are presented and discussed, for values of d and n from the intervals [2, 40] and [200, 10 000], respectively.