Let
(
x
i
,
Y
i
)
;
x
1
<
x
2
<
⋯
<
x
n
,
i
∈
Z
be a sequence of observations, modeled by the relation
Y
i
=
f
(
x
i
)
. The smoothing spline estimate
f
^
of the function
f
is defined to be the minimizer (over the class of twice differentiable functions) of
∑
i
=
1
n
(
Y
i
−
f
^
(
x
i
)
)
2
+
λ
∫
x
1
x
n
f
^
″
(
x
)
2
d
x
.
Remarks:
λ
≥
0
is a smoothing parameter, controlling the trade-off between fidelity to the data and roughness of the function estimate.
The integral is evaluated over the range of the
x
i
.
As
λ
→
0
(no smoothing), the smoothing spline converges to the interpolating spline.
As
λ
→
∞
(infinite smoothing), the roughness penalty becomes paramount and the estimate converges to a linear least squares estimate.
The roughness penalty based on the second derivative is the most common in modern statistics literature, although the method can easily be adapted to penalties based on other derivatives.
In early literature, with equally-spaced
x
i
, second or third-order differences were used in the penalty, rather than derivatives.
When the sum-of-squares term is replaced by a log-likelihood, the resulting estimate is termed penalized likelihood. The smoothing spline is the special case of penalized likelihood resulting from a Gaussian likelihood.
It is useful to think of fitting a smoothing spline in two steps:
- First, derive the values
f
^
(
x
i
)
;
i
=
1
,
…
,
n
.
- From these values, derive
f
^
(
x
)
for all x.
Now, treat the second step first.
Given the vector
m
^
=
(
f
^
(
x
1
)
,
…
,
f
^
(
x
n
)
)
T
of fitted values, the sum-of-squares part of the spline criterion is fixed. It remains only to minimize
∫
f
^
″
(
x
)
2
d
x
, and the minimizer is a natural cubic spline that interpolates the points
(
x
i
,
f
^
(
x
i
)
)
. This interpolating spline is a linear operator, and can be written in the form
f
^
(
x
)
=
∑
i
=
1
n
f
^
(
x
i
)
f
i
(
x
)
where
f
i
(
x
)
are a set of spline basis functions. As a result, the roughness penalty has the form
∫
f
^
″
(
x
)
2
d
x
=
m
^
T
A
m
^
.
where the elements of A are
∫
f
i
″
(
x
)
f
j
″
(
x
)
d
x
. The basis functions, and hence the matrix A, depend on the configuration of the predictor variables
x
i
, but not on the responses
Y
i
or
m
^
.
Now back to the first step. The penalized sum-of-squares can be written as
∥
Y
−
m
^
∥
2
+
λ
m
^
T
A
m
^
,
where
Y
=
(
Y
1
,
…
,
Y
n
)
T
. Minimizing over
m
^
gives
m
^
=
(
I
+
λ
A
)
−
1
Y
.
De Boor's approach exploits the same idea, of finding a balance between having a smooth curve and being close to the given data.
p
∑
i
=
1
n
(
Y
i
−
f
^
(
x
i
)
δ
i
)
2
+
(
1
−
p
)
∫
(
f
^
(
m
)
(
x
)
)
2
d
x
where
p
is a parameter called smooth factor and belongs to the interval
[
0
,
1
]
, and
δ
i
;
i
=
1
,
…
,
n
are the quantities controlling the extent of smoothing (they represent the weight
δ
i
−
2
of each point
Y
i
). In practice, since cubic splines are mostly used,
m
is usually
2
. The solution for
m
=
2
was proposed by Reinsch in 1967. For
m
=
2
, when
p
approaches
1
,
f
^
converges to the "natural" spline interpolant to the given data. As
p
approaches
0
,
f
^
converges to a straight line (the smoothest curve). Since finding a suitable value of
p
is a task of trial and error, a redundant constant
S
was introduced for convenience.
S
is used to numerically determine the value of
p
so that the function
f
^
meets the following condition:
∑
i
=
1
n
(
Y
i
−
f
^
(
x
i
)
δ
i
)
2
≤
S
The algorithm described by de Boor starts with
p
=
0
and increases
p
until the condition is met. If
δ
i
is an estimation of the standard deviation for
Y
i
, the constant
S
is recommended to be chosen in the interval
[
n
−
2
n
,
n
+
2
n
]
. Having
S
=
0
means the solution is the "natural" spline interpolant. Increasing
S
means we obtain a smoother curve by getting farther from the given data.
Given the constraint from the definition formula
x
1
<
x
2
<
⋯
<
x
n
we can conclude that the algorithm doesn't work for all sets of data. If we plan to use this algorithm for random points in a multidimensional space, to find a solution we need to give, as input to the algorithm, sets of data where these constraints are met. A solution for this is to introduce a parameter so that the input data would be represented as single-valued functions depending on that parameter; after this the smoothing will be performed for each function. In a bidimensional space a solution would be to parametrize
x
and
y
so that they would become
x
(
t
)
and
y
(
t
)
where
t
1
<
t
2
<
⋯
<
t
n
. A convenient solution for
t
is the cumulating distance
t
i
+
1
=
t
i
+
(
x
i
+
1
−
x
i
)
2
+
(
y
i
+
1
−
y
i
)
2
where
t
1
=
0
.
A more detailed analysis on parametrization is done by E.T.Y. Lee.
Smoothing splines are related to, but distinct from:
Regression splines. In this method, the data is fitted to a set of spline basis functions with a reduced set of knots, typically by least squares. No roughness penalty is used. (See also multivariate adaptive regression splines.)
Penalized Splines. This combines the reduced knots of regression splines, with the roughness penalty of smoothing splines.
Elastic maps method for manifold learning. This method combines the least squares penalty for approximation error with the bending and stretching penalty of the approximating manifold and uses the coarse discretization of the optimization problem; see thin plate splines.
Source code for spline smoothing can be found in the examples from Carl de Boor's book A Practical Guide to Splines. The examples are in the Fortran programming language. The updated sources are available also on Carl de Boor's official site [1].
