Mem. Gra. Sci. Eng. Shimane Univ.
Series B: Mathematics
47 (2014), pp. 63–71

ASYMPTOTICS FOR PENALIZED SPLINES IN ADDITIVE
MODELSTAKUMA YOSHIDA
Communicated by Kanta Naito(Received: December 19, 2013)
Abstract. This paper is based on the author’s thesis, “Asymptotic theory ofpenalized spline regression”. The focus of the present paper is on the penalized
spline estimators obtained by the backfitting algorithm in additive models. The
convergence of the algorithm as well as the uniqueness of its solution are shown.
Asymptotic equivalence between the penalized spline estimators by the backfitting algorithm and the convenient estimators proposed by Marx and Eilers [9] isaddressed. Asymptotic normality of the estimators is also developed.

1. Introduction
The additive model is a typical regression model with multidimensional covariates and is usually expressed asyi = f1 (xi1 ) + · · · + fD (xiD ) + εi ,for given data {(yi , xi1 , · · · , xiD ) : i = 1, · · · , n}, where each fd (d = 1, · · · , D) is a
univariate function with a certain degree of smoothness.
The additive model has become a popular smoothing technique and its fundamental properties have been summarized in literature such as Buja et al. [2] andHastie and Tibshirani [7]. Buja et al. [2] proposed the so-called backfitting algorithm, which is eﬃcient for nonparametric estimation of fd (d = 1, · · · , D). The
backfitting algorithm is a repetition update algorithm and its convergence and theuniqueness of its solution are not always assured. Buja et al. [2] showed the suﬃcient condition for convergence of the backfitting algorithm and the uniqueness of
its solution for some smoothing methods.
In this paper, we discuss the asymptotic properties of the penalized spline estimator for the additive model. Penalized spline estimators have been discussed
in O’Sullivan [10], Eilers and Marx [5], Marx and Eilers [9], Aerts et al. [1] and2000 Mathematics Subject Classification. 62G08, 62G20, 65D07.
Key words and phrases. Additive model, Asymptotic normality, Backfitting algorithm, Penalized spline.
63

64

T. YOSHIDA

Ruppert et al. [11]. Despite its richness of application, the asymptotics for splinesmoothing seem to have not yet been suﬃciently developed.
For the univariate model (D = 1), Hall and Opsomer [6] gave the mean squared
error and proved consistency of the penalized spline estimator. The asymptotic
bias and variance of the penalized spline estimator were obtained in Claeskens etal. [3]. Kauermann et al. [8] worked with the generalized linear model. Wang et al.
[14] showed that the penalized spline estimator is asymptotically equivalent to a
Nadaraya-Watson estimator. Thus, it seems that developments of the asymptotic
theories of the penalized splines are relatively recent events and we note that thoseworks are mainly regarding the univariate model.
The penalized spline estimators for the additive models are obtained by using
the penalized least squares method. However there is one problem in that the loss
function L is not strictly convex. Therefore it is diﬃcult to find the minimizerof L since its Hessian is not invertible. Though the backfitting algorithm yields
a solution that makes the gradient of L equal to zero, there is no guarantee thatthe solution obtained by this method minimizes L. It has been known that the
backfitting algorithm does converge, but the uniqueness of the obtained solution
cannot be proved in general.
On the other hand, Marx and Eilers [9] proposed a new loss function to avoid thissingularity problem. They proposed to use Lγ which consists of L plus an additional
small ridge penalty. This Lγ is strictly convex, hence the global minimum of Lγ
is equivalent to a unique local minimum which can be easily obtained. Of course
there is a gap between the minimizer of Lγ and that of L. In this paper, theestimator obtained by minimizing Lγ is called the ridge corrected penalized spline
estimator (RCPS) and the penalized spline estimator obtained by the backfitting
algorithm is denoted as the backfitting penalized spline estimator (BPS). Because
we are interested in the estimator obtained by using L, we mainly focus on BPS.
The aim of this paper is to derive the asymptotic distribution of the BPS in
the general D-variate additive model. First we show the asymptotic distributionof the RCPS. Next, it is shown that the diﬀerence of the RCPS and the BPS
asymptotically vanishes. As a result, it demonstrates that the asymptotic normalityof the BPS and that the solution of the backfitting algorithm is asymptoticallyunique. As will be seen in the subsequent section, although the closed form of the
BPS can not be written, its asymptotic properties can be shown in each iterationof the algorithm. The properties of the band matrices play an important role as amathematical tool in asymptotic considerations.
This paper is organized as follows. In Section 2, our model settings and estimating equation in the penalized least squares method are discussed, and the RCPS
and the BPS are constructed. Section 3 provides the asymptotic bias and varianceof the RCPS, and then its asymptotic normality is developed. Furthermore, weshow that the BPS is asymptotically equivalent to the RCPS and the solution of
the backfitting algorithm is asymptotically unique. In Section 4, we give somecomments. Proofs of all mathematical results are omitted.

ASYMPTOTICS FOR PENALIZED SPLINES IN ADDITIVE MODELS

65

2. Model setting and penalized spline estimator
2.1. Additive spline model. Consider a D-variate additive regression model(1)

yi = f1 (xi1 ) + · · · + fD (xiD ) + εi

for the data {(yi , xi ) : i = 1, · · · , n}, where xi = (xi1 , · · · , xiD ) is the D-variate
explanatory variables, fj (·) is an unknown regression function and εi ’s are independent random errors with E[εi |X i = xi ] = 0 and V [εi |X i = xi ] = σ 2 (xi ) < ∞. We
assume E[fj (Xj )] = 0(j = 1, · · · , D) to ensure identifiability of fj . Let qj (xj ) be thedensity of Xj and q(x) = q(x1 , · · · , xD ) be the joint density of X = (X1 , · · · , XD ).
We assume without loss of generality that xi ∈ (0, 1)D for all i ∈ {1, · · · , n}, where(0, 1)D is the D-variate unit cube.
Now we consider the B-spline modelKn∑sj (xj ) =[p]

Bk (xj )bj,k

k=−p+1[p]

as an approximation to fj (xj ) at any xj ∈ (0, 1) for j = 1, · · · , D. Here, Bk (x)(k =−p + 1, · · · , Kn ) are the pth degree B-spline basis functions defined recursively as{
1, κk−1 < x ≤ κk ,[0]
Bk (x) =0,
otherwise,x − κk−1
κk+p − x [p−1]
[p]
[p−1]
Bk (x) =Bk (x) +B(x),κk+p−1 − κk−1
κk+p − κk k+1
where κk = k/Kn (k = −p + 1, · · · , Kn + p) are knots and bj,k (j = 1, · · · , D, k =[p]−p + 1, · · · , Kn ) are unknown parameters. We denote Bk (x) as Bk (x) in what
follows since only the pth degree is treated. The details and many properties of
the B-spline function are clarified in de Boor [4]. We aim to obtain an estimatorof fj via the B-spline additive regression model(2)

yi = s1 (xi1 ) + · · · + sD (xiD ) + εiinstead of model (1). Model (2) can be expressed as
y = Z1 b1 + · · · + ZD bD + ε = Zb + ε
by using the notations y = (y1 · · · yn )′ , bd = (bd,−p+1 · · · bd,Kn )′ , b = (b′1 · · · b′D )′ ,Zd = (B−p+j (xid ))ij , Z = [Z1 · · · ZD ] and ε = (ε1 · · · εn )′ . The estimator b̂ of b
is defined as the minimizer ofL(b) = (y − Zb)′ (y − Zb) +
(3)

′

D
∑

λjn b′j Qm bj

j=1′= (y − Zb) (y − Zb) + b Qm (λn )b,

where λjn (j = 1, · · · , D) are the smoothing parameters, Qm is the mth order
diﬀerence matrix and Q(λn ) = diag[λ1n Qm · · · λDn Qm ]. This estimation method

66

T. YOSHIDA

is called the penalized least squares method and it has been frequently utilized in
spline regression. For a fixed point xj ∈ (0, 1), the estimator fˆj (xj ) of fj (xj ) isKn∑fˆj (xj ) =Bk (xj )b̂j,kk=−p+1and is called the penalized spline estimator of fj (xj ). The predictor of y at a fixed
point x ∈ (0, 1)D is defined as(4)ŷ = fˆ1 (x1 ) + · · · + fˆD (xD ).

Since E[fj (Xj )] = 0 is assumed for fj , the estimator of each component fj isusually centered. Hence, fˆj (xj ) is rewritten as
n
1∑ˆ
ˆ
ˆ
fj,c (xj ) = fj (xj ) −fj (xij ),n i=1

as discussed in Wang and Yang [13]. In this paper, however, we do not examine fˆj,cbecause our interests are in asymptotics for fˆj and ŷ, and asymptotic distributionsof fˆj (xj ) and fˆj,c (xj ) become equivalent.
2.2. The Ridge Corrected Penalized Spline Estimator. In general, b̂ =′′(b̂1 · · · b̂D )′ is a solution of
∂L(b)
= 0.
∂b
However, this method has one defect: the L(b) is not strictly convex as a function of
b in general. Hence, the solution of (5) does not necessarily become the minimizerof (3). Actually, because each column sum of Zj equal to 1, Z ′ Z + Qm (λn ), theHessian matrix of L(b), has eigenvalue 0. Marx and Eilers [9] also noted this point
as a typical problem of the additive spline regression. They studied a new method
such that the loss function has strict convexity for obtaining the estimator of b.
Let Lγ (b) = L(b) + γb′ b, where γ > 0 is very small constant. Since Lγ (b) is strictly′′
convex, the solution b̂γ = (b̂1,γ · · · b̂D,γ )′ of(5)∂Lγ (b)=0∂b
can be obtained uniquely asb̂γ = (Z ′ Z + Q(λn ) + γI)−1 Z ′ y.
The RCPS of fj (xj ) can be written asfˆj,γ (xj ) = B(xj )′ b̂j,γ ,where B(xj ) = (B−p+1 (xj ) · · · BKn (xj ))′ . If γ = 0, b̂γ cannot be calculated
because Z ′ Z + Q(λn ) is not invertible.

ASYMPTOTICS FOR PENALIZED SPLINES IN ADDITIVE MODELS

67

2.3. The Backfitting Penalized Spline Estimator. The merit and usage ofthe backfitting algorithm are clarified in Hastie and Tibshirani [7]. The ℓ-stage(ℓ)
backfitting estimator bj of bj is defined as′
bj = Λ−1
j Zj (y − Z1 b1 − · · · − Zj−1 bj−1 − Zj+1 bj+1 − · · · − ZD bD
(ℓ)

(ℓ)

(ℓ)

(ℓ−1)

(ℓ−1)

(0)

),

(ℓ)

where bj ’s are initial values. The ℓ-stage backfitting estimator fj (xj ) of fj (xj )at xj ∈ (0, 1) is defined as
Kn
∑

(ℓ)
fj (xj ) =Bk (xj )bj,k = B(xj )′ bj , j = 1, · · · , D.
(ℓ)

(ℓ)

k=−p+1(ℓ)

For D = 2, the explicit form of bj can be obtained (see Yoshida and Naito(ℓ)[15]). However, for general D, the exact form of bj is too complicated to be
written down. A mathematical property of the backfitting algorithm is that b(∞) =
(∞)
(∞)
(ℓ)
(ℓ)
((b1 )′ , · · · , (bD )′ )′ ≡ limℓ→∞ ((b1 )′ , · · · , (bD )′ )′ satisfies
∂L(b)
= 0.
∂b b=b(∞)(ℓ)(∞)

It is shown by Theorem 9 of Buja et al. [2] that bj converges to bj , but those(∞)(0)
bj (j = 1, ..., D) are depending on initial values bj (j = 1, ..., D). This meansthat the convergence property of the backfitting estimator is guaranteed, but theuniqueness of the solutions is not trivial. We will study the asymptotic behavior of
(∞)
(∞)
(∞)
fˆj (xj ) = fj (xj ) = B(xj )′ bj , as well as the relationship between b̂j,γ and bj .
3. Asymptotic theory
We prepare some symbols and notations to be used hereafter. Define the (Kn +p) × (Kn + p) square matrix Gk = (Gk,ij )ij with its (i, j)-component∫ 1
Gk,ij =
B−p+i (x)B−p+j (x)qk (x)dx0

for k = 1, · · · , D and the (Kn + p) × (Kn + p) square matrix Σk = (Σk,ij )ij having
the (i, j)-component∫
Σk,ij =σ 2 (x)B−p+i (xk )B−p+j (xk )q(x)dx[0,1]Dfor k = 1, · · · , D. Let b∗j be a best L∞ approximation to the true function fj . Thismeans that b∗j satisfies
sup fj (x) + bj,a (x) − B(x)′ b∗j = o(Kn−(p+1) ),x∈(0,1)

where
(p+1)

bj,a (x) = −

fj

(x)

Kn
∑

Knp+1 (p + 1)! k=1(
I(κk−1 ≤ x < κk )Bp+1

x − κk−1
Kn−1

)
,

68

T. YOSHIDA

I(a < x < b) is the indicator function of an interval (a, b) and Bp (x) is the pth
Bernoulli polynomial, see Zhou et al. [17]. For a random sequence Un , E[Un |X n ]and V [Un |X n ] designate the conditional expectation and the conditional varianceof Un given (X 1 , · · · , X n ) = (x1 , · · · , xn ), respectively.
In spline smoothing, the smoothing parameter λjn is usually selected as λjn → ∞with n → ∞ because a spline curve often yields overfitting for large n. In thefollowing, we assume that λjn = o(nKn−1 ).
In this section, first we discuss the asymptotic distribution of the RCPS. Next weshow that the diﬀerence of the BPS and RCPS asymptotically vanishes, by which
we finally obtain the asymptotic distribution of the BPS.
3.1. Asymptotic distribution of the RCPS. We will show the asymptotic
property of fˆj,γ (xj )(j = 1, · · · , D). By using the result of the partitioned matrix of Z ′ Z + Qm (λn ) + γI and its asymptotic property, the form of fˆj,γ (xj ) and
its asymptotic expression can be clarified. As a result, we obtain the following
Theorem.
Theorem 3.1. Let fj ∈ C p+1 (j = 1, · · · , D). Suppose that Kn = o(n1/2 ) and
λjn = o(n/Kn )(j = 1, · · · , D). Then, for j, k = 1, · · · , D,E[fˆj,γ (xj )|X n ] − fj (xj ) = bj,λ (xj ) + bj,γ (xj ) + oP (Kn−1 ) + oP (λjn Kn n−1 ),1−1−1V [fˆj,γ (xj )|X n ] =B(xj )′ G−1
j Σj Gj B(xj )(1 + oP (1)) = OP (Kn n ),nˆˆCov(fj,γ (xj ), fk,γ (xk )) = OP (n−1 ), j ̸= k
wherebj,λ (x) = −
and

λjn
∗
−1
B(xj )′ G−1
j Qm bj = O(λjn Kn n )
n

γ
∗
−1
bj,γ (xj ) = − B(xj )′ G−1
j bj = O(Kn n ).
n

In Theorem 3.1, the influence of γ appears only in bj,γ (xj ), which is in fact ofnegligible order. Furthermore, compared to Theorem 2 of Claeskens et al. [3],fˆj,γ (xj ) is asymptotically equivalent to the penalized spline estimator based on
the dataset {(yi , xij ) : i = 1, · · · , n} in the univariate regression model. By using
Theorem 3.1 and Lyapunov’s condition of the central limit theorem, we obtain theasymptotic joint distribution of [fˆ1,γ (x1 ) · · · fˆD,γ (xD )]′ .
Theorem 3.2. Suppose that there exists δ ≥ 2 such that E[|εi |2+δ |Xi = xi ] < ∞and fj ∈ C p+1 . Furthermore, Kn and λjn satisfy Kn = o(n1/2 ), n1/3 = o(Kn ) and
λjn = o((nKn−1 )1/2 ). Then, for any fixed point x ∈ (0, 1)D , as n → ∞,fˆ1,γ (x1 ) − f1 (x1 )
√
n 
 d
..

 −→ ND (0, Ψ) ,.
KnˆfD,γ (xD ) − fD (xD )ASYMPTOTICS FOR PENALIZED SPLINES IN ADDITIVE MODELS

69

where Ψ = diag[ψ1 (x1 ) · · · ψD (xD )] and1−1
B(xj )′ G−1
j Σj Gj B(xj ), j = 1, · · · , D.
n→∞ Knψj (xj ) = lim

3.2. Asymptotic distribution of the BPS. It is easily confirmed that the BPS
fˆj (xj ) can be expressed as(1)fˆj (xj ) = fj (xj ) + lim

ℓ→∞

ℓ−1
∑

(k+1)

{fj

(k)

(xj ) − fj (xj )}.

k=1

(k+1)
(k)
We show that for all k ∈ N, fj(xj ) − fj (xj ) asymptotically vanishes, from(1)
which we find that fˆj (xj ) is asymptotically dominated by fj (xj ) for j = 1, · · · , D.

These properties are summarized in the following two Propositions.
Proposition 3.3. Let fj ∈ C p+1 . Suppose that Kn = o(n1/2 ) and λjn = o(Kn n−1 ).
Then, as n → ∞,
(
)
(k+1)
(k)
fj
(xj ) − fj (xj ) = OP (Kn−k ) + OP Kn−(k−1) (Kn n)−1/2 , k = 1, 2, · · · .
Proposition 3.4. Under the same assumption as Proposition 3.3, as n → ∞,
)
(√
K
n
(1)
.
fˆj (xj ) = fj (xj ) + OP (Kn−1 ) + oP
n
(1)

(1)

By its simple form of fj (xj ), it is easy to show that fj (xj ) is asymptoticallyequivalent to fˆj,γ (xj ). Thus, Propositions 3.3 and 3.4 yield that the asymptotic
equivalence between fˆj (xj ) and fˆj,γ (xj ). Consequently, Theorem 3.2 implies theasymptotic distribution of BPS summarized as follows:
Theorem 3.5. Under the same assumption as Theorem 3.2, for any fixed pointx ∈ (0, 1)D , as n → ∞,fˆ1 (x1 ) − f1 (x1 )
√
n 
 d
..

 −→ ND (0, Ψ) ,.
KnfˆD (xD ) − fD (xD )where Ψ is that given in Theorem 3.2.
From Theorem 3.5, for i ̸= j, fˆi (xi ) and fˆj (xj ) are asymptotically independent.
Asymptotic normality and the independence of fˆi (xi ) and fˆj (xj ) in kernel smoothing also hold, as shown in Wand [12]. Thus, the penalized spline estimator andthe kernel estimator for the additive model have the same asymptotic property.
Asymptotic normality of ŷ in (4) can be shown as a direct consequence of Theorem
3.5. Though the BPS depends on the initial value, the eﬀect of the initial value on
the distribution of the BPS vanishes as n → ∞, which means that the uniquenessof the BPS is asymptotically satisfied. Furthermore, Theorem 3.6 indicates thatb̂ = b(∞) minimizes L(b).

70

T. YOSHIDA

Theorem 3.6. Let H(L) be the Hessian matrix of L(b). Then, H(L) is asymptotically positive definite.
We now give the optimal order of Kn and λjn in the context of minimizationof MSE of fˆj (xj ). Note that Kn controls the trade-oﬀ between the bias and thevariance of the estimator, both of which have been obtained by Proposition 3.4 andTheorem 3.1.
Corollary 3.7. Under the same assumption as Proposition 3.3, it follows that
({
(
}2 )
)
K
λ
Kn
1
n
jn
2
ˆ
ˆ
MSE(fj (xj )) = E[{fj (xj ) − fj (xj )} |X n ] = OP
+
+ OP
.
n
Kn
n
Furthermore taking Kn = O(n1/3 ) and λjn = O(nν ), ν ≤ 1/3 leads to the rate ofconvergence, MSE(fˆj (xj )) = OP (n−2/3 ).
Asymptotic normality in Theorem 3.2 also holds even if the optimal orders Kn =O(n1/3 ) and λjn = O(nν ), ν ≤ 1/3 in Corollary 3.7 are utilized, however the meanof the asymptotic distribution is not zero in such cases, that is, the bias terms given
in Theorem 3.1 do not vanish. We see that the centered asymptotic distribution canbe obtained by the assumptions for the orders of Kn and λjn as given in Theorem
3.2.
4. Discussion
In this paper, the asymptotic behavior of the penalized spline estimators in
the additive models was investigated. The BPS and the RCPS have been shown
to be asymptotically equivalent. For practical purposes, we compare the RCPS
with the BPS from the view point of computation. We have to calculate theinverse of MD = Z ′ Z + Qm (λn ) + γI in order to obtain the RCPS. Since the sizeof MD is {D(Kn + p)} × {D(Kn + p)}, the computation of MD−1 is O({D(Kn +p)}3 ). On the other hand, the ℓ-stage backfitting algorithm requires O(ℓD{Kn +p})
computations. Therefore when D is large, the BPS can be computed more quickly
than the RCPS. Even for D = 2, the RCPS requires O(8(Kn + p)3 ) computations
which is larger than the BPS. The same conclusion holds for a large sample size n,which is also detailed in Hastie and Tibshirani [7] as the advantage of the backfitting
algorithm. So it might be better to utilize the BPS in additive penalized splinesmoothing.
On the other hand, it is known that the BPS can be constructed by using a
blend of backfitting and scoring algorithms also in the generalized additive models(GAM). However, since the backfitting algorithm is required to be implemented
within each iteration of the scoring algorithm, the computation of the BPS finallybecomes a heavy task. Therefore in GAM, the RCPS, the direct method without
the backfitting cycle, might be better and its asymptotic theory should be investigated. Actually, Yoshida and Naito [16] showed that the RCPS has the asymptotic
normality in GAM. Such theoretical research is a generalization of the works bythis paper and Kauermann et al. [8].

ASYMPTOTICS FOR PENALIZED SPLINES IN ADDITIVE MODELS

71

Acknowledgements
Almost all of the work was done while the author was visiting the Graduate
School of Science and Engineering, Shimane University. The author would liketo thank Prof. Kanta Naito for his suggestions, encouragements and insightfulcomments.
References[1] Aerts, M., Claeskens, G. and Wand, M.P. Some theory for penalized spline generalized
additive models. Journal of Statistical Planning and Inference, 103 (2002), 455-470.
[2] Buja,A.Hastie,T. and Tibshirani,R. Linear smoothers and additive models (with discussion).
Ann. Statist. 17 (1989),453-555.
[3] Claeskens,G., Krivobokova,T. and Opsomer,J.D. Asymptotic properties of penalized spline
estimators. Biometrika. 96 (2009), 529-544.
[4] de Boor,C. A Practical Guide to Splines. Springer-Verlag, 2001.
[5] Eilers,P.H.C. and Marx,B.D. Flexible smoothing with B-splines and penalties(with Discussion). Statist.Sci. 11 (1996), 89-121.
[6] Hall,P. and Opsomer,J.D. Theory for penalized spline regression. Biometrika. 92 (2005),
105-118.
[7] Hastie,T. and Tibshirani,R. Generalized Additive Models. London: Chapman & Hall, 1990.
[8] Kauermann,G., Krivobokova,T., and Fahrmeir,L. Some asymptotic results on generalized
penalized spline smoothing. J. R. Statist. Soc. B 71 (2009), 487-503.
[9] Marx,B,D. and Eilers,P.H.C. Direct generalized additive modeling with penalized liklihood.
Comp. Statist & Data Anal. 28 (1998), 193-209.
[10] O’Sullivan,F. A statistical perspective on ill-posed inverse problems(with discussion). Statist.
Sci. 1 (1986), 505–27.
[11] Ruppert,D., Wand,M.P. and Carroll,R.J. Semiparametric Regression. Cambridge University
Press, 2003.
[12] Wand,M.P. A central limit theorem for local polynomial backfitting estimators. J. Mult.
Anal. 70 (1999), 57–65.
[13] Wang,L. and Yang,L. Spline-backfitted kernel smoothing of nonlinear additive auto regression
model. Ann. Statist. 35 (2007), 2474-2503.
[14] Wang,X., Shen,J. and Ruppert,D. On the asymptotics of penalized spline smoothing. Ele. J.
Statist. 5 (2011), 1-17.
[15] Yoshida,T. and Naito,K. Asymptotics for penalized additive B-spline regression. J. Japan.
Statist. Soc. 42 (2012), 81-107.
[16] Yoshida,T. and Naito,K. Asymptotics for penalized splines in generalized additive models.
Revised. (2013).
[17] Zhou,S., Shen,X. and Wolfe,D.A. Local asymptotics for regression splines and confidence
regions. Ann. Statist. 26 (1998), 1760-1782.
Takuma Yoshida: Graduate School of Science and Engineering, Kagoshima University, Kagoshima, Kagoshima, 890-8580, Japan
E-mail address: yoshida@sci.kagoshima-u.ac.jp