.p>.CMCH3.DOC!SPLINE_LSQ_CONSTRAINED;spline_lsq_constrained
Computes a least-squares constrained spline approximation.
Synopsis
#include <imsl.h>
Imsl_f_spline *imsl_f_spline_lsq_constrained (int ndata, float xdata[], float fdata[], int spline_space_dim, int num_con_pts, f_constraint_struct constraints[], ¼, 0)
The type Imsl_d_spline function is imsl_d_spline_lsq_constrained.
Required Arguments
int ndata
(Input)
Number of data points.
float xdata[]
(Input)
Array with ndata components
containing the abscissas of the least-squares problem.
float fdata[]
(Input)
Array with ndata components
containing the ordinates of the least-squares problem.
int
spline_space_dim (Input)
The linear dimension of the
spline subspace. It should be smaller than ndata and greater than
or equal to order (whose default
value is 4).
int
num_con_pts (Input)
The number of points in the vector
constraints.
f_constraint_struct
constraints[] (Input)
A structure containing the abscissas
at which the fit is to be constrained, the derivative of the spline that is to
be constrained, the type of constraints, and any lower or upper limits. A
description of the structure fields follows:
|
Field |
Description |
|
xval |
point at which fit is constrained |
|
der |
derivative value of the spline to be constrained |
|
type |
types of the general constraints |
|
bl |
lower limit of the general constraints |
|
bu |
upper limit of the general constraints |
Notes: If you want to constrain the integral of the spline over the closed interval (c, d), then set constraints[i].der = constraints [i+1].der = −1 and constraints[i].xval = c and constraints[i+1].xval = d. For consistency, insist that constraints[i].type = constraints[i+1].type ³ 0 and c ≤ d. Note that every der must be at least −1.
|
constraints [i].type |
i-th constraint |
|
1 |
|
|
2 |
|
|
3 |
|
|
4 |
|
|
5 |
|
|
6 |
|
|
7 |
|
|
8 |
|
|
20 |
periodic end conditions |
|
99 |
disregard this constraint |
In order to have two point constraints, must have
constraints[i].type = constraints[i+1].type
|
constraints [i]. type |
i-th constraint |
|
9 |
|
|
10 |
|
|
11 |
|
|
12 |
|
Return Value
A pointer to the structure that represents the spline fit. If a fit cannot be computed, then NULL is returned. To release this space, use free.
Synopsis with Optional Arguments
#include <imsl.h>
Imsl_f_spline
*imsl_f_spline_lsq_constrained (int
ndata,
float xdata[], float
fdata[],
int
spline_space_dim, int num_con_pts, f_constraint_struct
constraints[],
IMSL_NHARD, int
nhard,
IMSL_WEIGHTS, float
weights[],
IMSL_ORDER, int
order,
IMSL_KNOTS, float
knots[],
0)
Optional Arguments
IMSL_NHARD, int nhard
(Output)
The argument nhard is the number of
entries of constraints involved in the “hard” constraints. Note that
0 ≤ nhard ≤ num_con_pts. The
default, nhard = 0,
always results in a fit, while setting nhard = num_con_pts forces all
constraints to be met. The “hard” constraints must be met, or else the function
signals failure. The “soft” constraints need not be satisfied, but there will be
an attempt to satisfy the “soft” constraints. The constraints must be listed in
terms of priority with the most important constraints first. Thus, all of the
“hard” constraints must precede the “soft” constraints. If infeasibility is
detected among the “soft” constraints, we satisfy, in order, as many of the
“soft” constraints as possible.
Default: nhard = 0
IMSL_WEIGHTS, float weights[]
(Input)
This option requires the user to provide the weights.
Default: all
weights equal one
IMSL_ORDER, int order
(Input)
The order of the spline subspace for which the knots are desired.
This option is used to communicate the order of the spline subspace.
Default:
order = 4(i.e.,
cubic splines)
IMSL_KNOTS, float knots[]
(Input)
This option requires the user to provide the knots. The user must
provide a knot sequence of length spline_space_dimension + order.
Default: an
appropriate knot sequence is selected. See below for more details.
Description
The function imsl_f_spline_lsq_constrained produces a constrained, weighted least-squares fit to data from a spline subspace. Constraints involving one point, two points, or integrals over an interval are allowed. The types of constraints supported by the functions are of four types:
|
Ep[f] |
|
|
Or |
|
|
Or |
|
|
Or |
= periodic end conditions |
An interval, Ip (which may be a point, a finite interval, or a semi-infinite interval), is associated with each of these constraints.
The
input for this function consists of several items; first, the data set
(xi, fi) for i = 1, ¼, N (where N = NDATA), that is the
data which is to be fit. Second, we have the weights to be used in the
least-squares fit (w = WEIGHT, defaulting to
1). The vector constraints contains
the abscissas of the points involved in specifying the constraints, as well as
information relating the type of constraints and the constraint interval.
Let nf denote the number of feasible constraints as described above. Then, the function solved the problem
subject to
This linearly constrained least-squares problem is treated as a quadratic program and is solved by invoking the function imsl_f_quadratic_prog (See Chapter 8, “Optimization”)
The choice of weights depends on the data uncertainty in the problem. In some cases, there is a natural choice for the weights based on the estimates of errors in the data points.
Determining feasibility of linear constraints is a numerically sensitive task. If you encounter difficulties, a quick fix would be to widen the constraint intervals Ip.
Examples
Example 1
This is a simple application of imsl_f_lsq_constrained. Data is generated from the function
and contaminated with random noise and fit with cubic splines. The function is increasing, so least-squares fit should also be increasing. This is not the case for the unconstrained least-squares fit generated by imsl_f_spline_least_squares. Then, the derivative is forced to be greater than 0 at num_con_pts = 15 equally spaced points and imsl_f_lsq_constrained is called. The resulting curve is monotone. The error is printed for the two fits averaged over 100 equally spaced points.
#include <imsl.h>
#include
<math.h>
#define MXKORD 4
#define MXNCOF
20
#define MXNDAT 51
#define MXNXVL
15
main()
{
f_constraint_struct
constraint[MXNXVL];
int i, korder, ncoef,
ndata, nxval;
float *noise, errlsq, errnft, grdsiz,
x;
float fdata[MXNDAT],
xdata[MXNDAT];
Imsl_f_spline *sp,
*spls;
#define F1(x) (float)(.5*(x) + sin( .5*(x)
))
korder = 4;
ndata =
15;
nxval = 15;
ncoef =
8;
/*
* Compute original xdata
and fdata with random noise.
*/
imsl_random_seed_set (234579);
noise = imsl_f_random_uniform (ndata, 0);
grdsiz =
10.0;
for (i = 0; i < ndata; i++)
{
xdata[i] = grdsiz * ((float) (i)
/ (float) (ndata - 1));
fdata[i] =
F1 (xdata[i]) + (noise[i] - .5);
}
/* Compute least-squares fit.
*/
spls = imsl_f_spline_least_squares (ndata, xdata,
fdata, ncoef, 0);
/*
*
Construct the constraints.
*/
for (i = 0; i < nxval; i++) {
constraint[i].xval = grdsiz * (float)(i) / (float)(nxval -
1);
constraint[i].type =
3;
constraint[i].der =
1;
constraint[i].bl =
0.0;
}
/* Compute constrained
least-squares fit. */
sp = imsl_f_spline_lsq_constrained
(ndata, xdata, fdata,
ncoef,
nxval, constraint, 0);
/*
*
Compute the average error of 100 points in the
interval.
*/
errlsq =
0.0;
errnft = 0.0;
for (i = 0; i
< 100; i++) {
x = grdsiz *
(float) (i) / 99.0;
errnft += fabs
(F1 (x) -
imsl_f_spline_value(x,sp,0));
errlsq += fabs (F1 (x) - imsl_f_spline_value(x,spls,0));
}
/* Print results */
printf ("
Average error with spline_least_squares fit:
%8.5f\n",
errlsq / 100.0);
printf (" Average error with
spline_lsq_constrained fit:
%8.5f\n",
errnft / 100.0);
}
Output
Average error with spline_least_squares
fit: 0.20250
Average error with
spline_lsq_constrained fit: 0.14334
Example 2
Now, try to recover the function
from noisy data. First, try the unconstrained least-squares fit using imsl_f_spline_least_squares. Finding that fit somewhat unsatisfactory, several constraints are applied using imsl_f_spline_lsq_constrained. First, notice that the unconstrained fit oscillates through the true function at both ends of the interval. This is common for flat data. To remove this oscillation, the cubic spline is constrained to have zero second derivative at the first and last four knots. This forces the cubic spline to reduce to a linear polynomial on the first and last three knot intervals. In addition, the fit is constrained (called s) as follows:
s(-7) ³ 0
s(-7) = s(7)
Notice that the last constraint was generated using the periodic option (requiring only the zero-th derivative to be periodic). The error is printed for the two fits averaged over 100 equally spaced points.
#include <imsl.h>
#include
<math.h>
#define KORDER 4
#define NDATA
51
#define NXVAL 12
#define NCOEF
13
main()
{
f_constraint_struct
constraint[NXVAL];
int
i;
float *noise, errlsq, errnft, grdsiz,
x;
float fdata[NDATA], xdata[NDATA],
xknot[NDATA+KORDER];
Imsl_f_spline *sp,
*spls;
#define F1(x)
(float)(1.0/(1.0+x*x*x*x))
/* Compute original
xdata and fdata with random noise */
imsl_random_seed_set (234579);
noise =
imsl_f_random_uniform (NDATA, 0);
grdsiz =
14.0;
for (i = 0; i < NDATA; i++)
{
xdata[i] = grdsiz *
((float)(i)/(float)(NDATA - 1))
- grdsiz/2.0;
fdata[i] = F1
(xdata[i]) + 0.125*(noise[i] - .5);
}
/* Generate knots. */
for (i = 0; i
< NCOEF-KORDER+2; i++) {
xknot[i+KORDER-1] = grdsiz *
((float)(i)/
(float)(NCOEF-KORDER+1)) - grdsiz/2.0;
}
for (i = 0; i < KORDER - 1; i++)
{
xknot[i] =
xknot[KORDER-1];
xknot[i+NCOEF+1]
= xknot[NCOEF];
}
/* Compute
spline_least_squares fit */
spls =
imsl_f_spline_least_squares (NDATA, xdata, fdata, NCOEF,
IMSL_KNOTS, xknot, 0);
/* Construct the
constraints for CONFT */
for (i = 0; i < 4; i++)
{
constraint[i].xval =
xknot[KORDER+i-1];
constraint[i+4].xval =
xknot[NCOEF-3+i];
constraint[i].itype = 1;
constraint[i+4].itype = 1;
constraint[i].ider = 2;
constraint[i+4].ider = 2;
constraint[i].bl = 0.0;
constraint[i+4].bl = 0.0;
}
constraint[8].xval = -7.0;
constraint[8].itype =
3;
constraint[8].ider = 0;
constraint[8].bl = 0.0;
constraint[9].xval =
-7.0;
constraint[9].itype = 6;
constraint[9].bu = 2.3;
constraint[10].xval =
7.0;
constraint[10].itype = 6;
constraint[10].bu = 2.3;
constraint[11].xval =
-7.0;
constraint[11].itype = 20;
constraint[11].ider = 0;
sp =
imsl_f_spline_lsq_constrained (NDATA, xdata, fdata, NCOEF,
NXVAL, constraint, IMSL_KNOTS, xknot, 0);
/* Compute
the average error of 100 points in the interval */
errlsq = 0.0;
errnft = 0.0;
for (i =
0; i < 100; i++) {
x = grdsiz *
(float) (i) / 99.0 - grdsiz/2.0;
errnft += fabs (F1 (x) -
imsl_f_spline_value(x,sp,0));
errlsq += fabs (F1 (x) - imsl_f_spline_value(x,spls,0));
}
/* Print results */
printf ("
Average error with BSLSQ fit: %8.5f\n",
errlsq / 100.0);
printf (" Average error with CONFT
fit:
%8.5f\n",
errnft / 100.0);
}
Output
Average error with BSLSQ fit:
0.01783
Average error with CONFT fit: 0.01339
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |