# Comparative Analysis for Robust Penalized  Spline Smoothing Methods

Authors: Bin Wang; Wenzhong Shi; Zelang Miao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Smoothing noisy data is commonly encountered in engineering domain, and currently robust penalized regression spline models are perceived to be the most promising methods for coping with this issue, due to their flexibilities in capturing the nonlinear trends in the data and effectively alleviating the disturbance from the outliers. Against such a background, this paper conducts a thoroughly comparative analysis of two popular robust smoothing techniques, theM-type estimator and S-estimation for penalized regression splines, both of which are reelaborated starting from their origins, with their derivation process reformulated and the corresponding algorithms reorganized under a unified framework. Performances of these two estimators are thoroughly evaluated from the aspects of fitting accuracy, robustness, and execution time upon the MATLAB platform. Elaborately comparative experiments demonstrate that robust penalized spline smoothing methods possess the capability of resistance to the noise effect compared with the nonrobust penalized LS spline regression method. Furthermore, the M-estimator exerts stable performance only for the observations with moderate perturbation error, whereas the S-estimator behaves fairly well even for heavily contaminated observations, but consuming more execution time. These findings can be served as guidance to the selection of appropriate approach for smoothing the noisy data.


---

## Body

## 1. Introduction

Penalized spline smoothing technique has already been perceived as a popular nonparametric smoothing approach for smoothing noisy data in engineering domain [1–3] during the past 20 years due to its ease of fitting, flexible choice of the inner knots, and smoothing parameter. The ideas of employing the combination of regression splines, with substantially smaller number of inner knots than the sample size, and the introduction of the penalized parameter, which compromises the fitting accuracy and the smoothness of the regression curve, can be traced back to at least O’Sullivan [4] who utilized a cubic B-spline basis for estimation in ill-posed inverse problems.During the early stage, Kelly and Rice [5] and Besse et al. [6] approximated the smoothing splines by the hybrid splines equipped with the inner knots equal to the data samples and a penalty parameter for determining the smoothing amount of regression curve, which was usually calculated according to the (generalized) cross-validation criterion. Other criterions included the Mallows’ C
p criterion [7] and the Akaike information criterion [8, 9]. Besides, Eilers and Marx [10] proposed flexible criteria for choosing an optimal penalty parameter of dealing with the B-splines coefficients. Wand [11] additionally derived a quick and simple formula approximation to the optimal penalty parameter for the penalized spline regression. Moreover, Ruppert and Carroll [12] investigated the spatially varying penalties and Ruppert [13] recommended the selection of designed inner knots. More discussion and examples on the penalized regression spline models can be referred to the excellent monograph written by Ruppert et al. [14].The subsequent explorations on penalized spline smoothing remained a hot area arousing great interests from both theory researchers and practice engineers for several years. For example, the theoretical aspects and practical applications of penalized spline smoothing were fully discussed in a dissertation composed by Krivobokova [15] from financial mathematics. Crainiceanu et al. [16] investigated the implementation of nonparametric Bayesian analysis for penalized spline regression using WinBUGS, in which the Markov Chain Monte Carlo (MCMC) mixing operator was substantially improved by employing the low rank thin plate splines basis. Besides, Marra and Radice [2] presented an accessible overview of generalized additive models based on the penalized regression splines and explored their potential applications into medical research. Generally, the penalized regression spline model is referred to the estimation of univariate smooth functions based on noisy data. However, some researchers also put forward its extension into higher dimension. Typical models include the thin plate regression splines [17] which are low rank smoothers obtained by truncating the basis of standard thin plate splines thereby requiring less computation to fit to large data sets and the penalized bivariate splines [18] with application on triangulations for irregular regions.The regression curve for penalized spline is established by minimizing the sum of squared residuals subject to a bound on the size of the spline coefficients. Accordingly, a closed-form formula can be derived from the penalized minimization problem. Apparently, such penalized least squares regression splines cannot be resistant to the disturbance from the outliers. A simple and direct idea is to replace the squared residuals by some slower-increasing loss function similarly as employed inM-regression estimator [19, 20] for alleviating the outlier effects from the atypical observations. Early research regarding the M-type smoothing technique started from Huber et al. [21] and Cox’s [22] work on cubic regression splines. Another robust model can be the S-estimation [23] which possesses high-breakdown point and minimizes the residual scale in an extremely robust manner. A fast algorithm has been created by Salibian-Barrera and Yohai [24] for computing the S-regression estimates. These two estimators have been elaborately compared in the linear regression by ÇetIn and Toka [25] with the results indicating that the S-estimator is more capable of resisting the outlier disturbance compared with the M-estimator but with lower efficiency.Up until almost a decade ago, however, little work has been delivered taking the penalized spline regression into the category of robust estimation. Oh et al. [26] suggested a robust smoothing approach implemented to the period analysis of variable stars by simply imposing Huber’s favorite loss function upon the generalized cross-validation criterion. And Lee and Oh [27] proposed an iterative algorithm for calculating the M-type estimator of penalized regression splines via introducing the empirical pseudodata. Afterwards, Finger [28] confirmed that such M-type estimator compromised between the very robust estimators and the very efficient least square type estimators for the penalized spline regressions when investigating the performance of different estimation techniques in crop insurance applications.   Meanwhile, by replacing the least squares estimation with a suitable S-estimator, Tharmaratnam et al. [29] put forward an S-estimation for penalized regression spline as well, which proves to be equivalent to a weighted penalized least squares regression and behaves well even for heavily contaminated observations.Against such a background, this paper performs a thoroughly comparative analysis mainly based on these two popular robust smoothing techniques, theM-type estimator and S-estimation for penalized regression splines, both of which are reelaborated starting from their origins, with their derivation process reformulated and the corresponding algorithms reorganized under a unified framework. The rest of the paper is organized as follows. Section 2 reelaborates the M-type estimator and S-estimation for penalized regression splines, starting from the least squares estimator and reorganizes the corresponding algorithms within a unified framework. Section 3 describes the comparative experiments upon both simulated synthetic data and one real weather balloon data set. Concluding remarks are included in the end of Section 4.


## 2. Robust Penalized Regression Splines

The following discussion is devoted to the reelaboration of penalized regression splines starting from the origin, with their derivation processes reformulated and the corresponding practical algorithms reorganized under a unified framework.The idea of the penalized regression spline begins with the following model:(1)
y
i
~
N
(
m
(
x
i
)
,
σ
ε
2
)
,
i
=
1
,
…
,
n
,
where m
(
x
) is supposed to be a smooth but unknown regression function which needs to be estimated based on the sample observations (
x
i
,
y
i
), i
=
1
,
…
,
n and σ
ε
2 denotes the constant variance of the random deviation error between the response variable and the regression function m
(
x
).It is a flexible concept as regards the penalized spline smoothing since different basis function may correspond to different penalized spline smoothing regression function. A common selection is the truncated polynomial bases, but other choices can also be explored straightway similarly(2)
m
(
x
;
β
)
=
β
0
+
β
1
x
+
⋯
+
β
p
x
p
+
∑
j
=
1
K
β
p
+
j
(
x
-
κ
j
)
+
p
,
where (
x
-
κ
j
)
+
=
max
⁡
⁡
(
x
-
κ
j
,
0
) and β
=
(
β
0
,
β
1
,
…
,
β
p
+
K
)
T denotes a vector of regression coefficients, κ
j
,
…
,
κ
K is the specified inner knots, and p indicates the exponential order for truncated power basis.


### 2.1. Penalized Least Squares Regression Splines

Given a group of sample observations(
x
1
,
y
1
)
,
…
,
(
x
n
,
y
n
), an increasingly popular way for obtaining the estimation of m
(
·
) is via transforming it into the category of a least squares problem, where we need to seek out the member from the class m
(
x
;
β
) which minimizes the sum of squared residuals. To prevent overfitting, a constraint on the spline coefficients is always imposed; that is,
(3)
min
⁡
β
∈
R
p
+
K
+
1
⁡
∑
i
=
1
n
(
y
i
-
m
(
x
i
,
β
)
)
2
,
subject
to
∑
j
=
1
K
β
p
+
j
2
≤
C
,
for some constant C
>
0 as discussed in Ruppert et al. [14]. Thus, by employing the Lagrange multiplier, the coefficients of the penalized least squares regression spline β
^
LS are equivalent to the minimizer of
(4)
∑
i
=
1
n
(
y
i
-
m
(
x
i
;
β
)
)
2
+
λ
∑
j
=
1
K
β
p
+
j
2
,
for some smoothing (penalty) parameter λ
≥
0 with respect to β, which compromises the bias and variance tradeoff of the fitted regression spline curve.IntroduceF
(
x
)
=
(
1
,
x
,
…
,
x
p
,
(
x
-
κ
1
)
+
p
,
…
,
(
x
-
k
K
)
+
p
)
T
, the spline design matrix
(5)
X
=
(
F
(
x
1
)
T
,
…
,
F
(
x
n
)
T
)
T
=
(
1
x
1
⋯
x
1
p
(
x
1
-
κ
1
)
+
p
⋯
(
x
1
-
κ
K
)
+
p
⋮
⋮
⋱
⋮
⋮
⋱
⋮
1
x
n
⋯
x
n
p
(
x
n
-
κ
1
)
+
p
⋯
(
x
n
-
κ
K
)
+
p
)
,
the vector of responses y
=
(
y
1
,
…
,
y
n
)
T, and the diagonal matrix D
=
diag
⁡
(
0
(
p
+
1
)
×
1
,
1
K
×
1
) which indicates that only the spline coefficients are penalized. Thus, the closed-form minimizer for the objective function (4) can arrive by direct calculations; namely,
(6)
β
^
LS
=
(
X
T
X
+
λ
D
)
-
1
X
T
y
.
Consequently, the corresponding estimation vector m
^
=
(
m
^
(
x
1
)
,
…
,
m
^
(
x
n
)
)
T generated by the penalized least square regression spline can be represented as
(7)
m
^
LS
=
X
β
^
LS
=
X
(
X
T
X
+
λ
D
)
-
1
X
T
y
.
For figuring out the penalty parameter λ in (7), the (generalized) cross-validation technique is employed hereafter, which is computed by leave-one-out of the residual from sum of squares so as to avoid overfitting in the regression spline. The cross-validation criterion is given by
(8)
CV
(
λ
)
=
1
n
∑
i
=
1
n
(
y
i
-
m
^
-
i
(
x
i
;
λ
)
)
2
,
where m
^
-
i
(
x
i
;
λ
) is the regression spline estimated by leaving out the ith observation point. The optimized penalty parameter λ can be calculated by minimizing CV
(
λ
) over λ
≥
0. For the linear smoothing matrix, H
λ
=
X
(
X
T
X
+
λ
D
)
-
1
X
T, and it can be verified [14] that
(9)
m
^
-
i
(
x
i
;
λ
)
=
∑
j
≠
i
(
H
λ
i
j
1
-
H
λ
i
j
y
i
)
.
Substitution of m
^
-
i
(
x
i
;
λ
) in (8) with the formula (9) results in
(10)
CV
(
λ
)
=
1
n
∑
i
=
1
n
(
y
i
-
m
^
(
x
i
;
λ
)
1
-
H
λ
i
i
)
2
.
The generalized cross-validation criterion is proposed [30] simply by replacing H
λ
i
i in (10) with their average tr
⁡
(
H
λ
)
/
n; namely,
(11)
GCV
(
λ
)
=
n
∑
i
=
1
n
(
y
i
-
m
^
(
x
i
;
λ
)
n
-
tr
⁡
(
H
λ
)
)
2
.


### 2.2.M-Type Estimator for Penalized Regression Splines

The above estimatem
^
LS of the penalized least squares regression splines may greatly suffer from various robustness issues especially when the data is contaminated with outliers. For alleviating this effect, a penalized regression estimator can be straightforwardly constructed by replacing the previous squared residual loss function with the following M-type criterion [27]:
(12)
∑
i
=
1
n
ρ
(
y
i
-
m
(
x
i
;
β
)
)
2
+
λ
∑
j
=
1
K
β
p
+
j
2
,
where ρ is even, nondecreasing in [
0
,
+
∞
) and ρ
(
0
)
=
0, for which a common choice is taken of the Huber loss function with cutoff c
>
0,
(13)
ρ
c
(
t
)
=
{
t
2
|
t
|
≤
c
2
c
|
t
|
-
c
2
|
t
|
>
c
.
Apparently, the above Huber’s function is a parabola in the vicinity of zero, and it increases linearly at a given level |
t
|
>
c. A default choice of the tuning constant c
=
1.345 aims for a 95% asymptotic efficiency with respect to the standard normal distribution.Denote the derivative ofρ
c as ψ
c
(
s
)
=
max
⁡
⁡
[
-
c
,
min
⁡
⁡
(
c
,
s
)
]. When a set of fitted residuals r
i
=
y
i
-
m
(
x
i
;
β
), i
=
1
,
…
,
n from the penalized spline regression are given, a robust M-scale estimator [19, 20] can be employed for estimating the standard deviation σ
ε of these residuals as follows:
(14)
∑
i
=
1
n
ψ
c
(
y
i
-
m
(
x
i
;
β
)
σ
^
ε
)
=
0
.
A common choice for a measure of scale with a high finite sample breakdown point can be the value of ω determined by
(15)
P
(
|
r
-
r
0.5
|
<
ω
)
=
1
2
,
where r
0.5 is the population median. Besides, the standard estimation of ω usually employs the median absolute deviation (MAD) statistic, which is defined as
(16)
MAD
=
Median
{
|
r
1
-
M
r
|
,
…
,
|
r
n
-
M
r
|
}
,
where M
r is the usual sample median of the fitted residuals and MAD is actually the sample median of the n values |
r
1
-
M
r
|
,
…
,
|
r
n
-
M
r
|, with its finite sample breakdown point approximately being 0.5 [20]. Suppose the fitted residuals are randomly sampled from a normal distribution; note that MAD does not estimate the corresponding standard deviation, σ
^
ε, but rather estimates Z
0.75
σ
^
ε (with Z
0.75 being the 0.75 quantile of the standard normal distribution). To put MAD in a more familiar context, it typically needs to be rescaled so as to estimate σ
^
ε especially when the residuals are sampled from a normal distribution. In particular,
(17)
MADN
=
MAD
z
0.75
≈
1.4836
MAD
.
By utilizing the aforementioned spline basis functions, an M-type penalized spline estimator using the standardized residuals can be computed as m
^
M
=
X
β
^
M, with its regression coefficients being as follows:
(18)
β
^
M
=
argmin
β
⁡
∑
i
=
1
n
ρ
c
(
y
i
-
m
(
x
i
;
β
)
σ
^
ε
)
+
λ
β
T
D
β
.
Due to the nonlinear nature of ρ
c and ψ
c, it is nontrivial to seek the minimizer of (18) while satisfying (14). Under such circumstances, an M-type iterative algorithm for calculating the penalized regression splines is proposed [27] by introducing the empirical pseudodata
(19)
z
=
m
^
p
+
ψ
c
(
y
-
m
^
p
)
2
,
where m
^
p is the penalized least square estimator corresponding to the above pseudodata
(20)
m
^
p
=
X
(
X
T
X
+
λ
D
)
-
1
X
T
z
.
As proved in Theorem 1 of Lee and Oh’s [27] work, the penalized least square estimator m
^
p converges to m
^
M asymptotically, which implies an M-type iterative algorithm consisting of the steps in Algorithm 1.Algorithm 1:
M-type iterative algorithm for penalized regression splines.
Step  0.  Input an initial curve estimate m
^
p
(
0
), termination tolerance ε, and the maximum iteration
number Itermax. Meanwhile, set k = 0 and conduct the following loop iterations.
Step  1.  Estimate σ
^
ε by using MADN in (17);
Step  2.  Produce the empirical pseudo data according to (19);
Step  3.  Calculate the penalized least-square estimator m
^
p
(
k
+
1
) defined in (20) for the pseudo data,
where the penalty parameter λ is chosen according to the generalized cross-validation
(GCV) criterion formulated by (11) in the foregoing subsection.
Step  4.  If ∥
m
^
p
(
k
+
1
)
-
m
^
p
(
k
)
∥
<
ε
∥
m
^
p
(
k
)
∥ or k = Itermax, terminate and output the robust
M-type penalized spline estimate m
^
M = m
^
p
(
k
+
1
); else set k = k + 1 and continue Step  1.Regarding theM-type estimator for the penalized regression splines (see Algorithm 1), the initial curve estimate m
^
p
(
0
) in Step 0 is suggested [27], choosing the nonrobust penalized least square regression spline presented by (7).


### 2.3.S-Estimation for Penalized Regression Splines

The robustness properties of the aforementionedM-type estimator for penalized regression splines are not completely satisfactory since they have low breakdown point [31], which was clarified by Maronna and Yohai [32] for the general M-estimators when considering both the regression and scale, simultaneously.Another more robust and relatively latest approach for smoothing the noisy data is theS-estimation for penalized regression splines [29], whose core idea is utilizing a flexible S-estimator to replace the traditional least squares estimator.  The coefficient vector of the penalized S-regression splines β
^
S is defined as follows:
(21)
β
^
S
=
argmin
β
⁡
{
n
σ
^
ε
2
(
β
)
+
λ
β
T
D
β
}
,
in which, for each β, σ
^
ε
2
(
β
) satisfies
(22)
1
n
∑
i
=
1
n
ρ
(
y
i
-
F
(
x
i
)
T
β
σ
^
ε
(
β
)
)
=
E
Φ
(
ρ
(
Z
)
)
,
where Φ is cumulative distribution function of the standard norm distribution, Z
~
N
(
0,1
), and the loss function is taken from Tukey’s bisquare family [33]
(23)
ρ
d
(
u
)
=
{
3
(
u
d
)
2
-
3
(
u
d
)
4
+
(
u
d
)
6
|
u
|
≤
d
1
|
u
|
>
d
.
It is evidently known from theoretical assumption that each fitted residual vector should follow a normal distribution with mean of zero; that is, r
i
~
N
(
0
,
σ
ε
2
), i
=
1
,
…
,
n, however its median statistic, Median(
r
), in which r
=
(
r
1
,
…
,
r
n
)
T, may greatly depart from zero to a large extent in the general practice. Considering this, we ameliorate the estimation of σ
^
ε
(
β
) by using MADN in (17) with the addition of the absolute value of residual median; namely,
(24)
σ
^
ε
(
β
)
=
MADN
+
Median
(
r
(
β
)
)
.
The penalized S-estimator for the regression spline model is computed as m
^
S
=
X
β
^
S, with its coefficients vector β
^
S satisfying both (21) and (22), simultaneously. The following formula derivations indicate that such penalized S-estimator can actually be equivalent to a weighted penalized least squares regression.Before proceeding with the derivation, it is forewarned that the symbolβ
^
S is simply written as β just for clarity. Firstly, when (21) reaches its minimum, the following equation should be satisfied:
(25)
n
σ
^
ε
(
β
)
∇
σ
^
ε
(
β
)
+
λ
D
β
=
0
.
On the other hand, taking the derivative of the M-scale function in (22) with respect to β leads to
(26)
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
σ
^
ε
(
β
)
)
(
-
F
(
x
i
)
T
σ
^
ε
(
β
)
-
r
i
(
β
)
∇
σ
^
ε
(
β
)
T
σ
^
ε
2
(
β
)
)
=
0
,
in which r
i
(
β
)
=
y
i
-
F
(
x
i
)
T
β, r
~
i
(
β
)
=
r
i
(
β
)
/
σ
^
ε
(
β
), and Λ
(
β
)
=
diag
⁡
(
ρ
′
(
r
~
i
(
β
)
)
/
r
~
i
(
β
)
),  i
=
1
,
…
,
n.Deformation of the above formula produces(27)
∇
σ
^
ε
(
β
)
=
-
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
F
(
x
i
)
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
=
-
(
σ
^
ε
(
β
)
X
T
Λ
(
β
)
r
(
β
)
r
(
β
)
T
Λ
(
β
)
r
(
β
)
)
.
Substituting ∇
σ
^
ε
(
β
) of (27) into (25),
(28)
-
n
σ
^
ε
(
β
)
2
r
(
β
)
T
Λ
(
β
)
r
(
β
)
X
T
Λ
(
β
)
r
(
β
)
+
λ
D
β
=
0
.
Recall that r
(
β
)
=
y
-
X
β and import new symbol τ
(
β
)
=
n
σ
^
ε
(
β
)
2
/
r
(
β
)
T
Λ
(
β
)
r
(
β
); we have
(29)
-
τ
(
β
)
X
T
Λ
(
β
)
(
y
-
X
β
)
+
λ
D
β
=
0
.
And, thus, a further simple deformation arrives at the coefficient vector β
^
S with an iteration form
(30)
β
=
(
X
T
Λ
(
β
)
X
+
λ
τ
(
β
)
D
)
-
1
X
T
Λ
(
β
)
y
,
which forms the basis of constructing an iterative algorithm for penalized S-regression spline.The above formula derivations verify thatβ
^
S is actually the solution to the following weighted penalized least squares regression problem:
(31)
min
⁡
β
⁡
∥
Λ
(
β
)
1
/
2
(
y
-
X
β
)
∥
2
+
λ
τ
(
β
)
β
T
D
β
.
And such representation indicates that the GCV criterion in (11) can be employed straightforwardly with response variable y
~
=
Λ
(
β
)
1
/
2
y, weighted design matrix X
~
=
Λ
(
β
)
1
/
2
X, and the penalty term λ
/
τ
(
β
). Meanwhile, due to some weight elements being zero, the penalty parameter λ is thus calculated based on the following Regularized GCV criterion:
(32)
RGCV
(
λ
)
=
n
w
∥
Λ
(
β
)
1
/
2
(
y
-
X
β
)
∥
2
(
n
w
-
tr
⁡
(
H
~
λ
)
)
2
,
where H
~
λ
=
X
~
(
X
~
T
X
~
+
(
λ
/
τ
(
β
)
)
D
)
-
1
X
~
T
=
Λ
(
β
)
1
/
2
X
(
X
~
T
Λ
(
β
)
X
+
(
λ
/
τ
(
β
)
)
D
)
-
1
X
T
Λ
(
β
)
1
/
2 and n
w denotes the number of nonzero elements in weight matrix Λ
(
β
).Consequently, based on the foregoing in-depth discussions, theS-estimator iterative algorithm for penalized regression splines can be constructed as in Algorithm 2.Algorithm 2:
S-estimator iterative algorithm for penalized regression splines.
Step  0. Input an initial regression coefficients β
^
(
0
), the termination tolerance ε, and the maximum
iteration number Itermax. Meanwhile, set k = 0 and enter into the following loop iterations.
Step  1. Estimate σ
^
ε based on (24), that is
σ
^
ε
(
β
^
(
k
)
)
=
MADN
+
Median
(
r
(
β
^
(
k
)
)
).
Step  2. Compute the intermediate variables τ
(
β
^
(
k
)
) and Λ
(
β
^
(
k
)
).
Step  3. Calculate the coefficient vector β
^
(
k
+
1
) according to (30) as
β
^
(
k
+
1
)
=
(
X
T
Λ
(
β
^
(
k
)
)
X
+
λ
τ
(
β
^
(
k
)
)
D
)
-
1
X
T
Λ
(
β
^
(
k
)
)
y.
Hereinto, the penalty parameter λ is figured out based on the regularized GCV criterion (32).
Step  4. If ∥
β
^
(
k
+
1
)
-
β
^
(
k
)
∥
<
ε
∥
β
^
(
k
)
∥ or k = Itermax, terminate and output the convergent
coefficient vector for the penalized regression spline of S-estimator β
^
★
=
β
^
(
k
+
1
); else
accumulate the counter by setting k = k + 1 and repeat Step  1.In particular, it should be pointed out that since the employed loss functionρ
d is bounded giving rise to the nonconvexity of σ
^
ε
(
β
), the objective function in (21) may have several critical points, each of which only corresponds to one local minima. Therefore, the convergent critical point of the iterations derived from the above algorithm completely depends on the prescribed initial regression coefficients vector, β
^
(
0
), in Step 0 of Algorithm 2.One possible solution to cope with this issue proposed by Tharmaratnam et al. [29] is based on the random sampling theory, which can be briefly described as, firstly, creating J random subsamples with a proper size, for example, one fifth of the sample number; namely, n
s
=
max
⁡
⁡
(
p
+
K
+
1
,
⌊
n
/
5
⌋
), from the original sample observations. Each subsample brings about an initial candidate of the coefficient vector, for instance, assigned by least squares fitting for the penalized regression splines or the M-type estimator with a small amount of additional computation cost, with the initial regression coefficient vector of the jth subsample being denoted by β
^
j
(
0
), staring from which the final iterative convergent vector produced through Algorithm 2 is being written as β
^
j
⋆, j
=
1
,
…
,
J. Eventually, the global optimum is taken among all these J convergent vectors to be the one corresponding to the minimum objective function in (21); namely,
(33)
β
^
S
⋆
=
argmin
1
≤
j
≤
J
⁡
{
n
σ
^
ε
2
(
β
^
j
⋆
)
+
λ
β
^
j
⋆
T
D
β
^
j
⋆
}
.
Therefore, the robust S-estimator for penalized regression splines can be computed as m
^
S
=
X
β
^
S
⋆.


## 2.1. Penalized Least Squares Regression Splines

Given a group of sample observations(
x
1
,
y
1
)
,
…
,
(
x
n
,
y
n
), an increasingly popular way for obtaining the estimation of m
(
·
) is via transforming it into the category of a least squares problem, where we need to seek out the member from the class m
(
x
;
β
) which minimizes the sum of squared residuals. To prevent overfitting, a constraint on the spline coefficients is always imposed; that is,
(3)
min
⁡
β
∈
R
p
+
K
+
1
⁡
∑
i
=
1
n
(
y
i
-
m
(
x
i
,
β
)
)
2
,
subject
to
∑
j
=
1
K
β
p
+
j
2
≤
C
,
for some constant C
>
0 as discussed in Ruppert et al. [14]. Thus, by employing the Lagrange multiplier, the coefficients of the penalized least squares regression spline β
^
LS are equivalent to the minimizer of
(4)
∑
i
=
1
n
(
y
i
-
m
(
x
i
;
β
)
)
2
+
λ
∑
j
=
1
K
β
p
+
j
2
,
for some smoothing (penalty) parameter λ
≥
0 with respect to β, which compromises the bias and variance tradeoff of the fitted regression spline curve.IntroduceF
(
x
)
=
(
1
,
x
,
…
,
x
p
,
(
x
-
κ
1
)
+
p
,
…
,
(
x
-
k
K
)
+
p
)
T
, the spline design matrix
(5)
X
=
(
F
(
x
1
)
T
,
…
,
F
(
x
n
)
T
)
T
=
(
1
x
1
⋯
x
1
p
(
x
1
-
κ
1
)
+
p
⋯
(
x
1
-
κ
K
)
+
p
⋮
⋮
⋱
⋮
⋮
⋱
⋮
1
x
n
⋯
x
n
p
(
x
n
-
κ
1
)
+
p
⋯
(
x
n
-
κ
K
)
+
p
)
,
the vector of responses y
=
(
y
1
,
…
,
y
n
)
T, and the diagonal matrix D
=
diag
⁡
(
0
(
p
+
1
)
×
1
,
1
K
×
1
) which indicates that only the spline coefficients are penalized. Thus, the closed-form minimizer for the objective function (4) can arrive by direct calculations; namely,
(6)
β
^
LS
=
(
X
T
X
+
λ
D
)
-
1
X
T
y
.
Consequently, the corresponding estimation vector m
^
=
(
m
^
(
x
1
)
,
…
,
m
^
(
x
n
)
)
T generated by the penalized least square regression spline can be represented as
(7)
m
^
LS
=
X
β
^
LS
=
X
(
X
T
X
+
λ
D
)
-
1
X
T
y
.
For figuring out the penalty parameter λ in (7), the (generalized) cross-validation technique is employed hereafter, which is computed by leave-one-out of the residual from sum of squares so as to avoid overfitting in the regression spline. The cross-validation criterion is given by
(8)
CV
(
λ
)
=
1
n
∑
i
=
1
n
(
y
i
-
m
^
-
i
(
x
i
;
λ
)
)
2
,
where m
^
-
i
(
x
i
;
λ
) is the regression spline estimated by leaving out the ith observation point. The optimized penalty parameter λ can be calculated by minimizing CV
(
λ
) over λ
≥
0. For the linear smoothing matrix, H
λ
=
X
(
X
T
X
+
λ
D
)
-
1
X
T, and it can be verified [14] that
(9)
m
^
-
i
(
x
i
;
λ
)
=
∑
j
≠
i
(
H
λ
i
j
1
-
H
λ
i
j
y
i
)
.
Substitution of m
^
-
i
(
x
i
;
λ
) in (8) with the formula (9) results in
(10)
CV
(
λ
)
=
1
n
∑
i
=
1
n
(
y
i
-
m
^
(
x
i
;
λ
)
1
-
H
λ
i
i
)
2
.
The generalized cross-validation criterion is proposed [30] simply by replacing H
λ
i
i in (10) with their average tr
⁡
(
H
λ
)
/
n; namely,
(11)
GCV
(
λ
)
=
n
∑
i
=
1
n
(
y
i
-
m
^
(
x
i
;
λ
)
n
-
tr
⁡
(
H
λ
)
)
2
.


## 2.2.M-Type Estimator for Penalized Regression Splines

The above estimatem
^
LS of the penalized least squares regression splines may greatly suffer from various robustness issues especially when the data is contaminated with outliers. For alleviating this effect, a penalized regression estimator can be straightforwardly constructed by replacing the previous squared residual loss function with the following M-type criterion [27]:
(12)
∑
i
=
1
n
ρ
(
y
i
-
m
(
x
i
;
β
)
)
2
+
λ
∑
j
=
1
K
β
p
+
j
2
,
where ρ is even, nondecreasing in [
0
,
+
∞
) and ρ
(
0
)
=
0, for which a common choice is taken of the Huber loss function with cutoff c
>
0,
(13)
ρ
c
(
t
)
=
{
t
2
|
t
|
≤
c
2
c
|
t
|
-
c
2
|
t
|
>
c
.
Apparently, the above Huber’s function is a parabola in the vicinity of zero, and it increases linearly at a given level |
t
|
>
c. A default choice of the tuning constant c
=
1.345 aims for a 95% asymptotic efficiency with respect to the standard normal distribution.Denote the derivative ofρ
c as ψ
c
(
s
)
=
max
⁡
⁡
[
-
c
,
min
⁡
⁡
(
c
,
s
)
]. When a set of fitted residuals r
i
=
y
i
-
m
(
x
i
;
β
), i
=
1
,
…
,
n from the penalized spline regression are given, a robust M-scale estimator [19, 20] can be employed for estimating the standard deviation σ
ε of these residuals as follows:
(14)
∑
i
=
1
n
ψ
c
(
y
i
-
m
(
x
i
;
β
)
σ
^
ε
)
=
0
.
A common choice for a measure of scale with a high finite sample breakdown point can be the value of ω determined by
(15)
P
(
|
r
-
r
0.5
|
<
ω
)
=
1
2
,
where r
0.5 is the population median. Besides, the standard estimation of ω usually employs the median absolute deviation (MAD) statistic, which is defined as
(16)
MAD
=
Median
{
|
r
1
-
M
r
|
,
…
,
|
r
n
-
M
r
|
}
,
where M
r is the usual sample median of the fitted residuals and MAD is actually the sample median of the n values |
r
1
-
M
r
|
,
…
,
|
r
n
-
M
r
|, with its finite sample breakdown point approximately being 0.5 [20]. Suppose the fitted residuals are randomly sampled from a normal distribution; note that MAD does not estimate the corresponding standard deviation, σ
^
ε, but rather estimates Z
0.75
σ
^
ε (with Z
0.75 being the 0.75 quantile of the standard normal distribution). To put MAD in a more familiar context, it typically needs to be rescaled so as to estimate σ
^
ε especially when the residuals are sampled from a normal distribution. In particular,
(17)
MADN
=
MAD
z
0.75
≈
1.4836
MAD
.
By utilizing the aforementioned spline basis functions, an M-type penalized spline estimator using the standardized residuals can be computed as m
^
M
=
X
β
^
M, with its regression coefficients being as follows:
(18)
β
^
M
=
argmin
β
⁡
∑
i
=
1
n
ρ
c
(
y
i
-
m
(
x
i
;
β
)
σ
^
ε
)
+
λ
β
T
D
β
.
Due to the nonlinear nature of ρ
c and ψ
c, it is nontrivial to seek the minimizer of (18) while satisfying (14). Under such circumstances, an M-type iterative algorithm for calculating the penalized regression splines is proposed [27] by introducing the empirical pseudodata
(19)
z
=
m
^
p
+
ψ
c
(
y
-
m
^
p
)
2
,
where m
^
p is the penalized least square estimator corresponding to the above pseudodata
(20)
m
^
p
=
X
(
X
T
X
+
λ
D
)
-
1
X
T
z
.
As proved in Theorem 1 of Lee and Oh’s [27] work, the penalized least square estimator m
^
p converges to m
^
M asymptotically, which implies an M-type iterative algorithm consisting of the steps in Algorithm 1.Algorithm 1:
M-type iterative algorithm for penalized regression splines.
Step  0.  Input an initial curve estimate m
^
p
(
0
), termination tolerance ε, and the maximum iteration
number Itermax. Meanwhile, set k = 0 and conduct the following loop iterations.
Step  1.  Estimate σ
^
ε by using MADN in (17);
Step  2.  Produce the empirical pseudo data according to (19);
Step  3.  Calculate the penalized least-square estimator m
^
p
(
k
+
1
) defined in (20) for the pseudo data,
where the penalty parameter λ is chosen according to the generalized cross-validation
(GCV) criterion formulated by (11) in the foregoing subsection.
Step  4.  If ∥
m
^
p
(
k
+
1
)
-
m
^
p
(
k
)
∥
<
ε
∥
m
^
p
(
k
)
∥ or k = Itermax, terminate and output the robust
M-type penalized spline estimate m
^
M = m
^
p
(
k
+
1
); else set k = k + 1 and continue Step  1.Regarding theM-type estimator for the penalized regression splines (see Algorithm 1), the initial curve estimate m
^
p
(
0
) in Step 0 is suggested [27], choosing the nonrobust penalized least square regression spline presented by (7).


## 2.3.S-Estimation for Penalized Regression Splines

The robustness properties of the aforementionedM-type estimator for penalized regression splines are not completely satisfactory since they have low breakdown point [31], which was clarified by Maronna and Yohai [32] for the general M-estimators when considering both the regression and scale, simultaneously.Another more robust and relatively latest approach for smoothing the noisy data is theS-estimation for penalized regression splines [29], whose core idea is utilizing a flexible S-estimator to replace the traditional least squares estimator.  The coefficient vector of the penalized S-regression splines β
^
S is defined as follows:
(21)
β
^
S
=
argmin
β
⁡
{
n
σ
^
ε
2
(
β
)
+
λ
β
T
D
β
}
,
in which, for each β, σ
^
ε
2
(
β
) satisfies
(22)
1
n
∑
i
=
1
n
ρ
(
y
i
-
F
(
x
i
)
T
β
σ
^
ε
(
β
)
)
=
E
Φ
(
ρ
(
Z
)
)
,
where Φ is cumulative distribution function of the standard norm distribution, Z
~
N
(
0,1
), and the loss function is taken from Tukey’s bisquare family [33]
(23)
ρ
d
(
u
)
=
{
3
(
u
d
)
2
-
3
(
u
d
)
4
+
(
u
d
)
6
|
u
|
≤
d
1
|
u
|
>
d
.
It is evidently known from theoretical assumption that each fitted residual vector should follow a normal distribution with mean of zero; that is, r
i
~
N
(
0
,
σ
ε
2
), i
=
1
,
…
,
n, however its median statistic, Median(
r
), in which r
=
(
r
1
,
…
,
r
n
)
T, may greatly depart from zero to a large extent in the general practice. Considering this, we ameliorate the estimation of σ
^
ε
(
β
) by using MADN in (17) with the addition of the absolute value of residual median; namely,
(24)
σ
^
ε
(
β
)
=
MADN
+
Median
(
r
(
β
)
)
.
The penalized S-estimator for the regression spline model is computed as m
^
S
=
X
β
^
S, with its coefficients vector β
^
S satisfying both (21) and (22), simultaneously. The following formula derivations indicate that such penalized S-estimator can actually be equivalent to a weighted penalized least squares regression.Before proceeding with the derivation, it is forewarned that the symbolβ
^
S is simply written as β just for clarity. Firstly, when (21) reaches its minimum, the following equation should be satisfied:
(25)
n
σ
^
ε
(
β
)
∇
σ
^
ε
(
β
)
+
λ
D
β
=
0
.
On the other hand, taking the derivative of the M-scale function in (22) with respect to β leads to
(26)
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
σ
^
ε
(
β
)
)
(
-
F
(
x
i
)
T
σ
^
ε
(
β
)
-
r
i
(
β
)
∇
σ
^
ε
(
β
)
T
σ
^
ε
2
(
β
)
)
=
0
,
in which r
i
(
β
)
=
y
i
-
F
(
x
i
)
T
β, r
~
i
(
β
)
=
r
i
(
β
)
/
σ
^
ε
(
β
), and Λ
(
β
)
=
diag
⁡
(
ρ
′
(
r
~
i
(
β
)
)
/
r
~
i
(
β
)
),  i
=
1
,
…
,
n.Deformation of the above formula produces(27)
∇
σ
^
ε
(
β
)
=
-
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
F
(
x
i
)
∑
i
=
1
n
ρ
′
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
(
r
i
(
β
)
/
σ
^
ε
(
β
)
)
=
-
(
σ
^
ε
(
β
)
X
T
Λ
(
β
)
r
(
β
)
r
(
β
)
T
Λ
(
β
)
r
(
β
)
)
.
Substituting ∇
σ
^
ε
(
β
) of (27) into (25),
(28)
-
n
σ
^
ε
(
β
)
2
r
(
β
)
T
Λ
(
β
)
r
(
β
)
X
T
Λ
(
β
)
r
(
β
)
+
λ
D
β
=
0
.
Recall that r
(
β
)
=
y
-
X
β and import new symbol τ
(
β
)
=
n
σ
^
ε
(
β
)
2
/
r
(
β
)
T
Λ
(
β
)
r
(
β
); we have
(29)
-
τ
(
β
)
X
T
Λ
(
β
)
(
y
-
X
β
)
+
λ
D
β
=
0
.
And, thus, a further simple deformation arrives at the coefficient vector β
^
S with an iteration form
(30)
β
=
(
X
T
Λ
(
β
)
X
+
λ
τ
(
β
)
D
)
-
1
X
T
Λ
(
β
)
y
,
which forms the basis of constructing an iterative algorithm for penalized S-regression spline.The above formula derivations verify thatβ
^
S is actually the solution to the following weighted penalized least squares regression problem:
(31)
min
⁡
β
⁡
∥
Λ
(
β
)
1
/
2
(
y
-
X
β
)
∥
2
+
λ
τ
(
β
)
β
T
D
β
.
And such representation indicates that the GCV criterion in (11) can be employed straightforwardly with response variable y
~
=
Λ
(
β
)
1
/
2
y, weighted design matrix X
~
=
Λ
(
β
)
1
/
2
X, and the penalty term λ
/
τ
(
β
). Meanwhile, due to some weight elements being zero, the penalty parameter λ is thus calculated based on the following Regularized GCV criterion:
(32)
RGCV
(
λ
)
=
n
w
∥
Λ
(
β
)
1
/
2
(
y
-
X
β
)
∥
2
(
n
w
-
tr
⁡
(
H
~
λ
)
)
2
,
where H
~
λ
=
X
~
(
X
~
T
X
~
+
(
λ
/
τ
(
β
)
)
D
)
-
1
X
~
T
=
Λ
(
β
)
1
/
2
X
(
X
~
T
Λ
(
β
)
X
+
(
λ
/
τ
(
β
)
)
D
)
-
1
X
T
Λ
(
β
)
1
/
2 and n
w denotes the number of nonzero elements in weight matrix Λ
(
β
).Consequently, based on the foregoing in-depth discussions, theS-estimator iterative algorithm for penalized regression splines can be constructed as in Algorithm 2.Algorithm 2:
S-estimator iterative algorithm for penalized regression splines.
Step  0. Input an initial regression coefficients β
^
(
0
), the termination tolerance ε, and the maximum
iteration number Itermax. Meanwhile, set k = 0 and enter into the following loop iterations.
Step  1. Estimate σ
^
ε based on (24), that is
σ
^
ε
(
β
^
(
k
)
)
=
MADN
+
Median
(
r
(
β
^
(
k
)
)
).
Step  2. Compute the intermediate variables τ
(
β
^
(
k
)
) and Λ
(
β
^
(
k
)
).
Step  3. Calculate the coefficient vector β
^
(
k
+
1
) according to (30) as
β
^
(
k
+
1
)
=
(
X
T
Λ
(
β
^
(
k
)
)
X
+
λ
τ
(
β
^
(
k
)
)
D
)
-
1
X
T
Λ
(
β
^
(
k
)
)
y.
Hereinto, the penalty parameter λ is figured out based on the regularized GCV criterion (32).
Step  4. If ∥
β
^
(
k
+
1
)
-
β
^
(
k
)
∥
<
ε
∥
β
^
(
k
)
∥ or k = Itermax, terminate and output the convergent
coefficient vector for the penalized regression spline of S-estimator β
^
★
=
β
^
(
k
+
1
); else
accumulate the counter by setting k = k + 1 and repeat Step  1.In particular, it should be pointed out that since the employed loss functionρ
d is bounded giving rise to the nonconvexity of σ
^
ε
(
β
), the objective function in (21) may have several critical points, each of which only corresponds to one local minima. Therefore, the convergent critical point of the iterations derived from the above algorithm completely depends on the prescribed initial regression coefficients vector, β
^
(
0
), in Step 0 of Algorithm 2.One possible solution to cope with this issue proposed by Tharmaratnam et al. [29] is based on the random sampling theory, which can be briefly described as, firstly, creating J random subsamples with a proper size, for example, one fifth of the sample number; namely, n
s
=
max
⁡
⁡
(
p
+
K
+
1
,
⌊
n
/
5
⌋
), from the original sample observations. Each subsample brings about an initial candidate of the coefficient vector, for instance, assigned by least squares fitting for the penalized regression splines or the M-type estimator with a small amount of additional computation cost, with the initial regression coefficient vector of the jth subsample being denoted by β
^
j
(
0
), staring from which the final iterative convergent vector produced through Algorithm 2 is being written as β
^
j
⋆, j
=
1
,
…
,
J. Eventually, the global optimum is taken among all these J convergent vectors to be the one corresponding to the minimum objective function in (21); namely,
(33)
β
^
S
⋆
=
argmin
1
≤
j
≤
J
⁡
{
n
σ
^
ε
2
(
β
^
j
⋆
)
+
λ
β
^
j
⋆
T
D
β
^
j
⋆
}
.
Therefore, the robust S-estimator for penalized regression splines can be computed as m
^
S
=
X
β
^
S
⋆.


## 3. Experimental Performance Evaluation

This section is devoted to the comparative experiments upon both simulated synthetic data and one real weather balloon data set for the foregoing explored penalized spline smoothing methods, including the nonrobust penalized least square regression spline (LS) presented by (7) and two robust smoothing approaches, the M-type estimator [27] and S-estimation [29], for penalized regression splines as described in Algorithms 1 and 2.


### 3.1. Configuration and Criteria

It should be incidentally mentioned that the exponential orderp
=
3 for truncated power basis is always employed for all these three methods, and the knot locations are selected in accordance with the recommended choice by Ruppert et al. [14], κ
k
=
(
(
k
+
1
)
/
(
K
+
2
)
)th sample quantile of the Unique(x),  k
=
1
,
…
,
K; together with a simple choice of the knot number K
=
min
⁡
⁡
(
35
,
(
1
/
4
)
×
Number
of
Unique
(
x
)
), where x
=
(
x
1
,
x
2
,
…
,
x
n
)
T. Besides, the corresponding penalty parameter λ is determined according to the generalized cross-validation (GCV) criterion for the LS fitting and M-estimation methods and the regularized GCV criterion for the S-estimation method. Meanwhile, the termination tolerance ε
=
10
-
6 and the maximum iteration number of 100 are assigned for the initial input parameters in Algorithms 1 and 2 and J
=
5 subsamples are randomly generated for exploring the ultimate optimum.Meanwhile, the hardware equipment and software environment established for the experiments below are Dell XPS 8700 Desktop (CPU: Intel Core i7-4770 @ 3.40 GHz), MATAB R2014a installed upon Windows 7 SP1 Ultimate (64 Bit).


### 3.2. Synthetic Data

The simulated study is designed the same as in Lee and Oh [27] with the testing function
(34)
y
i
=
sin
⁡
(
2
π
(
1
-
x
i
)
2
)
+
0.5
ε
i
,
i
=
1
,
…
,
n
.
The designed samples are generated with their horizontal axes being uniformly extracted from the open interval (
0,1
) and the corresponding error distribution ε employing nine possible scenarios with the probability distribution function (pdf) specified if necessary, including (1) uniform distribution, U
(
-
1,1
), (2) standard normal distribution, N
(
0,1
), (3) logistic distribution with its pdf as e
x
/
(
1
+
e
x
)
2, (4) double exponential distribution with its pdf as (
(
1
/
2
)
e
-
|
x
|
), (5) contaminated normal distribution I, 0.95
N
(
0,1
)
+
0.05
N
(
0,10
), that is the error constituted with 95% from N
(
0,1
) and 5% from N
(
0,10
),  (6) contaminated normal distribution II, 0.9
N
(
0,1
)
+
0.1
N
(
0,10
), (7) slash distribution, defined from two statistically independent random variables following (1) and (2), that is, N
(
0,1
)
/
U
(
0,1
), (8) Cauchy distribution with the pdf as 1
/
π
(
1
+
x
2
), whose simulated outliers can be generated from uniform distribution by inverse transformation technique [34], and (9) an asymmetric normal mixture distribution, 0.9
N
(
0,1
)
+
0.1
N
(
30,1
). For the sake of reducing the simulation variability, the sample size is taken as a fixed value, n
=
200. It should be mentioned incidentally that the above (1)~(8) error scenarios are symmetrically distributed and provided with the ascending order according to the heaviness of their tails, and scenario (9) denotes the standard normal distribution mixed with 10% asymmetric outliers from N
(
30,1
).Figure1 presents the spline regression curves from the penalized LS fitting, M-estimator, and S-estimator for one time simulated scattered samples in nine possible scenarios, each of which with a subfigure panel demonstrating the comparative regression results for these three penalized spline smoothing methods.Spline regression curves (with true testing function in black solid line) from the penalized LS fitting (blue dotted),M-estimator (green dot-dashed), and S-estimator (red dashed) for the simulated samples (black circle-dot) in 9 possible scenarios.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)For proceeding further quantitative evaluation of the comparative analysis for these three penalized spline smoothing methods, the statistical indicator of average squared error (ASE) is utilized here to qualify the goodness of fitting between the regression splinem
^
j
(
x
) for the jth simulation replicate from each of the above outlier scenarios and the true simulated function m
(
x
)
=
sin
⁡
(
2
π
(
1
-
x
)
2
)
(35)
ASE
j
=
1
n
∑
i
=
1
n
(
m
(
x
i
)
-
m
^
j
(
x
i
)
)
2
,
j
=
1
,
…
,
J
,
with the specified replicate number J
=
10
3 for each distribution.Performances of the penalized spline smoothing methods are thoroughly evaluated from the aspects of the fitting accuracy, robustness, and execution time. In the following, Figures2 and 3 illustrate the box plots of both ASEs and code running times using (a) penalized LS-estimation, (b) penalized M-estimation, and (c) penalized S-estimation for the spline smoothing methods, as well as their median value comparisons of these two indicators for the simulated samples from each of the 9 distribution scenarios.Box plots of ASEs using (a) penalized LS-estimation, (b) penalizedM-estimation, (c) penalized S-estimation, and (d) the comparison of Median ASE for the simulated samples from each of the 9 distribution scenarios with replicate number J
=
10
3.
(a)
(b)
(c)
(d)Box plots of run times using (a) penalized LS-estimation, (b) penalizedM-estimation, (c) penalized S-estimation, and (d) the comparison of median run time for the simulated samples from each of the 9 distribution scenarios with replicate number J
=
10
3.
(a)
(b)
(c)
(d)From the ASE comparisons in Figure2, it can be found that the penalized LS-estimation behaves well for the symmetric error scenarios with small tail weights, such as the simple uniform and standard normal distribution. However, it cannot cope with the complex conditions such as the three latter error distributions possessing heavy tails, thereby binging about large ASEs which are specially rescaled to the additive vertical coordinate on the right so as to be accommodated with the ASEs results for other error sources in the same window (refer to panels (a) and (d) in Figure 2). In comparison, the robust penalized M-type estimation presents satisfactory regression fittings for the distributions with moderate tail weights like (2)~(5) scenarios, which are located at the middle in each panel. Meanwhile, the robust S-estimator provides the comparatively meritorious smoothing results especially for the complicated circumstances like slash, Cauchy, or the asymmetric distributions with heavy tails.In terms of the robustness, LS-estimator simply does not have the qualifications to be mentioned from this aspect. Instead, theM-estimator demonstrates stable performance for the designed samples with moderate perturbation error and the S-estimator behaves dominantly stable for all kinds of error source scenarios.Unfortunately, as for the execution efficiency, the most time-elapsed one is theS-estimator, followed by the M-estimator, with the least consuming estimator being the LS-estimator since it has been treated as the initial curve estimate for the two former approaches. Moreover, although both S-estimator and M-estimator are composed of similar iterative process, the random sampling theory is further involved in the former model for exploring the ultimate solution for a nonconvex optimization problem.


### 3.3. Balloon Data

This section is devoted to the comparative experiments of the aforementioned penalized spline smoothing methods upon a set of weather balloon data [35]. The data consist of 4984 observations collected from a balloon at about 30 kilometers altitude, with the outliers caused due to the fact that the balloon slowly rotates, causing the ropes from which the measuring instrument is suspended to cut off the direct radiation from the sun.As illustrated above, Figure4 displays these balloon observational data together with the regression fitting curves by the penalized LS method, the robust penalized M-estimation method, and the robust S-estimation method. Besides, it can also be explicitly discovered by contrast that the non-robust LS regression curve suffers from the disturbance of the outliers, and this phenomenon seems even more significantly around the position of horizontal axis being 0.8, while the other two robust splines pass through approximately the majority of the observations. Meanwhile, the S-estimator seems to behave slightly more superior to the M-estimator since the former accords with the data trend better especially at the ridge and two terminal nodes.Figure 4
Regression fitted curves for the balloon data using penalized spline smoothing models including the penalized LS fitting (blue dashed),M-estimator (green dot-dashed), and S-estimator (red solid line).At the same time, under the platform specified in Section3.1, the algorithm execution times of this weather balloon data set with 4984 observations are 2.3 s, 35.1 s, and 58.9 s for the penalized LS regression, the penalized M-type estimator, and the penalized S-estimator, respectively.


## 3.1. Configuration and Criteria

It should be incidentally mentioned that the exponential orderp
=
3 for truncated power basis is always employed for all these three methods, and the knot locations are selected in accordance with the recommended choice by Ruppert et al. [14], κ
k
=
(
(
k
+
1
)
/
(
K
+
2
)
)th sample quantile of the Unique(x),  k
=
1
,
…
,
K; together with a simple choice of the knot number K
=
min
⁡
⁡
(
35
,
(
1
/
4
)
×
Number
of
Unique
(
x
)
), where x
=
(
x
1
,
x
2
,
…
,
x
n
)
T. Besides, the corresponding penalty parameter λ is determined according to the generalized cross-validation (GCV) criterion for the LS fitting and M-estimation methods and the regularized GCV criterion for the S-estimation method. Meanwhile, the termination tolerance ε
=
10
-
6 and the maximum iteration number of 100 are assigned for the initial input parameters in Algorithms 1 and 2 and J
=
5 subsamples are randomly generated for exploring the ultimate optimum.Meanwhile, the hardware equipment and software environment established for the experiments below are Dell XPS 8700 Desktop (CPU: Intel Core i7-4770 @ 3.40 GHz), MATAB R2014a installed upon Windows 7 SP1 Ultimate (64 Bit).


## 3.2. Synthetic Data

The simulated study is designed the same as in Lee and Oh [27] with the testing function
(34)
y
i
=
sin
⁡
(
2
π
(
1
-
x
i
)
2
)
+
0.5
ε
i
,
i
=
1
,
…
,
n
.
The designed samples are generated with their horizontal axes being uniformly extracted from the open interval (
0,1
) and the corresponding error distribution ε employing nine possible scenarios with the probability distribution function (pdf) specified if necessary, including (1) uniform distribution, U
(
-
1,1
), (2) standard normal distribution, N
(
0,1
), (3) logistic distribution with its pdf as e
x
/
(
1
+
e
x
)
2, (4) double exponential distribution with its pdf as (
(
1
/
2
)
e
-
|
x
|
), (5) contaminated normal distribution I, 0.95
N
(
0,1
)
+
0.05
N
(
0,10
), that is the error constituted with 95% from N
(
0,1
) and 5% from N
(
0,10
),  (6) contaminated normal distribution II, 0.9
N
(
0,1
)
+
0.1
N
(
0,10
), (7) slash distribution, defined from two statistically independent random variables following (1) and (2), that is, N
(
0,1
)
/
U
(
0,1
), (8) Cauchy distribution with the pdf as 1
/
π
(
1
+
x
2
), whose simulated outliers can be generated from uniform distribution by inverse transformation technique [34], and (9) an asymmetric normal mixture distribution, 0.9
N
(
0,1
)
+
0.1
N
(
30,1
). For the sake of reducing the simulation variability, the sample size is taken as a fixed value, n
=
200. It should be mentioned incidentally that the above (1)~(8) error scenarios are symmetrically distributed and provided with the ascending order according to the heaviness of their tails, and scenario (9) denotes the standard normal distribution mixed with 10% asymmetric outliers from N
(
30,1
).Figure1 presents the spline regression curves from the penalized LS fitting, M-estimator, and S-estimator for one time simulated scattered samples in nine possible scenarios, each of which with a subfigure panel demonstrating the comparative regression results for these three penalized spline smoothing methods.Spline regression curves (with true testing function in black solid line) from the penalized LS fitting (blue dotted),M-estimator (green dot-dashed), and S-estimator (red dashed) for the simulated samples (black circle-dot) in 9 possible scenarios.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)For proceeding further quantitative evaluation of the comparative analysis for these three penalized spline smoothing methods, the statistical indicator of average squared error (ASE) is utilized here to qualify the goodness of fitting between the regression splinem
^
j
(
x
) for the jth simulation replicate from each of the above outlier scenarios and the true simulated function m
(
x
)
=
sin
⁡
(
2
π
(
1
-
x
)
2
)
(35)
ASE
j
=
1
n
∑
i
=
1
n
(
m
(
x
i
)
-
m
^
j
(
x
i
)
)
2
,
j
=
1
,
…
,
J
,
with the specified replicate number J
=
10
3 for each distribution.Performances of the penalized spline smoothing methods are thoroughly evaluated from the aspects of the fitting accuracy, robustness, and execution time. In the following, Figures2 and 3 illustrate the box plots of both ASEs and code running times using (a) penalized LS-estimation, (b) penalized M-estimation, and (c) penalized S-estimation for the spline smoothing methods, as well as their median value comparisons of these two indicators for the simulated samples from each of the 9 distribution scenarios.Box plots of ASEs using (a) penalized LS-estimation, (b) penalizedM-estimation, (c) penalized S-estimation, and (d) the comparison of Median ASE for the simulated samples from each of the 9 distribution scenarios with replicate number J
=
10
3.
(a)
(b)
(c)
(d)Box plots of run times using (a) penalized LS-estimation, (b) penalizedM-estimation, (c) penalized S-estimation, and (d) the comparison of median run time for the simulated samples from each of the 9 distribution scenarios with replicate number J
=
10
3.
(a)
(b)
(c)
(d)From the ASE comparisons in Figure2, it can be found that the penalized LS-estimation behaves well for the symmetric error scenarios with small tail weights, such as the simple uniform and standard normal distribution. However, it cannot cope with the complex conditions such as the three latter error distributions possessing heavy tails, thereby binging about large ASEs which are specially rescaled to the additive vertical coordinate on the right so as to be accommodated with the ASEs results for other error sources in the same window (refer to panels (a) and (d) in Figure 2). In comparison, the robust penalized M-type estimation presents satisfactory regression fittings for the distributions with moderate tail weights like (2)~(5) scenarios, which are located at the middle in each panel. Meanwhile, the robust S-estimator provides the comparatively meritorious smoothing results especially for the complicated circumstances like slash, Cauchy, or the asymmetric distributions with heavy tails.In terms of the robustness, LS-estimator simply does not have the qualifications to be mentioned from this aspect. Instead, theM-estimator demonstrates stable performance for the designed samples with moderate perturbation error and the S-estimator behaves dominantly stable for all kinds of error source scenarios.Unfortunately, as for the execution efficiency, the most time-elapsed one is theS-estimator, followed by the M-estimator, with the least consuming estimator being the LS-estimator since it has been treated as the initial curve estimate for the two former approaches. Moreover, although both S-estimator and M-estimator are composed of similar iterative process, the random sampling theory is further involved in the former model for exploring the ultimate solution for a nonconvex optimization problem.


## 3.3. Balloon Data

This section is devoted to the comparative experiments of the aforementioned penalized spline smoothing methods upon a set of weather balloon data [35]. The data consist of 4984 observations collected from a balloon at about 30 kilometers altitude, with the outliers caused due to the fact that the balloon slowly rotates, causing the ropes from which the measuring instrument is suspended to cut off the direct radiation from the sun.As illustrated above, Figure4 displays these balloon observational data together with the regression fitting curves by the penalized LS method, the robust penalized M-estimation method, and the robust S-estimation method. Besides, it can also be explicitly discovered by contrast that the non-robust LS regression curve suffers from the disturbance of the outliers, and this phenomenon seems even more significantly around the position of horizontal axis being 0.8, while the other two robust splines pass through approximately the majority of the observations. Meanwhile, the S-estimator seems to behave slightly more superior to the M-estimator since the former accords with the data trend better especially at the ridge and two terminal nodes.Figure 4
Regression fitted curves for the balloon data using penalized spline smoothing models including the penalized LS fitting (blue dashed),M-estimator (green dot-dashed), and S-estimator (red solid line).At the same time, under the platform specified in Section3.1, the algorithm execution times of this weather balloon data set with 4984 observations are 2.3 s, 35.1 s, and 58.9 s for the penalized LS regression, the penalized M-type estimator, and the penalized S-estimator, respectively.


## 4. Conclusions

This paper conducts a comprehensively comparative analysis of the current two popular robust penalized spline smoothing methods, the penalizedM-type estimator and the penalized S-estimation, both of which are reelaborated starting from their origins, with their derivation process reformulated and the ultimate algorithms reorganized under a unified framework.Experiments consist of firstly one simulated synthetic data set with the outliers in 9 possible scenarios and then a practical weather balloon data set for the aforementioned penalized regression spline, whose performances are thoroughly evaluated from the aspects of fitting accuracy, robustness, and execution time.Conclusions from these experiments are that robust penalized spline smoothing methods can be resistant to the influence of outliers compared with the nonrobust penalized LS spline fitting method. Further, theM-estimator can exert stable performance only for the observed samples with moderate perturbation error while the S-estimator behaves dominantly well even for the cases with heavy tail weights or higher percentage of contaminations, but necessary with more executive time. Therefore, taking both the solving accuracy and implementation efficiency of the fitted regression spline into account, M-type estimator is generally preferred, whereas if only the fitting credibility is of concern, then the S-estimator naturally deserves to become the first choice.



---

*Source: 642475-2014-07-16.xml*