NAG CL Interface
e02dac (dim2_spline_panel)
1
Purpose
e02dac forms a minimal, weighted least squares bicubic spline surface fit with prescribed knots to a given set of data points.
2
Specification
| void |
e02dac (Integer m,
const double x[],
const double y[],
const double f[],
const double w[],
const Integer point[],
double dl[],
double eps,
double *sigma,
Integer *rank,
Nag_2dSpline *spline,
NagError *fail) |
|
The function may be called by the names: e02dac, nag_fit_dim2_spline_panel or nag_2d_spline_fit_panel.
3
Description
e02dac determines a bicubic spline fit
to the set of data points
with weights
, for
. The two sets of internal knots of the spline,
and
, associated with the variables
and
respectively, are prescribed by you. These knots can be thought of as dividing the data region of the
plane into panels (see
Figure 1 in
Section 5). A bicubic spline consists of a separate bicubic polynomial in each panel, the polynomials joining together with continuity up to the second derivative across the panel boundaries.
has the property that
, the sum of squares of its weighted residuals
, for
, where
is as small as possible for a bicubic spline with the given knot sets. The function produces this minimized value of
and the coefficients
in the B-spline representation of
– see
Section 9.
e02dec,
e02dfc and
e02dhc are available to compute values and derivatives of the fitted spline from the coefficients
.
The least squares criterion is not always sufficient to determine the bicubic spline uniquely: there may be a whole family of splines which have the same minimum sum of squares. In these cases, the function selects from this family the spline for which the sum of squares of the coefficients
is smallest: in other words, the minimal least squares solution. This choice, although arbitrary, reduces the risk of unwanted fluctuations in the spline fit. The method employed involves forming a system of
linear equations in the coefficients
and then computing its least squares solution, which will be the minimal least squares solution when appropriate. The basis of the method is described in
Hayes and Halliday (1974). The matrix of the equation is formed using a recurrence relation for B-splines which is numerically stable (see
Cox (1972) and
de Boor (1972) – the former contains the more elementary derivation but, unlike
de Boor (1972), does not cover the case of coincident knots). The least squares solution is also obtained in a stable manner by using orthogonal transformations, viz. a variant of Givens rotation (see
Gentleman (1973)). This requires only one row of the matrix to be stored at a time. Advantage is taken of the stepped-band structure which the matrix possesses when the data points are suitably ordered, there being at most sixteen nonzero elements in any row because of the definition of B-splines. First the matrix is reduced to upper triangular form and then the diagonal elements of this triangle are examined in turn. When an element is encountered whose square, divided by the mean squared weight, is less than a threshold
, it is replaced by zero and the rest of the elements in its row are reduced to zero by rotations with the remaining rows. The rank of the system is taken to be the number of nonzero diagonal elements in the final triangle, and the nonzero rows of this triangle are used to compute the minimal least squares solution. If all the diagonal elements are nonzero, the rank is equal to the number of coefficients
and the solution obtained is the ordinary least squares solution, which is unique in this case.
4
References
Cox M G (1972) The numerical evaluation of B-splines J. Inst. Math. Appl. 10 134–149
de Boor C (1972) On calculating with B-splines J. Approx. Theory 6 50–62
Gentleman W M (1973) Least squares computations by Givens transformations without square roots J. Inst. Math. Applic. 12 329–336
Hayes J G and Halliday J (1974) The least squares fitting of cubic spline surfaces to general data sets J. Inst. Math. Appl. 14 89–103
5
Arguments
-
1:
– Integer
Input
-
On entry: , the number of data points.
Constraint:
.
-
2:
– const double
Input
-
3:
– const double
Input
-
4:
– const double
Input
-
On entry: the coordinates of the data point
, for
. The order of the data points is immaterial, but see the array
point.
-
5:
– const double
Input
-
On entry: the weight
of the
th data point. It is important to note the definition of weight implied by the equation
(1) in
Section 3, since it is also common usage to define weight as the square of this weight. In this function, each
should be chosen inversely proportional to the (absolute) accuracy of the corresponding
, as expressed, for example, by the standard deviation or probable error of the
. When the
are all of the same accuracy, all the
may be set equal to
.
-
6:
– const Integer
Input
-
Note: the dimension,
dim, of the array
point
must be at least
.
On entry: indexing information usually provided by
e02zac which enables the data points to be accessed in the order which produces the advantageous matrix structure mentioned in
Section 3. This order is such that, if the
plane is thought of as being divided into rectangular panels by the two sets of knots, all data in a panel occur before data in succeeding panels, where the panels are numbered from bottom to top and then left to right with the usual arrangement of axes, as indicated in
Figure 1.
A data point lying exactly on one or more panel sides is considered to be in the highest numbered panel adjacent to the point.
e02zac should be called to obtain the array
point, unless it is provided by other means.
-
7:
– double
Output
-
Note: the dimension,
dim, of the array
dl
must be at least
.
On exit: gives the squares of the diagonal elements of the reduced triangular matrix, divided by the mean squared weight. It includes those elements, less than
, which are treated as zero (see
Section 3).
-
8:
– double
Input
-
On entry: a threshold
for determining the effective rank of the system of linear equations. The rank is determined as the number of elements of the array
dl which are nonzero. An element of
dl is regarded as zero if it is less than
.
Machine precision is a suitable value for
in most practical applications which have only
or
decimals accurate in data. If some coefficients of the fit prove to be very large compared with the data ordinates, this suggests that
should be increased so as to decrease the rank. The array
dl will give a guide to appropriate values of
to achieve this, as well as to the choice of
in other cases where some experimentation may be needed to determine a value which leads to a satisfactory fit.
-
9:
– double *
Output
-
On exit:
, the weighted sum of squares of residuals. This is not computed from the individual residuals but from the right-hand sides of the orthogonally-transformed linear equations. For further details see page 97 of
Hayes and Halliday (1974). The two methods of computation are theoretically equivalent, but the results may differ because of rounding error.
-
10:
– Integer *
Output
-
On exit: the rank of the system as determined by the value of the threshold
.
- The least squares solution is unique.
- The minimal least squares solution is computed.
-
11:
– Nag_2dSpline *
-
Pointer to structure of type Nag_2dSpline with the following members:
- nx – IntegerInput
-
On entry: must specify the total number of knots associated with the variables . It is such that is the number of interior knots.
Constraint:
.
- lamda – doubleInput/Output
-
On entry:
must contain the
th interior knot
associated with the variable
, for
. The knots must be in nondecreasing order and lie strictly within the range covered by the data values of
. A knot is a value of
at which the spline is allowed to be discontinuous in the third derivative with respect to
, though continuous up to the second derivative. This degree of continuity can be reduced, if you require, by the use of coincident knots, provided that no more than four knots are chosen to coincide at any point. Two, or three, coincident knots allow loss of continuity in, respectively, the second and first derivative with respect to
at the value of
at which they coincide. Four coincident knots split the spline surface into two independent parts. For choice of knots see
Section 9.
On exit: the interior knots
to
are unchanged, and the segments
and
contain additional (exterior) knots introduced by the function in order to define the full set of B-splines required. The four knots in the first segment are all set equal to the lowest data value of
and the other four additional knots are all set equal to the highest value: there is experimental evidence that coincident end-knots are best for numerical accuracy. The complete array must be left undisturbed if
e02dec or
e02dfc is to be used subsequently.
- ny – IntegerInput
-
On entry:
must specify the total number of knots associated with the variable
.
It is such that is the number of interior knots.
Constraint:
.
- mu – doubleInput/Output
-
On entry: must contain the th interior knot associated with the variable , .
On exit: the same remarks apply to
as to
above, with
y replacing
x, and
replacing
.
- c – doubleOutput
-
On exit: gives the coefficients of the fit.
is the coefficient
of
Sections 3 and
9, for
and
. These coefficients are used by
e02dec or
e02dfc to calculate values of the fitted function.
In normal usage, the call to
e02dac follows a call to
e01dac,
e02dcc or
e02ddc, in which case, members of the structure
spline will have been set up correctly for input to
e02dac.
-
12:
– NagError *
Input/Output
-
The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
- NE_ALLOC_FAIL
-
Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
- NE_BAD_PARAM
-
On entry, argument had an illegal value.
- NE_CONSTRAINT
-
Constraint: .
Constraint: .
- NE_INT
-
On entry, .
Constraint: .
- NE_INTERNAL_ERROR
-
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
- NE_KNOTS_COINCIDE
-
More than four knots coincide at a single point.
- NE_KNOTS_CONS
-
At least one set of knots is not in nondecreasing order.
- NE_NO_LICENCE
-
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
- NE_PANEL_ORDER
-
Array
point does not indicate the data points in panel order.
- NE_WEIGHT_ZERO
-
All the weights are zero, or rank determined as zero.
7
Accuracy
The computation of the B-splines and reduction of the observation matrix to triangular form are both numerically stable.
8
Parallelism and Performance
e02dac is not threaded in any implementation.
The time taken is approximately proportional to the number of data points, , and to .
The B-spline representation of the bicubic spline is
summed over
and over
. Here
and
denote normalized cubic B-splines, the former defined on the knots
and the latter on the knots
. For further details, see
Hayes and Halliday (1974) for bicubic splines and
de Boor (1972) for normalized B-splines.
The choice of the interior knots, which help to determine the spline's shape, must largely be a matter of trial and error. It is usually best to start with a small number of knots and, examining the fit at each stage, add a few knots at a time in places where the fit is particularly poor. In intervals of or where the surface represented by the data changes rapidly, in function value or derivatives, more knots will be needed than elsewhere. In some cases guidance can be obtained by analogy with the case of coincident knots: for example, just as three coincident knots can produce a discontinuity in slope, three close knots can produce rapid change in slope. Of course, such rapid changes in behaviour must be adequately represented by the data points, as indeed must the behaviour of the surface generally, if a satisfactory fit is to be achieved. When there is no rapid change in behaviour, equally-spaced knots will often suffice.
In all cases the fit should be examined graphically before it is accepted as satisfactory.
The fit obtained is not defined outside the rectangle
The reason for taking the extreme data values of
and
for these four knots is that, as is usual in data fitting, the fit cannot be expected to give satisfactory values outside the data region. If, nevertheless, you require values over a larger rectangle, this can be achieved by augmenting the data with two artificial data points
and
with zero weight, where
,
defines the enlarged rectangle. In the case when the data are adequate to make the least squares solution unique (
), this enlargement will not affect the fit over the original rectangle, except for possibly enlarged rounding errors, and will simply continue the bicubic polynomials in the panels bordering the rectangle out to the new boundaries: in other cases the fit will be affected. Even using the original rectangle there may be regions within it, particularly at its corners, which lie outside the data region and where, therefore, the fit will be unreliable. For example, if there is no data point in panel
of
Figure 1 in
Section 5, the least squares criterion leaves the spline indeterminate in this panel: the minimal spline determined by the function in this case passes through the value zero at the point
.
10
Example
This example reads a value for
, and a set of data points, weights and knot positions. If there are more
knots than
knots, it interchanges the
and
axes. It calls
e02zac to sort the data points into panel order,
e02dac to fit a bicubic spline to them, and
e02dec to evaluate the spline at the data points.
Finally it prints:
-
–the weighted sum of squares of residuals computed from the linear equations;
-
–the rank determined by e02dac;
-
–data points, fitted values and residuals in panel order;
-
–the weighted sum of squares of the residuals; and
-
–the coefficients of the spline fit.
The program is written to handle any number of datasets.
Note: the data supplied in this example is
not typical of a realistic problem: the number of data points would normally be much larger (in which case the array dimensions would have to be increased); and the value of
would normally be much smaller on most machines (see
Section 5; the relatively large value of
has been chosen in order to illustrate a minimal least squares solution when
; in this example
).
10.1
Program Text
10.2
Program Data
10.3
Program Results
