The problem of interpolating amongst data defined on a rectangular mesh in real N-space can be conveniently solved by constructing suitable tensor product hypersurfaces. In particular the problem can be neatly formulated with the use of Kronecker products of matrices. Greville (1961) has observed that the concept of the pseudo-inverse of a matrix allows this formulation to be extended to the corresponding least-squares problem. The purpose of this note is to review this approach and to relate it to the method of surface fitting described by Clenshaw and Hayes (1965), in which the surface is obtained by two stages of curve fitting. There is a dearth of algorithms for solving the N-variable problem and two new ones are given here in outline. They are similar to those developed independently by Pereyra and Scherer (1973).