Sparse multivariate polynomial interpolation via the quotient-difference algorithm

To reconstruct a black box multivariate sparse polynomial from its floating point evaluations, the existing algorithms need to know upper bounds for both the number of terms in the polynomial and the partial degree in each of the variables [2]. On the other hand, Rutishauser's quotient-difference algorithm [7, 3], or the qd-algorithm, is an iterative method that can be use to determine the poles of a meromorphic function from its Taylor coefficients. Combining the relation between the qd-algorithm, Prony's method and their connections to the Ben-Or/Tiwari algorithm [6], we present a new algorithm for reconstructing sparse multivariate polynomials in floating point arithmetic, in which neither the number of terms nor any partial degrees are required in the input [1]. For example, consider the bivariate black box polynomial function, <i>f</i>(<i>x</i><sub>1</sub>,<i>x</i><sub>2</sub>) = 1/3<i>x</i>6/1+4<i>x</i>4/2--2.1<i>x</i>4/1--4<i>x</i>2/2+<i>x</i><sub>1</sub><i>x</i><sub>2</sub>+4<i>x</i>2/1, termed the six-hump camel back function. Evaluate f (x1 , x2 ) at the non-equidistant s-th powers of (1/3, 1/5), in which 3 and 5 are relatively prime, and proceed with the qd-algorithm on the data <i>f</i> (1/3, 1/5), <i>f</i> (1/3<sup>2</sup>, 1/5<sup>2</sup>), <i>f</i> (1/3<sup>3</sup>, 1/5<sup>3</sup>), ....The magnitude of the values in the first five <i>e</i>-columns drops from 10<sup>--2</sup> to machine precision, while all values in the sixth <i>e</i>-column are of the order of machine precision. We conclude that the number of terms is 6 and obtain the multi-indices in the support directly from the <i>q</i>-values: <ul><li>1/q<sub>1</sub><sup>(S)</sup> → 9 = 3<sup>2</sup>.</li> <li>1/q<sub>2</sub><sup>(S)</sup> → 25 = 5<sup>2</sup>.</li> <li>1/q<sub>3</sub><sup>(S)</sup> → 625 = 5<sup>4</sup>.</li> <li>1/q<sub>2</sub><sup>(S)</sup> → 15 = 3<sup>1</sup>5<sup>1</sup>.</li> <li>1/q<sub>1</sub><sup>(S)</sup> → 81 = 3<sup>4</sup>.</li> <li>1/q<sub>1</sub><sup>(S)</sup> → 729 = 3<sup>6</sup>.</li></ul> After all six multi-indices in <i>f</i> are known, the coefficients in <i>f</i> can be recovered by a linear system formed by (at least six) evaluations of <i>f</i>. The convergence to zero of some e-columns in the qd-algorithm can be interpreted in an analogous way as the zero discrepancies in the early termination strategy [5, 4]. This is illustrated in Maple. Our sparse interpolation can be extended to polynomials represented in certain non-standard bases [4]. We show an extension to sparse interpolation of polynomials in the Pochhammer basis and present our current research on sparse polynomial interpolation in the Chebyshev basis, both in floating point arithmetic.