Thin-plate smoothing spline

is the stform of a thin-plate smoothing spline `st`

= tpaps(`x`

,`y`

)*f* for the
given data sites `x(:,j)`

and the given data values
`y(:,j)`

. The `x(:,j)`

must be distinct
points in the plane, the values can be scalars, vectors, matrices, even
ND-arrays, and there must be exactly as many values as there are sites.

The thin-plate smoothing spline *f* is the unique minimizer
of the weighted sum

$$pE(f)+(1-p)R(f)$$

with *E*(*f*) the error measure

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

and *R*(*f*) the roughness measure

$$R(f)={{\displaystyle \int (\left|{D}_{1}{D}_{1}f\right|}}^{2}+2{\left|{D}_{1}{D}_{2}f\right|}^{2}+{\left|{D}_{2}{D}_{2}f\right|}^{2})$$

Here, the integral is taken over all of
*R*^{2},
|*z*|^{2} denotes the sum of
squares of all the entries of *z*, and
*D _{i}f* denotes the partial
derivative of

`p`

so that `(1-p)/p`

equals the average
of the diagonal entries of the matrix `A`

, with ```
A +
(1-p)/p*eye(n)
```

the coefficient matrix of the linear system for the
`n`

coefficients of the smoothing spline to be determined.
This ensures staying in between the two extremes of interpolation (when
`p`

is close to `1`

and the coefficient
matrix is essentially `A`

) and complete smoothing (when
`p`

is close to `0`

and the coefficient
matrix is essentially a multiple of the identity matrix). This serves as a good
first guess for `p`

.

also inputs the `st`

= tpaps(`x`

,`y`

,`p`

)*smoothing parameter*, `p`

,
a number between 0 and 1. As the smoothing parameter varies from 0 to 1,
the smoothing spline varies, from the least-squares approximation to the data by
a linear polynomial when `p`

is `0`

, to the
thin-plate spline interpolant to the data when `p`

is
`1`

.

`[...,`

also returns the value of the smoothing parameter used in the final spline
result whether or not you specify `P`

] = tpaps(...) `p`

. This syntax is useful
for experimentation in which you can start with ```
[pp,P] =
tpaps(x,y)
```

and obtain a reasonable first guess for
`p`

.

The determination of the smoothing spline involves the solution of a linear system
with as many unknowns as there are data points. Since the matrix of this linear
system is full, the solving can take a long time even if, as is the case here, an
iterative scheme is used when there are more than 728 data points. The convergence
speed of that iteration is strongly influenced by `p`

, and is
slower the larger `p`

is. So, for large problems, use
interpolation, i.e., `p`

equal to 1, only if you can afford
the time.