Spline Regression with Penaties

1.前言

前面介绍spline的基本性质和最小二乘回归，这里进入正题，主要总结Spline Regression的内容，主要的复现参考论文中的P-spline,针对的是Cardinal B-spline，就是分割节点是等间隔的

2.基本理论

高阶基函数$B_{i,k}(x)$的递推公式
$$B_{i,k}(x):=\frac{x-t_i}{t_{i+k}-t_i}{B}_{i,k-1}(x)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}}{B}_{i+1,k-1}(x).$$
$B_{i,k}(x)$的一阶导数
$$\frac{d{B}_{i,k}(x)}{dx}=k\left(\frac{B_{i,k-1}(x)}{t_{i+k}-t_i}-\frac{B_{i+1,k-1}(x)}{t_{i+k+1}-t_{i+1}}\right)$$
在Cardinal B-spline时，假设间隔为$h$,则有
$$\frac{d{B}_{i,k}(x)}{dx}=\frac{1}{h}\left(B_{i,k-1}(x)-B_{i+1,k-1}(x)\right)$$
B-spline的加权最小二乘法
$$\hat{\alpha }=\underset{\alpha }{\mathrm{argmin}}\sum _{\text{all}x}{\left\{W(x)\left[y(x)-\sum _i{\alpha }_i{B}_{i,k,t}(x)\right]\right\}}^2$$
$W(x)$是一个权值系数，这里取单位矩阵$I$,${\alpha }_i$为待定的最小二乘系数

3惩罚项

3.1惩罚项介绍

为了改进这些最小二乘系数呢，在最小二乘的目标式子中，加入惩罚项。不同的惩罚项有不同的效果。
从前面的最小二乘法我们可以得到B-spline的曲线拟合$\hat{y}$
$$\hat{y}=\sum {\alpha }_i{B}_{i,k}(x)$$
惩罚项可以从$\hat{y}$的导数连续程度，以及$\alpha$的连续程度两个来入手
$\hat{y}$的一阶导数
$$\frac{d\hat{y}}{dx}=\frac{d}{dx}\sum _i{\alpha }_i{B}_{i,k}=\sum _i{\alpha }_i{B}_{i,k}^{\prime }=-\frac{1}{h}\sum _j\Delta {\alpha }_{j+1}{B}_j(x;q-1)$$
$\hat{y}$的二阶导数

$$\begin{gathered} \sum _j{a}_j{B}_j^{\prime \prime }(x;q)=h^{-2}\sum _j{\Delta }^2{a}_j{B}_j(x;q-2) \\ {\Delta }^2{a}_j=\Delta \Delta a_j=a_j-2a_{j-1}+a_{j-2} \end{gathered}$$

3.2惩罚项

3.2.1 二阶导数

$$\begin{array}{rl}S=& \sum _{i=1}^m{\left\{y_i-\sum _{j=1}^n{a}_j{B}_j\left(x_i\right)\right\}}^2+\lambda {\int }_{x_{\min }}^{x_{\max }}{\left\{\sum _{j=1}^n{a}_j{B}_j^{\prime \prime }(x)\right\}}^{2}dx\end{array}$$

对二阶导的积分，其实只是对二导求一个和。最终会以离散的形式计算，而不是以积分的形式。这个和越大，代表曲线越不光滑。
举例来说，以3阶B-spline的二阶导数作为惩罚项
$$\begin{gathered} h^{2}P=\lambda {\int }_{x_{\min }}^{x_{\max }}{\left\{\sum _j{a}_j{B}_j^{\prime \prime }(x;3)\right\}}^{2}dx=\lambda {\int }_{x_{\min }}^{x_{\max }}{\left\{\sum _j{\Delta }^2{a}_j{B}_j(x;1)\right\}}^{2}dx \\ =\lambda {\int }_{x_{\min }}^{x_{\max }}\sum _j\sum _k{\Delta }^2{a}_j{\Delta }^2{a}_k\cdot B_j(x;1)B_k(x;1)dx \end{gathered}$$
对于一阶B-spline的基函数$B_j(x;1)$,它的作用范围为只在$[t_j,t_{j+1})$,那么，$k$只有在$k\in [j-1,j+1]$，两者才存在overlap，其余皆为0.上式可以继续化简，但是最终的结果依旧比较复杂。这就是文献1，为什么要提出了系数高阶导的方法

3.2.2 系数高阶导

$$S=\sum _{i=1}^m{\left\{y_i-\sum _{j=1}^n{a}_j{B}_j\left(x_i\right)\right\}}^2+\lambda \sum _{j=k+1}^n{\left({\Delta }^k{a}_j\right)}^2$$
写成矩阵的矩阵形式如下
$$B^{T}y=\left(B^{T}B+\lambda D_k^T{D}_k\right)a$$
B基函数矩阵，每一列代表一个基函数产生的序列，$D_k$代表了系数$\mathbf{a}$的k阶差分矩阵，且有 $\Delta a^k=D_{k}a$.在计算$\Delta a^k$ 时，这个看似不好表达的惩罚项目，亏了优美的线性代数，这项还是容易表达的。只要对单位矩阵$I_k$做k次按行差分即可，matlab和R语言都有内置diff函数可以直接用。需要要注意的是，差分项随着差分阶数的增加而减少。如果没有$D_k$，这个方程就是最普通的最小二乘回归。
在论文采用的是系数的2阶导数，到了这里，我们已经可以可以利用$\lambda$控制拟合的平滑程度。如何选择一个合适的$\lambda$，需要有了评价指标。这篇文章选择了$Akaikeinformationcriterion(AIC)$指标，属于信息论里的概念，有许多种定义，目测还不等价，我引用了wiki中的定义。
$$\begin{gathered} AIC=2k-2L \\ L=-ln(\sqrt{(\mathrm{R}\mathrm{S}\mathrm{S}/\mathrm{n})}) \\ \Rightarrow AIC=2k+ln(RSS/n) \end{gathered}$$
$k$是参数数量，书中指的是有效参数数量，计算比想象得复杂得多，$RSS$残差，$n$代表样本数量。意义很好理解，显然AIC越小越好。
论文中的定义更加复杂，不知道是不是扯上信息论，都比较复杂
书中最终采用是GCV而不是直接使用AIC，其计算过程
$$H=B{\left(B^T\tilde{W}B+\lambda D_k^T{D}_k\right)}^{-1}{B}^T\tilde{W}$$
$\tilde{W}$为估计的加权系数
GCV计算过程
$$\mathrm{GCV}(\lambda )=\sum _{i=1}^m\frac{{\left(y_i-{\hat{y}}_i\right)}^2}{{\left(m-{\sum }_{i=1}^m{h}_{ii}\right)}^2}$$
这公式分子部分代表了残差，分母代表了误差的方差估计。论文假定误差是正态分布的。这公式中，$h_{ii}$体现出有效参数数量k.因为$\lambda$比较大的时候，系数会被压缩，它的有效长度会减少。为什么用这个来表征k呢，这个细节我没弄明白。
$$T(\lambda )=\sum _{i=1}^m{h}_{ii}=\mathrm{tr}(H)=\mathrm{tr}\left\{(I+\lambda L{)}^{-1}\right\}=\sum _{j=1}^n\frac{1}{1+\lambda {\gamma }_j}$$
${\gamma }_j$是$L$矩阵的特征向量，$L$矩阵是$H$的一部分，具体就不展开了
我感觉作者总想使算法变得更加generalized，所以无形中，感觉简单的部分变得复杂起来。

4.复现代码

论文中，附加了S-PLUS和matlab代码。所以，也谈不上说复现，他用到的测试数据我也没有去找，自己生成一个数据就完事了。
我就用他的S-PLUS代码，稍作修改一下，在R环境里面就可以跑了
作者的代码比 文章要简单太多，也暴露自己基础不扎实的问题

bspline<-function(x,xl,xr,ndx,bdeg){
  dx<-(xr-xl)/ndx
  knots<-seq(xl-bdeg*dx,xr+bdeg*dx,by=dx)
  B<-spline.des(knots,x,bdeg+1,0*x)$design
  B
}
pord <-2
n <- 100
x <- seq(0, 1, length=n)
knots <-seq(0, 1, length=10)
bdeg <-3
lambda <-100
dgp <- sin(2*pi*x)
y <- dgp + rnorm(n, sd=.2)
B <- spline.des(knots,x, bdeg,0,outer.ok = TRUE)
B <- B$design
D <- diag(ncol(B))
yhat<-NULL
diff(D,differences = pord)
gcv <- 1:10
lambdaArray <- 1:10

for(i in 1:10){
  lambda = (2^i-2)/50
  lambdaArray[i] = lambda
a <- solve(t(B) %*% B + lambda * t(D) %*% D,t(B) %*% y)
yhat<-cbind(yhat,B %*% a)
s <- sum((y - yhat)^2)
Q <- solve(t(B) %*% B + lambda * t(D) %*% D)
# matrix inversion
t <- sum(diag(Q %*% (t(B) %*% B)))
gcv[i] <- s / (nrow(B) - t)^2
}
par(mfrow = c(3,1))
matplot(x, B, type="l", lwd=2)
plot(x, y, cex=.25, col="grey")

lines(x, dgp,   lty=2)
lines(x, yhat[,1], col="red",lwd=0.5)
lines(x, yhat[,5], col="blue",lwd=0.5)
lines(x, yhat[,9], col="green", lwd=0.5)
legend("topright",legend=c("s+n","s","lb=0","lb=0.6","lb=10.2"),    col=c("grey","black","red","blue","green"),  lty=1,lwd=2)  
plot(lambdaArray, gcv,type = "b")

第一幅图是基函数，第二幅图是不同$\lambda$的拟合情况情况，$\lambda$越大，曲线越光整，简直到了平整，绿色的线到了往水平直线方向趋近了。第三幅图可以看到gcv曲线是单调上升的。也就是说这个数据其实并不需要做惩罚，本身已经达到很好的拟合效果。

5.小结

这篇论文写得不错，引用量很高。spline的回归和惩罚讲得比较清楚，特别是背后的原理部分，有许多值得借鉴。作者在写作的时候，感觉还是个小年轻，措辞能力上感觉还有待改进。但是可以看到作者的理论功力还是不错，有许多细节值得深挖下去，论文后面带的评价也比较高。由于研究的方向和作者的不太相同，不想花太多时间。不过花了这么几天功夫去学习，还是有一些收获。

References

1. Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing with B-splines and penalties (with comments and rejoinder). Statistical Science 11(2): 89-121.