Least-squares spline approximation

`spap2(knots,k,x,y) `

spap2(l,k,x,y)

sp = spap2(...,x,y,w)

spap2({knorl1,...,knorlm},k,{x1,...,xm},y)

spap2({knorl1,...,knorlm},k,{x1,...,xm},y,w)

`spap2(knots,k,x,y) `

returns
the B-form of the spline *f* of order `k`

with
the given knot sequence `knots`

for which

(*) y(:,j) = f(x(j)), all j

in the weighted mean-square sense, meaning that the sum

$$\sum _{j}w(j)|y(:,j)-f\left(x(j)\right){|}^{2}$$

is minimized, with default weights equal to 1.
The data values `y(:,j)`

may be scalars, vectors,
matrices, even ND-arrays, and |*z*|^{2} stands
for the sum of the squares of all the entries of *z*.
Data points with the same site are replaced by their average.

If the sites `x`

satisfy the (Schoenberg-Whitney)
conditions

$$\begin{array}{l}\text{knots}(j)x(j)\text{knots}(j+k)\\ (**)\text{}j=1,\mathrm{...},\text{length}(x)=\text{length(knots)}-k\end{array}$$

then there is a unique spline (of the given order and knot sequence)
satisfying (*) exactly. No spline is returned unless (**) is satisfied
for some subsequence of `x`

.

`spap2(l,k,x,y) `

, with `l`

a positive
integer, returns the B-form of a least-squares spline approximant, but with the knot
sequence chosen for you. The knot sequence is obtained by applying `aptknt`

to an appropriate
subsequence of `x`

. The resulting piecewise-polynomial consists of
`l`

polynomial pieces and has `k-2`

continuous
derivatives. If you feel that a different distribution of the interior knots might
do a better job, follow this up with

sp1 = spap2(newknt(sp),k,x,y));

`sp = spap2(...,x,y,w) `

lets you specify the weights `w`

in the error measure
(given above). `w`

must be a vector of the same size
as `x`

, with nonnegative entries. All the weights
corresponding to data points with the same site are summed when those
data points are replaced by their average.

`spap2({knorl1,...,knorlm},k,{x1,...,xm},y) `

provides a least-squares spline approximation to *gridded* data.
Here, each `knorli`

is either a knot sequence or
a positive integer. Further, `k`

must be an `m`

-vector,
and `y`

must be an (`r+m`

)-dimensional
array, with `y(:,i1,...,im)`

the datum to be fitted
at the `site`

`[x{1}(i1),...,x{m}(im)]`

,
all `i1`

, ..., `im`

. However, if
the spline is to be scalar-valued, then, in contrast to the univariate
case, `y`

is permitted to be an `m`

-dimensional
array, in which case `y(i1,...,im)`

is the datum
to be fitted at the `site`

`[x{1}(i1),...,x{m}(im)]`

,
all `i1`

, ..., `im`

.

`spap2({knorl1,...,knorlm},k,{x1,...,xm},y,w) `

also lets you specify the weights. In this `m`

-variate
case, `w`

must be a cell array with `m`

entries,
with `w{i}`

a nonnegative vector of the same size
as `xi`

, or else `w{i}`

must be
empty, in which case the default weights are used in the `i`

th
variable.

sp = spap2(augknt([a,xi,b],4),4,x,y)

is the least-squares approximant to the data `x`

, `y`

,
by cubic splines with two continuous derivatives, basic interval
[`a`

..`b`

], and interior breaks `xi`

,
provided `xi`

has all its entries in `(a..b)`

and
the conditions (**) are satisfied in some fashion. In that case, the
approximant consists of `length(xi)+1`

polynomial
pieces. If you do not want to worry about the conditions (**) but
merely want to get a cubic spline approximant consisting of `l`

polynomial
pieces, use instead

sp = spap2(l,4,x,y);

If the resulting approximation is not satisfactory, try using
a larger `l`

. Else use

sp = spap2(newknt(sp),4,x,y);

for a possibly better distribution of the knot sequence. In fact, if that helps, repeating it may help even more.

As another example, `spap2(1, 2, x, y); `

provides
the least-squares straight-line fit to data `x`

,`y`

,
while

w = ones(size(x)); w([1 end]) = 100; spap2(1,2, x,y,w);

forces that fit to come very close to the first data point and to the last.

`spcol`

is called
on to provide the almost block-diagonal collocation matrix (*B _{j}*,

`slvblk`

solves the linear system (*)
in the (weighted) least-squares sense, using a block QR factorization. Gridded data are fitted, in tensor-product fashion, one variable at a time, taking advantage of the fact that a univariate weighted least-squares fit depends linearly on the values being fitted.