5.4 平滑样条 | 王冠嵩

本节介绍一个样条基函数方法，它将所有样本点作为结点从而避免了结点选择问题。它通过正则化来控制拟合的复杂度。考虑以下这个问题：从所有的二次导数连续的函数 $f(x)$ 中，寻找那个可使惩罚（penalized）残差平方和最小化的函数：

$$\operatorname{RSS}(f, \lambda) = \sum_{i=1}^N \{y_i - f(x_i)\}^2 + \lambda\int\{f^{\prime\prime}(t)\}^2 dt \tag{5.9}$$

其中的 $\lambda$ 为一个固定的 平滑参数（smoothing parameter）。第一项衡量模型了（拟合结果）与数据的接近程度，而第二项是对函数的曲度进行惩罚，$\lambda$ 的作用则是对这两者的权衡。有两种极端的情况：

$\lambda=0$：最优解 $f$ 可以为任意穿过所有样本数据的曲线。
$\lambda=\infty$：由于不允许非零的二阶导数，最优解就是简单的最小二乘的线性拟合。

两个极端情况中，最优解的函数从非常平滑到非常粗糙，我们希望可以通过改变 $\lambda\in(0,\infty)$ 的值得到介于两者之间的一族有意义的函数。

准则函数（式 5.9）的定义域是一个无限维度的函数空间；实际上，它是所有在式 5.9 中第二项中有定义的函数的索伯列夫（Sobolev）空间。值得注意的是，可证明式 5.9 存在一个显式的、有限维度的、并且唯一的最小化解，即为以输入变量的所有唯一取值 $x_i,i=1,\dots,N$ 作为结点的自然三次样条函数（练习 5.7）。从表面上看起来这类函数仍然是过参数化（over-parametrized）的，其中可包含了 $N$ 个结点，意味着有 $N$ 个自由度。不过惩罚项会对样条函数的系数加以约束，使模型朝着线性拟合的方向有某种程度的收缩。

最小化解为自然样条函数，所以可将其写为：

$$f(x) = \sum_{j=1}^N N_j(x) \theta_j \tag{5.10}$$

其中的 $N_j(x)$ 是一组 $N$ 维的基函数，它们构建了这一族的自然样条函数（参考第 5.2.1 节和练习 5.4）。最小化准则函数可简化为：

$$\operatorname{RSS}(\theta, \lambda) = (\mathbf{y} - \mathbf{N}\theta)^T(\mathbf{y} - \mathbf{N}\theta) + \lambda \theta^T\mathbf{\Omega}_N\theta \tag{5.11}$$

其中的 $\{\mathbf{N}\}_{ij}=N_j(x_i)$，$\{\mathbf{\Omega_N}\}_{jk}=\int N’’_j(t)N’’_k(t) dt$。则容易看出最小化解可被视为一个广义的岭回归：

$$\hat{\theta} = (\mathbf{N}^T\mathbf{N} + \lambda \mathbf{\Omega}_N)^{-1} \mathbf{N}^T \mathbf{y} \tag{5.12}$$

拟合的平滑样条函数为：

$$\hat{f}(x) = \sum_{j=1}^N N_j(x) \hat{\theta}_j \tag{5.13}$$

本章附录中介绍了求解平滑样条的快速计算方法。

**图 5.6**：输出变量为青春期的脊椎骨质密度的相对变化，作为年龄的函数。对男性和女性分别拟合平滑样条模型，平滑参数 $\lambda = 0.00022$。其对应了大概 12 个自由度。

图 5.6 展示了在青春期的骨质密度（BMD）数据上拟合的平滑样条结果。输出变量为两次连续检查的脊椎骨质密度的相对变化（百分比），通常间隔一年。数据的颜色代表了性别，并分别拟合了两条曲线。这个简单的结果印证了女性的生长突增（growth spurt）比男性提前两年的结论。两个拟合中的平滑参数 $\lambda$ 都大概为 0.00022；下面会介绍如何选定这个参数。

5.4.1 自由度和平滑器矩阵

上文尚未说明如何选择平滑样条的参数 $\lambda$。本章的稍后章节会介绍使用诸如交叉验证的方法自动选择参数。本节阐述如何从直观含义来决定平滑程度参数。

在确定 $\lambda$ 后，平滑样条可视为一个线性平滑器（或线性算子）。这是因为在式 5.12 中的估计参数是 $y_i$ 的一个线性组合。记 $\hat{\mathbf{f}}$ 为在训练样本点 $x_i$ 处的拟合值 $\hat{f}(x_i)$ 的 $N$ 维向量，则：

$$\begin{align}\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} \tag{5.14} \end{align}$$

可见拟合值同样是 $\mathbf{y}$ 的线性组合，并且这个有限维度的线性算子 $\mathbf{S}_\lambda$ 被称为 平滑器矩阵（smoother matrix）。这个线性的性质，意味着从 $\mathbf{y}$ 到 $\hat{f}$ 的变换过程不依赖于 $\mathbf{y}$ 本身；也就是说 $\mathbf{S}_\lambda$ 只依赖于 $x_i$ 和 $\lambda$。

在传统的最小二乘拟合中这类线性算子被人熟知。假设 $\mathbf{B}_\xi$ 为 M 个三次样条基函数在 $N$ 个训练样本点 $x_i$ 上的取值矩阵，维度为 $N\times M$，$\xi$ 为一系列结点的集合，并且 $M \ll N$。则样条函数的拟合值向量为：

$$\begin{align}\hat{\mathbf{f}} &= \mathbf{B}_\xi(\mathbf{B}_\xi^T\mathbf{B}_\xi)^{-1} \mathbf{B}_\xi^T\mathbf{y} \\ &= \mathbf{H}_\xi \mathbf{y} \tag{5.15} \end{align}$$

其中的线性算子 $\mathbf{H}_\xi$ 为一个投影算子，在统计学中也称为 帽子矩阵（hat matrix）。两个矩阵 $\mathbf{H}_\xi$ 和 $\mathbf{S}_\lambda$ 之间有一些重要的相似和不同点：

两者都是对称的半正定矩阵。
$\mathbf{H}_\xi\mathbf{H}_\xi=\mathbf{H}_\xi$（幂等，idempotent）；而 $\mathbf{S}_\lambda\mathbf{S}_\lambda\preceq\mathbf{S}_\lambda$，意味着右手边超过左手边一个半正定矩阵。稍后会说明，这是由 $\mathbf{S}_\lambda$ 的收缩性质导致的。
$\mathbf{H}_\xi$ 的秩为 $M$，而 $\mathbf{S}_\lambda$ 的秩为 $N$。

表达式 $M=\operatorname{trace}(\mathbf{H}_\xi)$ 为投影空间的维度，同时也是基函数的个数，因此也是拟合中参数的个数。以此为类比，定义平滑样条的 有效自由度（effective degrees of freedom） 为：

$$\text{df}_\lambda = \operatorname{trace}(\mathbf{S}_\lambda) \tag{5.16}$$

即 $\mathbf{S}_\lambda$ 的对角线元素的和。这个很有用的定义以更直观和一致的方式来参数化平滑样条，以及其他的平滑器。例如图 5.6 中，若指定每个曲线的自由度为 $\text{df}_\lambda=12$，则通过 $\operatorname{trace}(\mathbf{S}_\lambda)=12$ 可从数值上解出对应的 $\lambda\approx 0.00022$。此外，这个定义有很多其他好处，下面介绍其中一部分。

由于 $\mathbf{S}_\lambda$ 为对称并且半正定矩阵，其存在实数的特征分解。方便起见，将 $\mathbf{S}_\lambda$ 写为 Reinsch 形式：

$$\mathbf{S}_\lambda = (\mathbf{I} + \lambda\mathbf{K})^{-1} \tag{5.17}$$

其中的 $\mathbf{K}$ 不依赖于 $\lambda$ （练习 5.9）。由于 $\hat{f}=\mathbf{S}_\lambda \mathbf{y}$ 为下式的解

$$\min_{\mathbf{f}} (\mathbf{y} - \mathbf{f})^T(\mathbf{y} - \mathbf{f}) + \lambda\mathbf{f}^T\mathbf{K}\mathbf{f} \tag{5.18}$$

故 $\mathbf{K}$ 也被称为 惩罚矩阵（penalty matrix），而实际上 a quadratic form in K has a representation in terms of a weighted sum of squared (divided) second differences. $\mathbf{S}_\lambda$ 的特征分解为：

$$ \mathbf{S}_\lambda = \sum_{k=1}^N \rho_k(\lambda)\mathbf{u}_k\mathbf{u}_k^T \tag{5.19}$$

其中

$$\rho_k(\lambda) = \frac{1}{1+\lambda d_k} \tag{5.20}$$

$d_k$ 为 $\mathbf{K}$ 的相应特征值。图 5.7 中上图为一些空气污染数据（128 个样本）¹的三次平滑样条结果。其中给出了两个拟合结果：一个对应着较大 $\lambda$ 的平滑拟合，和一个对应着较小 $\lambda$ 的不平滑拟合。下图中为相应的平滑器矩阵的特征值（左图）和部分特征向量²（右图）

**图 5.7**：上图为臭氧浓度对达盖特（Daggett）地区气压梯度的平滑样条曲线。两条拟合曲线分别对应了达到 5 和 7 个有效自由度的平滑参数，$\text{df}\_\lambda=\operatorname{trace}(\mathbf{S}\_\lambda)$。左下图为两个平滑样条矩阵的前 25 个特征值。前两个特征值都为 1，其他的特征值都大于等于 0。右下图为样条平滑器矩阵的第 3～6 个特征向量，其中将 $\mathbf{u}\_k$ 视作 $\mathbf{x}$ 的函数而做出曲线。图底部的地毯图标志了样本的出现位置。图中波动较弱的曲线为经过平滑处理后的函数（使用 5 个自由度的平滑器）。 — **图 5.7**：上图为臭氧浓度对达盖特（Daggett）地区气压梯度的平滑样条曲线。两条拟合曲线分别对应了达到 5 和 7 个有效自由度的平滑参数，$\text{df}_\lambda=\operatorname{trace}(\mathbf{S}_\lambda)$。左下图为两个平滑样条矩阵的前 25 个特征值。前两个特征值都为 1，其他的特征值都大于等于 0。右下图为样条平滑器矩阵的第 3～6 个特征向量，其中将 $\mathbf{u}_k$ 视作 $\mathbf{x}$ 的函数而做出曲线。图底部的地毯图标志了样本的出现位置。图中波动较弱的曲线为经过平滑处理后的函数（使用 5 个自由度的平滑器）。

以下为一些特征分解的表达方式的好处：

特征向量独立于 $\lambda$ 的取值，因此（对某特定的样本 $\mathbf{x}$ 所产生的）一组由 $\lambda$ 控制的平滑样条函数，有同样的特征向量。
$\mathbf{S}_\lambda\mathbf{y}=\sum_{k=1}^N\mathbf{u}_k\rho_k(\lambda)\langle \mathbf{u}_k,\mathbf{y}\rangle$，因此平滑样条实际上将 $\mathbf{y}$ 对（完整的）基 $\{\mathbf{u}_k\}$，再用 $\rho_k(\lambda)$ 对系数分别进行缩减。与之相比，基回归（basis-regression）方法让每一项或不做处理或缩减至 0；其诸如上文中 $\mathbf{H}_\xi$ 的投影矩阵的特征值中 $M$ 个为 1，其余为 0。因此平滑样条也被称为 收缩（shrinking）平滑器，而回归样条为 投影（projection）平滑器（见第 3.5 节：第 80 页，图 3.17）。
一系列 $\mathbf{u}_k$，对应着递减的 $\rho_k(\lambda)$，从图中看其复杂度是递增的。它们确实有高阶多项式的多次穿过 0 点的性质（zero-crossing behavior）。从 $\mathbf{S}_\lambda \mathbf{u}_k=\rho_k(\lambda)\mathbf{u}_k$ 中可见特征向量被平滑样条收缩的方式：复杂度越高，收缩的程度越大。如果 $X$ 的定义域为周期性的，则 $\mathbf{u}_k$ 为不同频率的正弦和余弦函数。
前两个特征值总为 1，它们对应着 $x$ 的线性函数特征空间的两个维度（练习 5.11），不被收缩。
特征值 $\rho_k(\lambda)=1/(1+\lambda d_k)$ 为惩罚矩阵 $\mathbf{K}$ 的特征值 $d_k$ 的倒数函数 ³，并含有系数 $\lambda$；$\lambda$ 控制着 $\rho_k(\lambda)$ 衰减至 0 的速率。同样地，$d_1=d_2=0$ 意味着线性函数不会加入惩罚。
可以用基向量 $\mathbf{u}_k$（Demmler–Reinsch 基）重写平滑样条的参数。这时平滑样条为下式的解： $$\min_\mathbf{\theta}\|\mathbf{y} - \mathbf{U}\mathbf{\theta}\|^2+ \lambda\mathbf{\theta}^T\mathbf{D}\mathbf{\theta} \tag{5.21}$$ 其中 $\mathbf{U}$ 的列为 $\mathbf{u}_k$，$\mathbf{D}$ 为对角矩阵，对角线元素为 $d_k$。
$\text{df}_\lambda=\operatorname{trace}(\mathbf{S}_\lambda)=\sum_{k=1}^N\rho_k(\lambda)$。在投影平滑器中，所有的特征值都为 1，分别对应着投影子空间的一个维度。⁴

**图 5.8**：平滑样条的平滑器矩阵呈近乎带状，相当于定义在局部的等价核函数。左图为矩阵 $\mathbf{S}$ 的热力图。右图为指定行对应的的等价核函数或加权函数。

图 5.8 描绘了一个平滑样条矩阵，在行的方向上按 $x$ 的大小顺寻排列。其中展示出的带状特性，意味着平滑样条是一个局部拟合方法，与第六章中的局部加权回归方法类似。右图为 $\mathbf{S}$ 的某几行，这个曲线也被称为 等价核函数（equivalent kernel）。随着 $\lambda\rightarrow 0$，$\text{df}_\lambda\rightarrow N$，而且 $\mathbf{S}_\lambda\rightarrow\mathbf{I}$（$N$ 维的单位矩阵）。随着 $\lambda\rightarrow\infty$，$\text{df}_\lambda\rightarrow 2$，而且 $\mathbf{S}_\lambda\rightarrow\mathbf{H}$ （对 $x$ 的线性回归帽子矩阵）。

本节练习

练习 5.7

平滑样条的推导（Green and Silverman, 1994）

假设 $N\geq2$，并且 $g$ 是点对集合 $\{x_i,z_i\}_1^N$ 的一个自然三次样条的插值函数，并有 $a<x_1<\dots<x_N<b$。这是一个将每个 $x_i$ 都作为节点的自然样条函数；这是一个 $N$ 维的函数空间，所以我们能够得到可以使样条函数的插值完全等于 $z_i$ 序列的一组系数。令 $\tilde{g}$ 代表 $N$ 个点对插值的定义在 $[a,b]$ 区间上的任意其他可导函数。

令 $h(x)=\tilde{g}̃(x)−g(x)$. 利用分部积分，以及 $g$ 是一个自然三次样条函数的已知条件，证明： $$\begin{align} \int_a^b g''(x) h''(x) dx &= - \sum_{j=1}^{N-1} g'''(x_j^+) \{ h(x_{j+1}) - h(x_j) \} \\ &= 0 \end{align}\tag{5.72}$$
所以可以证明： $$\int_a^b \tilde{g}''(t)^2 dt \geq \int_a^b g''(t)^2 dt$$ 并且只有当 $h$ 在 $[a,b]$ 区间上都等于零时，这个等号才能够成立。
考虑下面这个带惩罚的最小二乘问题： $$\min_f \left[ \sum_{i=1}^N (y_i - f(x_i))^2 + \lambda \int_a^b f''(t)^2 dt \right]$$ 利用（2）中的结论来说明它的最小值解必然是以每个 $x_i$ 为节点的一个三次样条函数。

练习 5.9

Derive the Reinsch form $\mathbf{S}_\lambda=(\mathbf{I}+\lambda\mathbf{K})^{-1}$ for the smoothing spline.

练习 5.11

Prove that for a smoothing spline the null space of $\mathbf{K}$ is spanned by functions linear in $\mathbf{X}$.

原文为“Daggot pressure gradient”，但译者无法搜索到“Daggot”，但经过关联找到了“Daggett”，为加利福尼亚州的一个气象站所在地。（链接） ↩︎
$\mathbf{u}_k$ 和 $\mathbf{x}$ 均是长度为 $N$ 的向量。 ↩︎
原文为“… an inverse function of …”，译者认为应理解为“倒数”，而不是“反函数”。 ↩︎
投影平滑器（基回归）的特征值为 1 或 0，投影子空间即为所有被选入模型的输入特征变量（特征值为 1）生成的空间。 ↩︎