ESL读书笔记：Chap 5 Basis Expansions and Regularization

Chap 5 Basis Expansions and Regularization

Introduction

首先回顾以下线性模型。线性模型建模便于估计参数，并且可以看成是真正模型的一阶泰勒展开，因此具有较好的近似性。本章的主要内容是将X进行变换改变输入，再利用相加的形式，从而拟合非线性模型。具体来说，对X进行$h_m(X):{\mathbf{R}}^p\mapsto \mathbb{R},m=1,...,M$，再拟合模型：
$$f(X)=\sum _{m=1}^M{\beta }_m{h}_m(X)$$
这样处理好处在于进行参数估计的时候还可以使用原来的一整套方法。

常用的变换有：

1. $h_m(X)=X_m,m=1,\mathrm{2...},p$，也就是使用选模型；
2. $h_m(X)=X_J^2$或者$h_m(X)=X_j{X}_k$也就是多项式的形式，从而作为对f(.)的高阶泰勒展开
3. $h_m(X)=\log (X_j),\sqrt{X_j}...$等常见的对于单变量的非线性变换
4. $h_m(X)=I(L_m\le X_k<U_m)$也就是表示$X_j$的取值区间

那么接下来的问题就在于如何选择合适的变换，并且使用一定的方法来控制模型的复杂度。常见的控制模型复杂度的方法：

1. 控制每一个变换$h_m(.)$涉及的变量数
2. 变量选择，例如选模型，使用stagewise greedy approaches（例如CART,boosting）
3. 正则化

Piecewise Polynomials and Splines

使用分段多项式作为basis function，其中的分段点称作knots。使用更高阶的多项式可以得到更加光滑的basis function，也就是保证在节点处的连续型。

定义：cubic spline【三次样条】
分段将区间[a,b]分成n个区间，$[(x_0,x_1),...,(x_{n-1},x_n)]$，共有n+1个点，其中两个端点$x_0=a,x_n=b$，每个小区间的曲线是一个三次多项式，并且满足条件：

1. 在每个分段小区间上$[x_i,x_{i+1}]$，$S(x)=S_i(x)$都是一个三次多项式
2. 在节点处是2阶导数连续的

【注】：假设将区间划分成n个区间，那么以下的cubic spline对应的basis function可以是
$$\begin{array}{rl}{h}_1(X)& =1,h_3(X)=X^2,h_5(X)=(X-{\xi }_i{)}_{+}^3,\\ h_2(X)& =X,h_4(X)=X^3,h_6(X)=(X-{\xi }_j{)}_{+}^3.\end{array}$$
所以需要估计的参数值为：$n\times 4-(n-1)\times 3$，前面是不加任何光滑性限制的参数估计个数，后面减去的是每一个节点由于是二阶导连续的，所以函数值、一阶导数值、二阶导数值都需要是相等的，即每一个节点都有三个方程限制，所以要减去$(n-1)\times 3$。

（定义）order-M spline:分段多项式是M-2阶导数连续的
【注】：

1. cubic spline是2阶导数连续的，M=4；如果是阶梯函数，则M=1，即order-1 spline；如果是分段线性函数，即0阶连续，M=2；
   2. M-order spline，对应的多项式的阶数为至少为M-1阶
   3. M-order有K个knots，自由度:
   $$df=(K+1)M-K(M-1)=K+M$$
   对应了K+M个basic truncated-power function：
   $$\begin{array}{rcl}{h}_j(X)& =& X^{j-1},j=1,\dots ,M,\\ h_{M+\ell }(X)& =& (X-{\xi }_{\ell }{)}_{+}^{M-1},\ell =1,\dots ,K.\end{array}$$
   实际上估计参数的时候就是抓住这K+M个function

一般来说，不会使用超过三次样条的多项式了，一般取M=1,2,4

固定konts的样条被称作regression splines,进行这样的操作需要选择样条的次数、节点的个数、节点的位置。

Natural Cubic Splines

使用多项式拟合数据在边界会出现方差较大的情况，这对于样条更为突出相较于全局的多项式。多项式拟合出现边界点方差较大，很大程度上是因为泰勒展开并不收敛，所以实际上在边界位置多项式拟合的方差往往会呈现阶数越大，方差越大的情况。所以一个很自然的想法是在靠近边界点，也就是边界的konts之外的部分使用线性函数，内部仍然使用多项式样条。这样做会减少4个参数，$df=K+4-4=K$，这个自由度可以安排进在内部放更多的konts。相对来说可以有以下几个basic truncated power functions:
$$N_1(X)=1,N_2(X)=X,N_{k+2}(X)=d_k(X)-d_{K-1}(X)$$
$$d_k(X)=\frac{(X-{\xi }_k{)}_{+}^3-(X-{\xi }_K{)}_{+}^3}{{\xi }_K-{\xi }_k}.$$
由于natural cubic splines还需要保持三次导数连续所以本质上仍然是：
$$f(x)=\sum _{j=0}^3{\beta }_j{x}^j+\sum _{k=1}^K{\theta }_k(x-{\xi }_k{)}_{+}^3$$
只是再加上$x<x_1,x>x_K$的时候需要是线性的函数。

**使用样条插值即可以作为增加复杂度，使得拟合效果更好；也可以作为降维的方法，进一步提取信息。**例如对于函数型数据的线性建模：$Y={\sum }_{j}x(t_j)\beta (t_j)$，为了使得$\beta (t)\in {\Re }^{T\times 1}$关于$t$是光滑的，可以使用spline的方法进行建模，那么$\beta ={\sum }_i{h}_i(t){\theta }_i=H\theta ,H\in {\Re }^{T\times M}$（H只和已知的T有关），由此转化为$Y=XH\theta$进行建模，将$X^{\ast }=XH$视作整体即可达到降维的效果。

Smoothing Splines

将所有点都作为konts，并且使用正则化的方法来减少复杂度，损失函数为：
RSS ⁡ ( f , λ ) = ∑ i = 1 N { y i − f ( x i ) } 2 + λ ∫ { f ′ ′ ( t ) } 2 d t \operatorname{RSS}(f,\lambda)=\sum_{i=1}^N\{y_i-f(x_i)\}^2+\lambda\int\{f''(t)\}^2dt RSS(f,λ)=i=1∑N​{yi​−f(xi​)}2+λ∫{f′′(t)}2dt
易知，当$\lambda =0$，则f可以是任何一个插值函数；当$\lambda \to +\infty$，则为任何一条直线。随着$\lambda$的增大，f越来越光滑。如果只看损失函数最小化，实际上我们找的是一个函数，参数空间应当是无穷维的，但是可以证明，最小化上式的一定是cubic spline，此时前一项为0，后一项是所有插值函数中最小的。

# 证明：

首先证明对于任何一个二阶可导的插值函数$\tilde{g}(x)$，令$h(x)=\tilde{g}(x)-g(x)$，则使用分部积分有：
∫ a b g ′ ′ ( x ) h ′ ′ ( x ) d x = − ∑ j = 1 N − 1 g ′ ′ ′ ( x j + ) { h ( x j + 1 ) − h ( x j ) } = 0 \begin{aligned} \int_{a}^{b}g^{\prime\prime}(x)h^{\prime\prime}(x)dx& =\quad-\sum_{j=1}^{N-1}g'''(x_j^+)\{h(x_{j+1})-h(x_j)\} \\ & =0 \end{aligned} ∫ab​g′′(x)h′′(x)dx​=−j=1∑N−1​g′′′(xj+​){h(xj+1​)−h(xj​)}=0​
从而
$${\int }_a^b{\tilde{g}}^{\prime \prime }(t{)}^{2}dt\ge {\int }_a^b{g}^{\prime \prime }(t{)}^{2}dt,$$
由此对于任何一个二阶可导的函数$\tilde{f}(x)$，可以选择在$x_i$与$\tilde{f}(x)$一致的cubic splines$f(x)$，由此
$$\begin{array}{r}\sum _{i=1}^N(y_i-f(x_i){)}^2\end{array}$$
是一致的，而前面的证明又说明：
$$J(f(x))\le J(\tilde{f}(x))$$
所以使用$x_i$作为knots的cubic splines是使得损失函数最小的。#

因此可以将$f(x)$表示成如下形式：
$$f(x)=\sum _{j=1}^N{N}_j(x){\theta }_j$$
相应的，原损失函数也可以写成：
$$\mathrm{RSS}(\theta ,\lambda )=(\mathbf{y}-\mathbf{N}\theta {)}^T(\mathbf{y}-\mathbf{N}\theta )+\lambda {\theta }^T{\mathbf{\Omega }}_N\theta ,$$
其中 { N } i j = N j ( x i ) ∈ ℜ N × N \{\mathbf{N}\}_{ij}=N_{j}(x_{i})\in \Re^{N\times N} {N}ij​=Nj​(xi​)∈ℜN×N是满秩矩阵, { Ω N } j k = ∫ N j ′ ′ ( t ) N k ′ ′ ( t ) d t . \{\boldsymbol{\Omega}_{N}\}_{jk}=\int N_{j}^{\prime\prime}(t)N_{k}^{\prime\prime}(t)dt. {ΩN​}jk​=∫Nj′′​(t)Nk′′​(t)dt.因此系数估计值为：
$$\hat{\theta }=({\mathbf{N}}^T\mathbf{N}+\lambda {\mathbf{\Omega }}_N{)}^{-1}{\mathbf{N}}^T\mathbf{y},$$
最终得到的样条函数为：
$$\hat{f}(x)=\sum _{j=1}^N{N}_j(x){\hat{\theta }}_j.$$

Degrees of Freedom and Smoother Matrices

介绍一种不使用交叉验证直接选择$\lambda$的方法：linear smoother

从前面的推到可以得到：
$$\begin{array}{rcl}\hat{\mathbf{f}}& =& \mathbf{N}({\mathbf{N}}^T\mathbf{N}+\lambda {\mathbf{\Omega }}_N{)}^{-1}{\mathbf{N}}^T\mathbf{y}\\ & =& {\mathbf{S}}_{\lambda }\mathbf{y}.\end{array}$$
再次可以看到估计值是Y的线性函数，其中$S_{\lambda }$被称作smoother matrix，它只和X、$\lambda$有关，和一般的最小二乘估计的帽子矩阵H对比来说有以下的特点：

1. 都是半正定对称矩阵
2. $H^2=H$是幂等阵，${\mathbf{S}}_{\lambda }{\mathbf{S}}_{\lambda }⪯{\mathbf{S}}_{\lambda }$,也就是$trace({\mathbf{S}}_{\lambda }{\mathbf{S}}_{\lambda })\le trace(S_{\lambda })$

和之前的岭回归一致，定义自由度：
$$df:=trace(S_{\lambda })$$
现在重新改写矩阵$S_{\lambda }$：
$$\begin{array}{rl}{S}_{\lambda }=& \mathbf{N}({\mathbf{N}}^T\mathbf{N}+\lambda {\mathbf{\Omega }}_N{)}^{-1}{\mathbf{N}}^T\\ =& N[N^T(I+\lambda (N^{-1}{)}^T{\mathbf{\Omega }}_N{N}^{-1})N{]}^{-1}{N}^T\\ =& (I+\lambda K{)}^{-1}\end{array}$$
其中$K=(N^{-1}{)}^T{\mathbf{\Omega }}_N{N}^{-1}$只和X有关和$\lambda$无关，此时可以把估计值写为$\hat{f}=(I+\lambda K{)}^{-1}Y$，是以下问题的解：
$$\underset{\mathbf{f}}{\min }(\mathbf{y}-\mathbf{f}{)}^T(\mathbf{y}-\mathbf{f})+\lambda {\mathbf{f}}^T\mathbf{Kf}$$
所以也把K称作penalty matrix。K是一个对称阵，对角化$K=UD{U}^T$，其中$U=[u_1,...,u_N]$，$D=diag(d_1,...,d_N)$，从而有：
$${\mathbf{S}}_{\lambda }=\sum _{k=1}^N{\rho }_k(\lambda ){\mathbf{u}}_k{\mathbf{u}}_k^T$$
$${\rho }_k(\lambda )=\frac{1}{1+\lambda d_k}$$
$$\hat{Y}=S_{\lambda }Y=\sum _{k=1}^N{\rho }_k(\lambda ){\mathbf{u}}_k{\mathbf{u}}_k^{T}Y=\sum _{k=1}^N{\rho }_k(\lambda ){\mathbf{u}}_k<{\mathbf{u}}_k,Y>$$
可以看到$u_k$和$\lambda$无关，本质上$\lambda$就是做了一个shrinkage，因此smoothing splines也被称作shrinking smoothers, 一般的splines被称作projection smoothers.针对$u_k$和${\rho }_k(\lambda )$有以下的讨论：

1. $u_k$并不会被$\lambda$影响(\emph{事实上K和$\lambda$无关})
2. ${\mathbf{S}}_{\lambda }\mathbf{y}={\sum }_{k=1}^N{\mathbf{u}}_k{\rho }_k(\lambda )⟨{\mathbf{u}}_k,\mathbf{y}⟩$可以把${\rho }_k(\lambda )$看成一个shrinkage，并且有$u_k$越复杂，则${\rho }_k(\lambda )$越小，也就是越复杂的shrink的越多
3. shrinkage的速率由$\lambda$控制，$\lambda$越大，则shrink的越多
4. 原始的参数估计等价于：
   $$\underset{\theta }{\min }\Vert \mathbf{y}-\mathbf{U}\theta {\Vert }^2+\lambda {\mathbf{\theta }}^T\mathbf{D}\mathbf{\theta },$$
   其中$D=diag(d_k)$，则有$\hat{\theta }=(U^{T}U+\lambda D{)}^{-1}{U}^{T}Y=U(I+\lambda D{)}^{-1}{U}^{T}Y$
5. ${\rho }_k(\lambda )$中一定有两个等于1，也就是K的特征值中一定有两个等于零，事实上：
   $$K=(N^{-1}{)}^T{\Omega }_N{N}^{-1}$$
   由于N是满秩的，现在看${\Omega }_N=\int N_j^{\prime \prime }{N}_k^{\prime \prime }ds$，有两列两行对应$N_1(x)=1,N_2(x)=x$所以为0，因此肯定有两个特征值为0，说明这两个是不会被shrink的
6. $\lambda \to 0,d{f}_{\lambda }\to N,S_{\lambda }\to I;\lambda \to \infty ,d{f}_{\lambda }\to 2,S_{\lambda }\to H$，H是线性回归的帽子矩阵

Automatic Selection of the Smoothing Parameters

承接上文如何选择$\lambda$的问题，一个和显然的想法是先定$d{f}_{\lambda }$
$$d{f}_{\lambda }=\sum _{i=1}^N\frac{1}{1+\lambda d_k}$$
由于是$\lambda$的单调函数，所以可以数值求解。而定$d{f}_{\lambda }$则可以使用最小化EPE(integrated squared prediction error)来决定。需要注意的是关于$d{f}_{\lambda }$存在一个trade-off，df越大偏差越小，但是方差越大。
$$\begin{array}{rl}\mathrm{EPE}({\hat{f}}_{\lambda })& =E(Y-{\hat{f}}_{\lambda }(X){)}^2\\ & =\mathrm{Var}(Y)+\mathrm{E}\left[{\mathrm{Bias}}^2({\hat{f}}_{\lambda }(X))+\mathrm{Var}({\hat{f}}_{\lambda }(X))\right]\\ & ={\sigma }^2+\mathrm{MSE}({\hat{f}}_{\lambda }).\end{array}$$

【注】：EPE并不是MSE！但由于${\sigma }^2$无法控制，所以本质上还是最小化MSE。

由于EPE无法计算，所以常常使用交叉验证的方法来估计，从而选择最小损失下的df，再确定$\lambda$即可。

Nonparametric Logistic Regression

非参数的方法实际上并不是和参数方法一样，实现给定包含参数的拟合函数形式，再根据最小化损失函数得到参数估计值。非参数方法是在实现给定的损失函数下，推导最小化损失函数的拟合函数是什么，并且往往这个函数需要具有一定的光滑性，常可以通过推导最小化损失函数得到这个拟合函数的形式，从而使用参数化的方法具象表示出来，再结合数据估计即可。例如对于natural cubic splines来说，就是在残差平方和以及$J=\int (f^{{}^{\prime \prime }}{)}^{2}dx$进行的，由此得到natural cubic splines的形式，在此基础上得到“自由度拉满”的basis function，$\tilde{X}=g(X)$转化为一般的线性回归建模。

例如本节的非参数二分类逻辑回归就是假设：
$$\log \frac{\mathrm{Pr}(Y=1|X=x)}{\mathrm{Pr}(Y=0|X=x)}=f(x),$$
$$\mathrm{Pr}(Y=1|X=x)=\frac{e^{f(x)}}{1+e^{f(x)}}.$$
损失函数一般包括两项，一项是衡量准确度，另一项是控制复杂程度：
ℓ ( f ; λ ) = ∑ i = 1 N [ y i log ⁡ p ( x i ) + ( 1 − y i ) log ⁡ ( 1 − p ( x i ) ) ] − 1 2 λ ∫ { f ′ ′ ( t ) } 2 d t = ∑ i = 1 N [ y i f ( x i ) − log ⁡ ( 1 + e f ( x i ) ) ] − 1 2 λ ∫ { f ′ ′ ( t ) } 2 d t \begin{array}{rcl}\ell(f;\lambda)&=&\sum_{i=1}^N[y_i\log p(x_i)+(1-y_i)\log(1-p(x_i))]-\frac12\lambda\int\{f''(t)\}^2dt\\\\&=&\sum_{i=1}^N\left[y_if(x_i)-\log(1+e^{f(x_i)})\right]-\frac12\lambda\int\{f''(t)\}^2dt\end{array} ℓ(f;λ)​==​∑i=1N​[yi​logp(xi​)+(1−yi​)log(1−p(xi​))]−21​λ∫{f′′(t)}2dt∑i=1N​[yi​f(xi​)−log(1+ef(xi​))]−21​λ∫{f′′(t)}2dt​
类似于smoothing natural splines的推导，这里也可以推导得到f(.)是smoothing natural splines，即$f(x)={\sum }_{i=1}^N{N}_j(x){\theta }_j$。由于一阶导数不是$\theta$的线性函数，无法直接得到解析解，可以使用Newton-Raphson迭代。

Multidimensional Splines

以上讨论的都是一维的splines，本届推广到多维度的场合。从一维的basis function出发，不难联想到以下两种最直接的推广：

1. tensor product basis。以两个维度为例：
   $$g_{jk}(X)=h_{1j}(X_1)h_{2k}(X_2),j=1,\dots ,M_1,k=1,\dots ,M_2$$
   从而进一步建模为：
   $$g(X)=\sum _{j=1}^{M_1}\sum _{k=1}^{M_2}{\theta }_{jk}{g}_{jk}(X).$$
2. additive splies。以两个维度为例:
   $$g(X)=\sum _{j=1}^{M_1}{h}_{1j}(X_1)+\sum _{i=1}^{M_2}{h}_{2i}(X_2)$$

对比来说以上两种方法的区别在于：product的方法是更加flexible的$d{f}_{prod}={\Pi }_{i=1}^p{M}_i,d{f}_{sum}={\sum }_{i=1}^p(M_i-1)+1$；但是随着维度的增加参数估计的数量是呈现指数倍增长的。

类似于前面的非参数逻辑回归，多维的、不同basis function下的、不同组合的spline也可以看成是某一个损失下的最优函数形式，例如thine-plate spline：
min ⁡ f ∑ i = 1 N { y i − f ( x i ) } 2 + λ J [ f ] , \min_f\sum_{i=1}^N\{y_i-f(x_i)\}^2+\lambda J[f], fmin​i=1∑N​{yi​−f(xi​)}2+λJ[f],
以二维为例：
$$J[f]=\int {\int }_{{\mathbb{R}}^2}\left[{\left(\frac{{\partial }^{2}f(x)}{\partial x_1^2}\right)}^2+2{\left(\frac{{\partial }^{2}f(x)}{\partial x_1\partial x_2}\right)}^2+{\left(\frac{{\partial }^{2}f(x)}{\partial x_2^2}\right)}^2\right]d{x}_{1}d{x}_2.$$
以上损失函数的解具有形式：
$$f(x)={\beta }_0+{\beta }^{T}x+\sum _{j=1}^N{\alpha }_j{h}_j(x),h_j(x)=||x-x_j|{|}^2\log ||x-x_j||$$
其中$h_j$被称作radial basis function。

而对于多维下的thin-plate spline的计算复杂度，一维的thin-plate spline是$O(N^3)$，高维度下可以不完全使用N个点作为Knots，减少为K个knots，计算复杂度为$O(N{K}^2+K^3)$。

对于additive spline models，也可以找到对应的正则函数，但更简单的表达是直接假设$f(X)=\alpha +f_1(X_1)+...+f_d(X_d)$，并设
$$\begin{array}{rcl}J[f]& =& J(f_1+f_2+\cdots +f_d)\\ & =& {\sum }_{j=1}^d\int f_j^{\prime \prime }(t_j{)}^{2}d{t}_j.\end{array}$$

ANOVA spline decompositions：
$$f(X)=\alpha +\sum _j{f}_j(X_j)+\sum _{j<k}{f}_{jk}(X_j,X_k)+\cdots ,$$
这就会涉及一些讨论：

1. 最大的交叉项阶数
   2. 需要包含哪些部分，哪些主要部分，哪些交互项
   3. 需要使用怎样的$f_{.}(.)$并且同时考虑到参数估计的数量

补充：Wavelet Smoothing

本质上是使用正交的（相乘积分等于0）的basis来表示函数，再使用shrink或选模型的方法来得到sparse representation.它可以有效地表示平滑和局部的波动，这一部分是选读部分，放一张图}来大概看一下basis function长什么样：