splineLsqConstrained¶

Computes a least-squares constrained spline approximation.

Synopsis¶

splineLsqConstrained (xdata, fdata, splineSpaceDim, numConPts, constraints)

Required Arguments¶

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 splineSpaceDim (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 numConPts (Input)

The number of points in the vector constraints.

structure 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	\(bl_i = f^{\left(d_i\right)} \left(x_i\right)\)
2	\(f^{\left(d_i\right)} \left(x_i\right) \leq bu_i\)
3	\(f^{\left(d_i\right)} \left(x_i\right) \geq bl_i\)
4	\(bl_i \leq f^{\left(d_i\right)} \left(x_i\right) \leq bu_i\)
5	\(bl_i = \int_{c}^{d} f(t) dt\)
6	\(\int_{c}^{d} f(t) dt \leq bu_i\)
7	\(\int_{c}^{d} f(t) dt \geq bl_i\)
8	\(bl_i \leq \int_{c}^{d} f(t) dt \leq bu_i\)
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	\(bl_i = f^{\left(d_i\right)} \left(x_i\right) - f^{\left(d_{i+1}\right)} \left(x_{i+1}\right)\)
10	\(f^{\left(d_i\right)} \left(x_i\right) \leq bu_i\)
11	\(f^{\left(d_i\right)} \left(x_i\right) - f^{\left(d_{i+1}\right)} \left(x_{i+1}\right) \geq bl_i\)
12	\(bli \leq f^{\left(d_i\right)} \left(x_i\right) - f^{\left(d_{i+1}\right)} \left(x_{i+1}\right) \leq bu_i\)

Return Value¶

The structure that represents the spline fit. If a fit cannot be computed, then None is returned.

Optional Arguments¶

nhard (Input)

The argument nhard is the number of entries of constraints involved in the “hard” constraints. Note that 0 ≤ nhard ≤ numConPts. The default, nhard = 0, always results in a fit, while setting nhard = numConPts 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

weights, float[] (Input)

This option requires the user to provide the weights.

Default: all weights equal one

order, int (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)

knots, float[] (Input)

This option requires the user to provide the knots. The user must provide a knot sequence of length splineSpaceDim + order.

Default: an appropriate knot sequence is selected. See below for more details.

Description¶

The function splineLsqConstrained 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:

\(E_p\)[f]	\(= f^{\left(j_p\right)} \left(y_p\right)\)
Or	\(= f^{\left(j_p\right)} \left(y_p\right) - f^{\left(j_{p+1}\right)} \left(y_{p+1}\right)\)
Or	\(= \int_{y_p}^{y_{p+1}} f(t) dt\)
Or	= periodic end conditions

An interval, \(I_p\) (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 \((x_i, f_i)\) 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 \(n_f\) denote the number of feasible constraints as described above. Then, the function solved the problem

\[\sum_{i=1}^{n} \left|f_i - \sum_{j=1}^{m} a_jB_j\left(x_i\right)\right|^2w_i\]

subject to

\[E_p \left[\sum_{j=1}^{m} a_jB_j\right] \in I_p \phantom{...} p=1, \dots, n_f\]

This linearly constrained least-squares problem is treated as a quadratic program and is solved by invoking the function quadraticProg (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 \(I_p\).

Examples¶

Example 1¶

This is a simple application of splineLsqConstrained. Data is generated from the function

\[\frac{x}{2} + \sin \left(\frac{x}{2}\right)\]

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 splineLeastSquares. Then, the derivative is forced to be greater than 0 at numConPts = 15 equally spaced points and splineLsqConstrained is called. The resulting curve is monotone. The error is printed for the two fits averaged over 100 equally spaced points.

from __future__ import print_function
from numpy import *
from pyimsl.math.mathStructs import d_constraint_struct
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform
from pyimsl.math.splineLeastSquares import splineLeastSquares
from pyimsl.math.splineLsqConstrained import splineLsqConstrained
from pyimsl.math.splineValue import splineValue


def F1(x):
    return .5 * x + sin(.5 * x)


korder = 4
ndata = 15
nxval = 15
ncoef = 8
constraint_ctype = d_constraint_struct * nxval
constraint = constraint_ctype()
fdata = empty(ndata)
xdata = empty(ndata)

# Compute original xdata and fdata with random noise
randomSeedSet(234579)
noise = randomUniform(ndata)
grdsiz = 10.0
for i in range(0, ndata):
    xdata[i] = grdsiz * (float(i) / (ndata - 1))
    fdata[i] = F1(xdata[i]) + (noise[i] - .5)

# Compute least-squares fit
spls = splineLeastSquares(xdata, fdata, ncoef)

# Construct the constraints.
for i in range(0, nxval):
    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 = splineLsqConstrained(xdata, fdata, ncoef,
                          nxval, constraint)

# Compute the average error of 100 points in the interval
errlsq = 0.0
errnft = 0.0
for i in range(0, 100):
    x = grdsiz * float(i) / 99.0
    errnft = errnft + abs(F1(x) - splineValue(x, sp))
    errlsq = errlsq + abs(F1(x) - splineValue(x, spls))

# Print results
print(" Average error with splineLeastSquares fit:   %8.5f"
      % (errlsq / 100.0))
print(" Average error with splineLsqConstrained fit: %8.5f"
      % (errnft / 100.0))

Output¶

 Average error with splineLeastSquares fit:    0.20250
 Average error with splineLsqConstrained fit:  0.14334

Example 2¶

Now, try to recover the function

\[\frac{1}{1+x^4}\]

from noisy data. First, try the unconstrained least-squares fit using splineLeastSquares. Finding that fit somewhat unsatisfactory, several constraints are applied using splineLsqConstrained. 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 \int_{-7}^{7} s(x)dx \le 2.3\]

\[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.

from __future__ import print_function
from numpy import *
from pyimsl.math.mathStructs import d_constraint_struct
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform
from pyimsl.math.splineLeastSquares import splineLeastSquares
from pyimsl.math.splineLsqConstrained import splineLsqConstrained
from pyimsl.math.splineValue import splineValue


def F1(x):
    return 1.0 / (1.0 + x * x * x * x)


# Compute original xdata and fdata with random noise
korder = 4
ndata = 51
nxval = 12
ncoef = 13
constraint_ctype = d_constraint_struct * nxval
constraint = constraint_ctype()
fdata = empty(ndata)
xdata = empty(ndata)
xknot = empty(ncoef + korder)

randomSeedSet(234579)
noise = randomUniform(ndata)
grdsiz = 14.0
for i in range(0, ndata):
    xdata[i] = grdsiz * (float(i) / (ndata - 1)) - grdsiz / 2.0
    fdata[i] = F1(xdata[i]) + 0.125 * (noise[i] - .5)

# Generate knots
for i in range(0, ncoef - korder + 2):
    xknot[i + korder - 1] = grdsiz * \
        (float(i) / (ncoef - korder + 1)) - grdsiz / 2.0
for i in range(0, korder - 1):
    xknot[i] = xknot[korder - 1]
    xknot[i + ncoef + 1] = xknot[ncoef]

# Compute splineLeastSquares fit
spls = splineLeastSquares(xdata, fdata, ncoef, knots=xknot)

# Construct the constraints for CONFT
for i in range(0, 4):
    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 = splineLsqConstrained(xdata, fdata, ncoef,
                          nxval, constraint, knots=xknot)

# Compute the average error of 100 points in the interval
errlsq = 0.0
errnft = 0.0
for i in range(0, 100):
    x = grdsiz * float(i) / 99.0 - grdsiz / 2.0
    errnft = errnft + abs(F1(x) - splineValue(x, sp))
    errlsq = errlsq + abs(F1(x) - splineValue(x, spls))

# Print results
print(" Average error with BSLSQ fit: %8.5f" % (errlsq / 100.0))
print(" Average error with CONFT fit: %8.5f" % (errnft / 100.0))

Output¶

 Average error with BSLSQ fit:  0.01783
 Average error with CONFT fit:  0.01339