BSLSQ

Computes the least-squares spline approximation, and return the B-spline coefficients.

Required Arguments

XDATA — Array of length NDATA containing the data point abscissas. (Input)

FDATA — Array of length NDATA containing the data point ordinates. (Input)

KORDER — Order of the spline. (Input)

KORDER must be less than or equal to NDATA.

KORDER must be less than or equal to NDATA.

XKNOT — Array of length NCOEF + KORDER containing the knot sequence. (Input)

XKNOT must be nondecreasing.

XKNOT must be nondecreasing.

NCOEF — Number of B-spline coefficients. (Input)

NCOEF cannot be greater than NDATA.

NCOEF cannot be greater than NDATA.

BSCOEF — Array of length NCOEF containing the B-spline coefficients. (Output)

Optional Arguments

NDATA — Number of data points. (Input)

Default: NDATA = size(XDATA, 1)

Default: NDATA = size(XDATA, 1)

WEIGHT — Array of length NDATA containing the weights. (Input)

Default: WEIGHT = 1.0.

Default: WEIGHT = 1.0.

FORTRAN 90 Interface

Generic: CALL BSLSQ (XDATA, FDATA, KORDER, XKNOT, NCOEF, BSCOEF [, …])

Specific: The specific interface names are S_BSLSQ and D_BSLSQ.

FORTRAN 77 Interface

Single: CALL BSLSQ (NDATA, XDATA, FDATA, WEIGHT, KORDER, XKNOT, NCOEF, BSCOEF)

Double: The double precision name is DBSLSQ.

Description

The routine BSLSQ is based on the routine L2APPR by de Boor (1978, page 255). The IMSL routine BSLSQ computes a weighted discrete L2 approximation from a spline subspace to a given data set (xi, fi) for i = 1, …, N (where N = NDATA). In other words, it finds B-spline coefficients,

a = BSCOEF, such that

a = BSCOEF, such that

is a minimum, where m = NCOEF and Bj denotes the j-th B-spline for the given order, KORDER, and knot sequence, XKNOT. This linear least squares problem is solved by computing and solving the normal equations. While the normal equations can sometimes cause numerical difficulties, their use here should not cause a problem because the B-spline basis generally leads to well-conditioned banded matrices.

The choice of weights depends on the problem. In some cases, there is a natural choice for the weights based on the relative importance of the data points. To approximate a continuous function (if the location of the data points can be chosen), then the use of Gauss quadrature weights and points is reasonable. This follows because BSLSQ is minimizing an approximation to the integral

The Gauss quadrature weights and points can be obtained using the IMSL routine GQRUL (see Chapter 4, “Integration and Differentiation”).

Comments

1. Workspace may be explicitly provided, if desired, by use of B2LSQ/DB2LSQ. The reference is:

CALL B2LSQ (NDATA, XDATA, FDATA, WEIGHT, KORDER, XKNOT, NCOEF, BSCOEF, WK1, WK2, WK3, WK4, IWK)

The additional arguments are as follows:

WK1 — Work array of length (3 + NCOEF) * KORDER.

WK2 — Work array of length NDATA.

WK3 — Work array of length NDATA.

WK4 — Work array of length NDATA.

IWK — Work array of length NDATA.

2. Informational errors

Type | Code | Description |
---|---|---|

4 | 5 | Multiplicity of the knots cannot exceed the order of the spline. |

4 | 6 | The knots must be nondecreasing. |

4 | 7 | All weights must be greater than zero. |

4 | 8 | The smallest element of the data point array must be greater than or equal to the KORDth knot. |

4 | 9 | The largest element of the data point array must be less than or equal to the (NCOEF + 1)st knot. |

3. The B-spline representation can be evaluated using BSVAL, and its derivative can be evaluated using BSDER.

Example

In this example, we try to recover a quadratic polynomial using a quadratic spline with one interior knot from two different data sets. The first data set is generated by evaluating the quadratic at 50 equally spaced points in the interval (0, 1) and then adding uniformly distributed noise to the data. The second data set includes the first data set, and, additionally, the values at 0 and at 1 with no noise added. Since the first and last data points are uncontaminated by noise, we have chosen weights equal to 105 for these two points in this second problem. The quadratic, the first approximation, and the second approximation are then evaluated at 11 equally spaced points. This example illustrates the use of the weights to enforce interpolation at certain of the data points.

USE IMSL_LIBRARIES

IMPLICIT NONE

INTEGER KORDER, NCOEF

PARAMETER (KORDER=3, NCOEF=4)

!

INTEGER I, NDATA, NOUT

REAL ABS, BSCOF1(NCOEF), BSCOF2(NCOEF), F,&

FDATA1(50), FDATA2(52), FLOAT, RNOISE, S1,&

S2, WEIGHT(52), X, XDATA1(50), XDATA2(52),&

XKNOT(KORDER+NCOEF), XT, YT

INTRINSIC ABS, FLOAT

!

DATA WEIGHT/52*1.0/

! Define function

F(X) = 8.0*X*(1.0-X)

! Set random number seed

CALL RNSET (12345679)

NDATA = 50

! Set up interior knots

DO 10 I=1, NCOEF - KORDER + 2

XKNOT(I+KORDER-1) = FLOAT(I-1)/FLOAT(NCOEF-KORDER+1)

10 CONTINUE

! Stack knots

DO 20 I=1, KORDER - 1

XKNOT(I) = XKNOT(KORDER)

XKNOT(I+NCOEF+1) = XKNOT(NCOEF+1)

20 CONTINUE

! Set up data points excluding

! the endpoints 0 and 1.

! The function values have noise

! introduced.

DO 30 I=1, NDATA

XDATA1(I) = FLOAT(I)/51.0

RNOISE = RNUNF()

RNOISE = RNOISE - 0.5

FDATA1(I) = F(XDATA1(I)) + RNOISE

30 CONTINUE

! Compute least squares B-spline

! representation.

CALL BSLSQ (XDATA1, FDATA1, KORDER, XKNOT, NCOEF, BSCOF1)

! Now use same XDATA values but with

! the endpoints included. These

! points will have large weights.

NDATA = 52

CALL SCOPY (50, XDATA1, 1, XDATA2(2:), 1)

CALL SCOPY (50, FDATA1, 1, FDATA2(2:), 1)

!

WEIGHT(1) = 1.0E5

XDATA2(1) = 0.0

FDATA2(1) = F(XDATA2(1))

WEIGHT(NDATA) = 1.0E5

XDATA2(NDATA) = 1.0

FDATA2(NDATA) = F(XDATA2(NDATA))

! Compute least squares B-spline

! representation.

CALL BSLSQ (XDATA2, FDATA2, KORDER, XKNOT, NCOEF, BSCOF2, &

WEIGHT=WEIGHT)

! Get output unit number

CALL UMACH (2, NOUT)

! Write heading

WRITE (NOUT,99998)

! Print the two interpolants

! at 11 points.

DO 40 I=1, 11

XT = FLOAT(I-1)/10.0

YT = F(XT)

! Evaluate splines

S1 = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOF1)

S2 = BSVAL(XT,KORDER,XKNOT,NCOEF,BSCOF2)

WRITE (NOUT,99999) XT, YT, S1, S2, (S1-YT), (S2-YT)

40 CONTINUE

!

99998 FORMAT (7X, 'X', 9X, 'F(X)', 6X, 'S1(X)', 5X, 'S2(X)', 7X,&

'F(X)-S1(X)', 7X, 'F(X)-S2(X)')

99999 FORMAT (' ', 4F10.4, 4X, F10.4, 7X, F10.4)

END

Output

X F(X) S1(X) S2(X) F(X)-S1(X) F(X)-S2(X)

0.0000 0.0000 0.0515 0.0000 0.0515 0.0000

0.1000 0.7200 0.7594 0.7490 0.0394 0.0290

0.2000 1.2800 1.3142 1.3277 0.0342 0.0477

0.3000 1.6800 1.7158 1.7362 0.0358 0.0562

0.4000 1.9200 1.9641 1.9744 0.0441 0.0544

0.5000 2.0000 2.0593 2.0423 0.0593 0.0423

0.6000 1.9200 1.9842 1.9468 0.0642 0.0268

0.7000 1.6800 1.7220 1.6948 0.0420 0.0148

0.8000 1.2800 1.2726 1.2863 -0.0074 0.0063

0.9000 0.7200 0.6360 0.7214 -0.0840 0.0014

1.0000 0.0000 -0.1878 0.0000 -0.1878 0.0000