5.a 样条函数的计算 😱

本附录小节将介绍用于表达多项式样条的 B 样条基函数，以及它们在平滑样条的计算中的应用。

5.10.1 B 样条

在开始之前，我们需要先对第 5.2 节中定义的节点（knot）序列做一个补充。令 ${\xi }_0<{\xi }_1$ 和 ${\xi }_K<{\xi }_{K+1}$ 为两个边界（boundary）节点，它们一般会定义出我们想要计算样条模型取值的定义域。现将扩大节点序列 $\tau$ 定义为：

• ${\tau }_1\le {\tau }_2\le \cdots \le {\tau }_M\le {\xi }_0$
• ${\tau }_{j+M}={\xi }_j,j=1,\dots ,K$
• ${\xi }_{K+1}\le {\tau }_{K+M+1}\le {\tau }_{K+M+2}\le \cdots \le {\tau }_{K+2M}$

这些超出边界的额外节点的实际取值是任意的，一般会让它们都相等并分别等于 ${\xi }_0$ 和 ${\xi }_{K+1}$ 上的取值。

将节点序列 $\tau$ 的第 $i$ 个 $m$ 阶的 B 样条基函数记为 $B_{i,m}(x)$，$m\le M$。 它们是以如下的差商形式递归地定义出来的。

以上这些也被称作哈尔（Haar）基函数。

所以当 $M=4$ 时，$B_{i,4}$（$i=1,\dots ,K+1$）是节点序列 $\xi$ 的 $K+4$ 个三次 B 样条基函数。继续这样递归下去，就会得出任意阶样条的 B 样条基函数。

图 5.20 展示了节点为 $\{0.0,0.1,\dots ,1.0\}$ 的一阶到四阶的 B 样条序列。由于生成了一些重复的节点，需要做一些处理来避免分母为零的情况。如果按惯例当 ${\tau }_i={\tau }_{i+1}$ 时令 $B_{i,1}=0$，那么根据归纳法当 ${\tau }_i={\tau }_{i+1}=\cdots ={\tau }_{i+m}$ 时 $B_{i,m}=0$。同时也需要注意到在上述的基函数构造中，节点 $\xi$ 的 $m<M$ 阶 B 样条基函数只需要 $B_{i,m},i=M-m+1,\dots ,M+K$ 这个子集。

要完全理解这些函数的性质，并证明它们的确张成了节点序列的三次样条空间，这需要额外的数学工具，包括了差商的一些性质。练习 5.2 给出了这方面的一些内容。

B 样条的适用范围实际上比这里介绍的要更宽泛，它的特点就在于节点复制。如果在上文 $\tau$ 序列的构建中复制了一个内部节点，然后在和之前一样生成 B 样条序列，那么得出的基函数张成了在被复制点处缺少一个连续导数的分段多项式空间。一般来说，如果除重复的边界节点之外也将内部节点 ${\xi }_j$ 重复了 $1\le r_j\le M$ 次，那么在 $x={\xi }_j$ 处的最低阶不连续导数会是 $M-r_j$ 阶。因此在没有重复的三次样条中，$r_j=1$，$j=1,\dots ,K$，并且在每个内部节点处第三阶导数（$4-1$）是不连续的。将节点 $j$ 重复三次将会得到一个不连续的一阶导数；将节点重复四次将会得到一个不练习的“零阶”导数，也就是函数本身在 $x={\xi }_j$ 处不连续。这正是在边节点处所出现的情况；将（边界）节点重复 $M$ 次，所以样条函数在边界节点处变得不连续（即在边界之外无法定义）。

B 样条的局部支撑集（support）在计算上有重要的意义，特别是当节点数量 $K$ 比较大时。有 $N$ 个观测样本和 $K+M$ 个特征变量（基函数）的最小二乘计算会消耗 $O(N(K+M{)}^2+(K+M{)}^3)$ 次浮点运算（flop）。若 $K$ 与 $N$ 成一个可观的比例，那么算法的消耗将是 $O(N^3)$，当 $N$ 比较大时这个计算量会变得不可行。如果将 $N$ 个观测样本排序，$N\times (K+M)$ 维度的回归矩阵包含了 $K+M$ 个 B 样条基函数在 $N$ 个点上的取值，其中会有很多的零，那么可利用这一点来将计算复杂度降低回 $O(N)$。下一节将对此展开介绍。

5.10.2 平滑样条的计算

虽然自然样条（第 5.2.1 节）给出了平滑样条的一组基函数，但是从计算量方面来看，在无约束的 B 样条的更大的空间上处理起来更方便。记 $f(x)=\sum _1^{N+4}{\gamma }_j{B}_j(x)$，其中 ${\gamma }_j$ 是系数而 $B_j$ 是三次 B 样条基函数。它的解和之前看起来一样：

只不过用 $N\times (N+4)$ 的矩阵 $\mathbf{B}$ 代替了之前 $N\times N$ 的矩阵 $\mathbf{N}$，以及类似地用 $(N+4)\times (N+4)$ 的惩罚矩阵 ${\mathbf{\Omega }}_B$ 代替了之前的 $N\times N$ 的矩阵 ${\mathbf{\Omega }}_N$。尽管从表面上看起来似乎并没有边界导数约束条件，但实际上惩罚项会给边界外任意非零导数赋予相当于无穷大的权重，从而自动地实施了那些约束条件。 在实际应用中，会将 $\hat{\gamma }$ 约束在一个使惩罚项总是有限的子空间上。

由于 $\mathbf{B}$ 的列是在从小到大排序的 $X$ 取值处的 B 样条函数取值，而且三次 B 样条有局部支撑集，所以 $\mathbf{B}$ 是一个带状（banded）矩阵（最下侧 4 个对角为零）。因而 $\mathbf{M}=({\mathbf{B}}^T\mathbf{B}+\lambda \mathbf{\Omega })$ 也是一个带状矩阵，并且可以容易地计算出它的 Cholesky 分解 $\mathbf{M}=\mathbf{L}{\mathbf{L}}^T$。然后就可通过回代（back-substitution）求解 $\mathbf{L}{\mathbf{L}}^T\mathbf{gamma}={\mathbf{B}}^T\mathbf{y}$，从而以 $O(N)$ 次运算得出 $\mathbf{f}$ 以及预测解 $\hat{f}$。

在实践中，当 $N$ 比较大时并不一定要使用所有的 $N$ 个内部节点，任意合理的“薄化（thinning）”策略都可能会降低计算量并对拟合有微乎其微的影响。例如 S-PLUS 中的 smooth.spline 函数就使用了一个近似对数的策略：当 $N<50$ 时使用所有的节点，但当 $N=5000$ 时只使用 204 个节点。

本节练习

练习 5.2

假设 $B_{i,M}(x)$ 是 186 页中式 5.77 至 5.78 中所定义的 $M$ 阶 B 样条。

1. 用归纳法证明当 $x\notin [{\tau }_i,{\tau }_{i+M}]$ 时 $B_{i,M}=0$。这就证明了例如三次 B 样条的支撑集最多是 5 个节点。
2. 用归纳法证明当 $x\in ({\tau }_i,{\tau }_{i+M})$ 时 $B_{i,M(x)>0}$。即 B 样条在它们的支撑集内部是正数。
3. 用归纳法证明 $\sum _{i=1}^{K+M}{B}_{i,M}(x)=1,\mathrm{\forall }x\in [{\xi }_0,{\xi }_{K+1}]$。
4. 证明 $B_{i,M}$ 是一个在 $[{\xi }_0,{\xi }_{K+1}]$ 区间上的 $M$ 阶（$M-1$ 次）的分段多项式函数，而且只在节点 ${\xi }_1,\dots ,{\xi }_K$ 处有间断。
5. 证明一个 $M$ 阶的 B 样条基函数是 $M$ 个均匀分布随机变量的卷积（convolution）的密度函数。