hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_fit_2dspline_panel (e02da)

Purpose

nag_fit_2dspline_panel (e02da) forms a minimal, weighted least squares bicubic spline surface fit with prescribed knots to a given set of data points.

Syntax

[lamda, mu, dl, c, sigma, rank, ifail] = e02da(x, y, f, w, lamda, mu, point, eps, 'm', m, 'px', px, 'py', py, 'npoint', npoint)
[lamda, mu, dl, c, sigma, rank, ifail] = nag_fit_2dspline_panel(x, y, f, w, lamda, mu, point, eps, 'm', m, 'px', px, 'py', py, 'npoint', npoint)

Description

nag_fit_2dspline_panel (e02da) determines a bicubic spline fit s(x,y)s(x,y) to the set of data points (xr,yr,fr)(xr,yr,fr) with weights wrwr, for r = 1,2,,mr=1,2,,m. The two sets of internal knots of the spline, {λ}{λ} and {μ}{μ}, associated with the variables xx and yy respectively, are prescribed by you. These knots can be thought of as dividing the data region of the (x,y)(x,y) plane into panels (see Figure 1 in Section [Parameters]). 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.
s(x,y)s(x,y) has the property that ΣΣ, the sum of squares of its weighted residuals ρrρr, for r = 1,2,,mr=1,2,,m, where
ρr = wr (s(xr,yr)fr)
ρr = wr ( s ( x r , y r ) - f r )
(1)
is as small as possible for a bicubic spline with the given knot sets. The function produces this minimized value of ΣΣ and the coefficients cijcij in the B-spline representation of s(x,y)s(x,y) – see Section [Further Comments]. nag_fit_2dspline_evalv (e02de), nag_fit_2dspline_evalm (e02df) and nag_fit_2dspline_derivm (e02dh) are available to compute values and derivatives of the fitted spline from the coefficients cijcij.
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 cijcij 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 mm linear equations in the coefficients cijcij 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 cijcij and the solution obtained is the ordinary least squares solution, which is unique in this case.

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

Parameters

Compulsory Input Parameters

1:     x(m) – double array
2:     y(m) – double array
3:     f(m) – double array
m, the dimension of the array, must satisfy the constraint m > 1m>1.
The coordinates of the data point (xr,yr,fr)(xr,yr,fr), for r = 1,2,,mr=1,2,,m. The order of the data points is immaterial, but see the array point.
4:     w(m) – double array
m, the dimension of the array, must satisfy the constraint m > 1m>1.
The weight wrwr of the rrth data point. It is important to note the definition of weight implied by the equation (1) in Section [Description], since it is also common usage to define weight as the square of this weight. In this function, each wrwr should be chosen inversely proportional to the (absolute) accuracy of the corresponding frfr, as expressed, for example, by the standard deviation or probable error of the frfr. When the frfr are all of the same accuracy, all the wrwr may be set equal to 1.01.0.
5:     lamda(px) – double array
px, the dimension of the array, must satisfy the constraint px8px8 and py8py8.
(They are such that px8px-8 and py8py-8 are the corresponding numbers of interior knots.) The running time and storage required by the function are both minimized if the axes are labelled so that py is the smaller of px and py.
lamda(i + 4)lamdai+4 must contain the iith interior knot λi + 4λi+4 associated with the variable xx, for i = 1,2,,px8i=1,2,,px-8. The knots must be in nondecreasing order and lie strictly within the range covered by the data values of xx. A knot is a value of xx at which the spline is allowed to be discontinuous in the third derivative with respect to xx, 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 xx at the value of xx at which they coincide. Four coincident knots split the spline surface into two independent parts. For choice of knots see Section [Further Comments].
6:     mu(py) – double array
py, the dimension of the array, must satisfy the constraint px8px8 and py8py8.
(They are such that px8px-8 and py8py-8 are the corresponding numbers of interior knots.) The running time and storage required by the function are both minimized if the axes are labelled so that py is the smaller of px and py.
mu(i + 4)mui+4 must contain the iith interior knot μi + 4μi+4 associated with the variable yy, for i = 1,2,,py8i=1,2,,py-8.
7:     point(npoint) – int64int32nag_int array
npoint, the dimension of the array, must satisfy the constraint npointm + (px7) × (py7)npointm+(px-7)×(py-7).
Indexing information usually provided by nag_fit_2dspline_sort (e02za) which enables the data points to be accessed in the order which produces the advantageous matrix structure mentioned in Section [Description]. This order is such that, if the (x,y)(x,y) 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.
Figure 1
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. nag_fit_2dspline_sort (e02za) should be called to obtain the array point, unless it is provided by other means.
8:     eps – double scalar
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 22 or 33 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.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The dimension of the arrays x, y, f, w. (An error is raised if these dimensions are not equal.)
mm, the number of data points.
Constraint: m > 1m>1.
2:     px – int64int32nag_int scalar
3:     py – int64int32nag_int scalar
Default: For px, the dimension of the array lamda. For py, the dimension of the array mu.
The total number of knots λλ and μμ associated with the variables xx and yy, respectively.
Constraint: px8px8 and py8py8.
(They are such that px8px-8 and py8py-8 are the corresponding numbers of interior knots.) The running time and storage required by the function are both minimized if the axes are labelled so that py is the smaller of px and py.
4:     npoint – int64int32nag_int scalar
Default: The dimension of the array point.
The dimension of the array point as declared in the (sub)program from which nag_fit_2dspline_panel (e02da) is called.
Constraint: npointm + (px7) × (py7)npointm+(px-7)×(py-7).

Input Parameters Omitted from the MATLAB Interface

nc ws nws

Output Parameters

1:     lamda(px) – double array
The interior knots lamda(5)lamda5 to lamda(px4)lamdapx-4 are unchanged, and the segments lamda(1 : 4)lamda1:4 and lamda(px3 : px)lamdapx-3:px 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 xx 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 nag_fit_2dspline_evalv (e02de) or nag_fit_2dspline_evalm (e02df) is to be used subsequently.
2:     mu(py) – double array
The same remarks apply to mu as to lamda above, with y replacing x, and yy replacing xx.
3:     dl(nc) – double array
nc = (px4) × (py4)nc=(px-4)×(py-4).
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 [Description]).
4:     c(nc) – double array
nc = (px4) × (py4)nc=(px-4)×(py-4).
Gives the coefficients of the fit. c((py4) × (i1) + j)c((py-4)×(i-1)+j) is the coefficient cijcij of Sections [Description] and [Further Comments], for i = 1,2,,px4i=1,2,,px-4 and j = 1,2,,py4j=1,2,,py-4. These coefficients are used by nag_fit_2dspline_evalv (e02de) or nag_fit_2dspline_evalm (e02df) to calculate values of the fitted function.
5:     sigma – double scalar
ΣΣ, 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.
6:     rank – int64int32nag_int scalar
The rank of the system as determined by the value of the threshold εε.
rank = ncrank=nc
The least squares solution is unique.
rankncranknc
The minimal least squares solution is computed.
7:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
At least one set of knots is not in nondecreasing order, or an interior knot is outside the range of the data values.
  ifail = 2ifail=2
More than four knots coincide at a single point, possibly because all data points have the same value of xx (or yy) or because an interior knot coincides with an extreme data value.
  ifail = 3ifail=3
Array point does not indicate the data points in panel order. Call nag_fit_2dspline_sort (e02za) to obtain a correct array.
  ifail = 4ifail=4
On entry,m1m1,
orpx < 8px<8,
orpy < 8py<8,
ornc(px4) × (py4)nc(px-4)×(py-4),
ornws is too small,
ornpoint is too small.
  ifail = 5ifail=5
All the weights wrwr are zero or rank determined as zero.

Accuracy

The computation of the B-splines and reduction of the observation matrix to triangular form are both numerically stable.

Further Comments

The time taken is approximately proportional to the number of data points, mm, and to (3 × (py4) + 4)2(3×(py-4)+4)2.
The B-spline representation of the bicubic spline is
s(x,y) =  cijMi(x)Nj(y)
i,j
s(x,y) = i,j cij Mi(x) Nj(y)
summed over i = 1,2,,px4i=1,2,,px-4 and over j = 1,2,,py4j=1,2,,py-4. Here Mi(x)Mi(x) and Nj(y)Nj(y) denote normalized cubic B-splines, the former defined on the knots λi,λi + 1,,λi + 4λi,λi+1,,λi+4 and the latter on the knots μj,μj + 1,,μj + 4μj,μj+1,,μj+4. 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 xx or yy 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
λ4 x λpx3 , μ4 y μpy3 .
λ4 x λpx-3 , μ4 y μpy-3 .
The reason for taking the extreme data values of xx and yy 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 (a,c,0)(a,c,0) and (b,d,0)(b,d,0) with zero weight, where axbaxb, cydcyd defines the enlarged rectangle. In the case when the data are adequate to make the least squares solution unique (rank = ncrank=nc), 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 11 of Figure 1 in Section [Parameters], 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 (λ4,μ4)(λ4,μ4).

Example

function nag_fit_2dspline_panel_example
x = [0.6;
     -0.95;
     0.87;
     0.84;
     0.17;
     -0.87;
     1;
     0.1;
     0.24;
     -0.77;
     0.32;
     1;
     -0.63;
     -0.66;
     0.93;
     0.15;
     0.99;
     -0.54;
     0.44;
     -0.72;
     0.63;
     -0.4;
     0.2;
     0.43;
     0.28;
     -0.24;
     0.86;
     -0.41;
     -0.05;
     -1];
y = [-0.52;
     -0.61;
     0.93;
     0.09;
     0.88;
     -0.7;
     1;
     1;
     0.3;
     -0.77;
     -0.23;
     -1;
     -0.26;
     -0.83;
     0.22;
     0.89;
     -0.8;
     -0.88;
     0.68;
     -0.14;
     0.67;
     -0.9;
     -0.84;
     0.84;
     0.15;
     -0.91;
     -0.35;
     -0.16;
     -0.35;
     -1];
f = [0.93;
     -1.79;
     0.36;
     0.52;
     0.49;
     -1.76;
     0.33;
     0.48;
     0.65;
     -1.82;
     0.92;
     1;
     8.88;
     -2.01;
     0.47;
     0.49;
     0.84;
     -2.42;
     0.47;
     7.15;
     0.44;
     -3.34;
     2.78;
     0.44;
     0.7;
     -6.52;
     0.66;
     2.32;
     1.66;
     -1];
w = [10;
     10;
     10;
     10;
     10;
     10;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1;
     1];
lamda = zeros(10,1);
lamda(5) = -0.5;
mu = zeros(8, 1);
point = [int64(3);6;4;5;7;10;8;9;11;13;12;15;14;18;16;17;19;20;21;30;23;26;24;25;27;28; ...
       0;29;0;0;2;22;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
epsilon = 1e-06;
[lamdaOut, muOut, dl, c, sigma, rank, ifail] = ...
    nag_fit_2dspline_panel(x, y, f, w, lamda, mu, point, epsilon)
 

lamdaOut =

   -1.0000
   -1.0000
   -1.0000
   -1.0000
   -0.5000
         0
    1.0000
    1.0000
    1.0000
    1.0000


muOut =

    -1
    -1
    -1
    -1
     1
     1
     1
     1


dl =

    1.0417
    0.0286
    0.0002
    0.0000
    0.0318
    0.0030
    0.0000
    0.0000
    0.0173
    0.0035
    0.0043
    0.0072
    0.0395
    0.0045
    0.0011
    0.0039
    0.0103
    0.0144
    0.0009
    0.0223
    0.0724
    0.0041
    0.0019
    0.0491


c =

   -1.0228
  115.4668
 -433.5558
  -68.1973
   24.8426
 -140.1485
  258.5042
   15.6756
  -29.4878
  132.2933
 -173.5103
   20.0983
    9.9575
  -51.6200
   67.6666
   -5.8765
   10.0577
    4.7543
  -15.3533
   -0.3260
    1.0835
   -2.7932
    7.7708
    0.6315


sigma =

   14.6671


rank =

                   22


ifail =

                    0


function e02da_example
x = [0.6;
     -0.95;
     0.87;
     0.84;
     0.17;
     -0.87;
     1;
     0.1;
     0.24;
     -0.77;
     0.32;
     1;
     -0.63;
     -0.66;
     0.93;
     0.15;
     0.99;
     -0.54;
     0.44;
     -0.72;
     0.63;
     -0.4;
     0.2;
     0.43;
     0.28;
     -0.24;
     0.86;
     -0.41;
     -0.05;
     -1];
y = [-0.52;
     -0.61;
     0.93;
     0.09;
     0.88;
     -0.7;
     1;
     1;
     0.3;
     -0.77;
     -0.23;
     -1;
     -0.26;
     -0.83;
     0.22;
     0.89;
     -0.8;
     -0.88;
     0.68;
     -0.14;
     0.67;
     -0.9;
     -0.84;
     0.84;
     0.15;
     -0.91;
     -0.35;
     -0.16;
     -0.35;
     -1];
f = [0.93;
     -1.79;
     0.36;
     0.52;
     0.49;
     -1.76;
     0.33;
     0.48;
     0.65;
     -1.82;
     0.92;
     1;
     8.88;
     -2.01;
     0.47;
     0.49;
     0.84;
     -2.42;
     0.47;
     7.15;
     0.44;
     -3.34;
     2.78;
     0.44;
     0.7;
     -6.52;
     0.66;
     2.32;
     1.66;
     -1];
w = [10;10;10;10;10;10;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
lamda = zeros(10,1);
lamda(5) = -0.5;
mu = zeros(8, 1);
point = [int64(3);6;4;5;7;10;8;9;11;13;12;15;14;18;16;17;19;20;21;30;23;26; ...
         24;25;27;28;0;29;0;0;2;22;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
epsilon = 1e-06;
[lamdaOut, muOut, dl, c, sigma, rank, ifail] = ...
    e02da(x, y, f, w, lamda, mu, point, epsilon)
 

lamdaOut =

   -1.0000
   -1.0000
   -1.0000
   -1.0000
   -0.5000
         0
    1.0000
    1.0000
    1.0000
    1.0000


muOut =

    -1
    -1
    -1
    -1
     1
     1
     1
     1


dl =

    1.0417
    0.0286
    0.0002
    0.0000
    0.0318
    0.0030
    0.0000
    0.0000
    0.0173
    0.0035
    0.0043
    0.0072
    0.0395
    0.0045
    0.0011
    0.0039
    0.0103
    0.0144
    0.0009
    0.0223
    0.0724
    0.0041
    0.0019
    0.0491


c =

   -1.0228
  115.4668
 -433.5558
  -68.1973
   24.8426
 -140.1485
  258.5042
   15.6756
  -29.4878
  132.2933
 -173.5103
   20.0983
    9.9575
  -51.6200
   67.6666
   -5.8765
   10.0577
    4.7543
  -15.3533
   -0.3260
    1.0835
   -2.7932
    7.7708
    0.6315


sigma =

   14.6671


rank =

                   22


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013