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 imsl_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 (Input)
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
int 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
int 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].type = 1;
constraint[i+4].type = 1;
constraint[i].der = 2;
constraint[i+4].der = 2;
constraint[i].bl = 0.0;
constraint[i+4].bl = 0.0;
}
constraint[8].xval = -7.0;
constraint[8].type = 3;
constraint[8].der = 0;
constraint[8].bl = 0.0;
constraint[9].xval = -7.0;
constraint[9].type = 6;
constraint[9].bu = 2.3;
constraint[10].xval = 7.0;
constraint[10].type = 6;
constraint[10].bu = 2.3;
constraint[11].xval = -7.0;
constraint[11].type = 20;
constraint[11].der = 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
