A measure of the above ‘distance’ between the set of data points and the function $f\left(x\right)$ is needed. The distance from a single data point $\left(x_r,y_r\right)$ to the function can simply be taken as and is called the residual of the point. (With this definition, the residual is regarded as a function of the coefficients contained in $f\left(x\right)$; however, the term is also used to mean the particular value of ${\epsilon }_r$ which corresponds to the fitted values of the coefficients.) However, we need a measure of distance for the set of data points as a whole. Three different measures are used in the different methods (which measure to select, according to circumstances, is discussed later in this sub-section). With ${\epsilon }_r$ defined in (1), these measures, or norms, are and respectively the ${\ell }_1$ norm, the ${\ell }_2$ norm and the ${\ell }_{\infty }$ norm.

Minimization of one or other of these norms usually provides the fitting criterion, the minimization being carried out with respect to the coefficients in the mathematical form used for $f\left(x\right)$: with respect to the $b_i$ for example if the mathematical form is the power series in (8) below. The fit which results from minimizing (2) is known as the ${\ell }_1$ fit, or the fit in the ${\ell }_1$ norm: that which results from minimizing (3) is the ${\ell }_2$ fit, the well-known least squares fit (minimizing (3) is equivalent to minimizing the square of (3), i.e., the sum of squares of residuals, and it is the latter which is used in practice), and that from minimizing (4) is the ${\ell }_{\infty }$, or minimax, fit.

Strictly speaking, implicit in the use of the above norms are the statistical assumptions that the random errors in the $y_r$ are independent of one another and that any errors in the $x_r$ are negligible by comparison. From this point of view, the use of the ${\ell }_2$ norm is appropriate when the random errors in the $y_r$ have a Normal distribution, and the ${\ell }_{\infty }$ norm is appropriate when they have a rectangular distribution, as when fitting a table of values rounded to a fixed number of decimal places. The ${\ell }_1$ norm is appropriate when the error distribution has its frequency function proportional to the negative exponential of the modulus of the normalized error – not a common situation.

However, it may be that you are indifferent to these statistical considerations, and simply seek a fit which can be assessed by inspection, perhaps visually from a graph of the results. In this event, the ${\ell }_1$ norm is particularly appropriate when the data are thought to contain some ‘wild’ points (since fitting in this norm tends to be unaffected by the presence of a small number of such points), though of course in simple situations you may prefer to identify and reject these points. The ${\ell }_{\infty }$ norm should be used only when the maximum residual is of particular concern, as may be the case for example when the data values have been obtained by accurate computation, as of a mathematical function. Generally, however, a method based on least squares should be preferred, as being computationally faster and usually providing more information on which to assess the results. In many problems the three fits will not differ significantly for practical purposes.

Some of the methods based on the ${\ell }_2$ norm do not minimize the norm itself but instead minimize some (intuitively acceptable) measure of smoothness subject to the norm being less than a user-specified threshold. These methods fit with cubic or bicubic splines (see (10) and (14) below) and the smoothing measures relate to the size of the discontinuities in their third derivatives. A much more automatic fitting procedure follows from this approach.

The use of the above norms also assumes that the data values $y_r$ are of equal (absolute) accuracy. Some of the methods enable an allowance to be made to take account of differing accuracies. The allowance takes the form of ‘weights’ applied to the $y$ values so that those values known to be more accurate have a greater influence on the fit than others. These weights, to be supplied by you, should be calculated from estimates of the absolute accuracies of the $y$ values, these estimates being expressed as standard deviations, probable errors or some other measure which has the same dimensions as $y$. Specifically, for each $y_r$ the corresponding weight $w_r$ should be inversely proportional to the accuracy estimate of $y_r$. For example, if the percentage accuracy is the same for all $y_r$, then the absolute accuracy of $y_r$ is proportional to $y_r$ (assuming $y_r$ to be positive, as it usually is in such cases) and so $w_{\mathit{r}}=K/y_{\mathit{r}}$, for $\mathit{r}=1,2,\dots ,m$, for an arbitrary positive constant $K$. (This definition of weight is stressed because often weight is defined as the square of that used here.) The norms (2), (3) and (4) above are then replaced respectively by and Again it is the square of (6) which is used in practice rather than (6) itself.

Two different forms for representing a polynomial are used in different methods. One is the usual power-series form The other is the Chebyshev series form where $T_i\left(x\right)$ is the Chebyshev polynomial of the first kind of degree $i$ (see page 9 of Cox and Hayes (1973)), and where the range of $x$ has been normalized to run from $-1$ to $+1$. The use of either form leads theoretically to the same fitted polynomial, but in practice results may differ substantially because of the effects of rounding error. The Chebyshev form is to be preferred, since it leads to much better accuracy in general, both in the computation of the coefficients and in the subsequent evaluation of the fitted polynomial at specified points. This form also has other advantages: for example, since the later terms in (9) generally decrease much more rapidly from left to right than do those in (8), the situation is more often encountered where the last terms are negligible and it is obvious that the degree of the polynomial can be reduced (note that on the interval $-1\le x\le 1$ for all $i$, $T_i\left(x\right)$ attains the value unity but never exceeds it, so that the coefficient $a_i$ gives directly the maximum value of the term containing it). If the power-series form is used it is most advisable to work with the variable $x$ normalized to the range $-1$ to $+1$, carrying out the normalization before entering the relevant method. This will often substantially improve computational accuracy.

A cubic spline is represented in the form where $N_{\mathit{i}}\left(x\right)$, for $\mathit{i}=1,2,\dots ,p$, is a normalized cubic B-spline (see Hayes (1974)). This form, also, has advantages of computational speed and accuracy over alternative representations.

The type of bivariate polynomial currently considered in the chapter can be represented in either of the two forms and where $T_i\left(x\right)$ is the Chebyshev polynomial of the first kind of degree $i$ in the argument $x$ (see page 9 of Cox and Hayes (1973)), and correspondingly for $T_j\left(y\right)$. The prime on the two summation signs, following standard convention, indicates that the first term in each sum is halved, as shown for one variable in equation (9). The two forms (11) and (12) are mathematically equivalent, but again the Chebyshev form is to be preferred on numerical grounds, as discussed in [Representation of polynomials].

The bicubic spline is defined over a rectangle $R$ in the $\left(x,y\right)$ plane, the sides of $R$ being parallel to the $x$ and $y$ axes. $R$ is divided into rectangular panels, again by lines parallel to the axes. Over each panel the bicubic spline is a bicubic polynomial, that is it takes the form Each of these polynomials joins the polynomials in adjacent panels with continuity up to the second derivative. The constant $x$ values of the dividing lines parallel to the $y$-axis form the set of interior knots for the variable $x$, corresponding precisely to the set of interior knots of a cubic spline. Similarly, the constant $y$ values of dividing lines parallel to the $x$-axis form the set of interior knots for the variable $y$. Instead of representing the bicubic spline in terms of the above set of bicubic polynomials, however, it is represented, for the sake of computational speed and accuracy, in the form where $M_{\mathit{i}}\left(x\right)$, for $\mathit{i}=1,2,\dots ,p$, and $N_{\mathit{j}}\left(y\right)$, for $\mathit{j}=1,2,\dots ,q$, are normalized B-splines (see Hayes and Halliday (1974) for further details of bicubic splines and Hayes (1974) for normalized B-splines).