5.9 小波平滑

目前介绍了对基函数字典的两种不同处理方式。在回归样条中，模型选择基函数的一个子集，这个选择基于对问题的主观认知或某种自动方法。自适应性更强的方法，比如第九章中的 MARS，可以同时捕捉到平滑和不平滑的数据关系。在平滑样条中，模型使用全部基函数，但会将系数向更平滑的曲线方向收缩。

小波平滑方法通常用完整的标准正交基函数来表达函数，再向一个稀疏的函数表达式对系数进行选择和收缩。就如同一个平滑函数可表达成几个样条基函数的组合，一个大部分为水平但有几个孤立的凸起的函数可表达成几个（含有凸起的）基函数的组合。小波基函数在信号处理和压缩中非常流行，因为它们能高效率地表述平滑并且/或者局部凸起的函数，这种函数也称为 时频局部化（time and frequency localization）。与之相比，传统的傅立叶基函数只允许频率局部化。

在进入具体细节之前，先通过图 5.16 中左边的 哈尔（Harr） 小波函数 来直观地理解小波平滑的原理。纵轴表示小波函数的尺度（频率），在底部取值低，在顶部取值高。图中只展示了几个尺度；为了方便展示，每个尺度下的小波函数被放置到不同的时间横轴位置上。小波平滑用最小二乘拟合这些基函数的系数，然后舍弃（过滤）其中比较小的系数。由于每个尺度下都有很多基函数，所以拟合结果会纳入其所需的小波基函数，舍弃不需要的基函数，从而达到时频局部化。哈尔小波函数易于理解，但在大多场景中不够平滑。图 5.16 中右边的 Symlet¹ 小波函数则更平滑，也同时具有标准正交性。

图 5.17 为核磁共振（NMR, nuclear magnetic resonance）信号，看上去其中包含了一些平滑的成分和孤立的尖峰，以及一些噪声。下左图为对其利用 symlet 基函数的小波转换。小波函数的系数按行排列，底部的尺度低，顶部的尺度高。每个（竖）线段的长度代表着系数的大小。下右图为阈值筛选后的小波函数系数。下面的式 5.69 给出了阈值筛选过程，其与线性回归的套索算法中的“软阈值”规则是一样的（第 3.4.2 节）。注意到其中很多小系数被设置为零。上图中的绿色曲线为用阈值筛选过的系数还原的信号，这是对原信号的平滑处理。下一节会详细介绍这个构建小波函数和阈值筛选规则的过程。

5.9.1 小波函数基底和小波转换 😱

本节详细介绍小波函数的构建和过滤。 对一个尺度函数 $\phi(x)$ （也被称为“父”小波）的平移和伸缩，就生成了小波基底。

图 5.18 中红色曲线为哈尔（Haar）和 symlet-8 的尺度函数。对于了解方差分析和树模型的读者，哈尔基函数非常易于理解；它其实是一个分段常数函数。因此，若 $\phi(x)=I(x\in[0,1])$，则 $\phi_{0,k}(x)=\phi(x-k)$（$k$ 为整数）即为在整数处跳动的函数的空间上的一组标准正交基。将此记为 参考（reference） 空间 $V_0$。伸缩后的函数 $\phi_{1,k}(x)=\sqrt{2}\phi(2x-k)$，即为在长度 $\frac{1}{2}$ 的区间上分段常数的函数空间的一组标准正交基，其生成的函数空间 $V_1\supset V_0$。实际上，可广泛地定义 $\cdots\supset V_1\supset V_0\supset V_{-1}\supset\cdots$，其中每个 $V_j$ 为 $\phi_{j,k}=2^{j/2}\phi(2^jx-k)$ 生成的空间。

现在回到小波函数的定义。在方差分析中常常将两个均值 $\mu_1$ 和 $\mu_2$ 表达为它们的总均值 $\mu=\frac{1}{2}(\mu_1+\mu_2)$，和对比度（contrast）$\alpha=\frac{1}{2}(\mu_1-\mu_2)$。当对比度 $\alpha$ 很小时，可以将其简化地设置为零。于此类似，空间 $V_{j+1}$ 上的一个函数，可表达为空间 $V_j$ 中的一个元素与一个 $V_j$ 在 $V_{j+1}$ 中的正交补集 $W_j$ 中一个元素的组合，可写为 $V_{j+1}=V_j\bigoplus W_j$。空间 $W_j$ 中的元素为 细节（detail），其中的一部分可能会被设为零。可证明由 母小波（mother wavelet） $\psi(x)=\phi(2x)-\phi(2x-1)$ 生成的函数 $\psi(x-k)$ 形成了哈尔函数的 $W_0$ 空间的一组标准正交基。同样地，$\psi_{j,k}=2^{j/2}\psi(2^jx-k)$ 为 $W_j$ 的一组标准正交基。

可见 $V_{j+1}=V_j\bigoplus W_j=V_{j-1}\bigoplus W_{j-1}\bigoplus W_j$，故一个函数可被表达为其 $j$ 级的细节（detatil）和 $j$ 级的粗略（rough）成分，而后者又可被拆分为 $(j-1)$ 级的细节和 $(j-1)$ 级的粗略成分，以此类推。最终，可得到表达式 $V_J=V_0\bigoplus W_0\bigoplus W_1\cdots\bigoplus W_{J-1}$。图 5.16 展示了几个特定的小波函数 $\psi_{j,k}(x)$。

由于这些空间互相正交，所有的基函数都是标准正交。实际上，若函数的定义域为离散的 $N=2^J$ 个点（时间），上述的分解已经到了极限。在 $j$ 级，共有 $2^j$ 个基函数元素，加总后 $W_j$ 中共有 $2^J-1$ 个元素，而 $V_0$ 中有一个元素。这种结构的标准正交基对应着下面会接续介绍的 多解析度分析（multiresolution analysis, MRA）。

哈尔（Haar）基函数虽然可帮助理解上述的构建方式，但在实际场景中过于粗糙。幸好还有其他设计巧妙的小波基。图 5.16 和 5.18 中包含了 多贝西（Daubechies）symlet-8 基函数。与哈尔（Haar）基函数相比，这个基函数更平滑，但也存在取舍：

• 每个小波函数的定义域覆盖着 15 个连续的时间区间，而哈尔基函数只有一个。一般地，symmlet-p 小波函数的定义域为 $2p-1$ 个连续的区间。定义域越宽，为了可得到平滑的曲线，小波函数在更多的时间点上要降至零。注意最后的有效定义域可能会更狭小。
• symlet-p 小波函数 $\psi(x)$ 有 $p$ 个零矩，即： $$\int\psi(x)x^jdx=0,j=0,\dots,p-1$$ 这意味着在 $N=2^J$ 个时间点上任意 $p$ 阶多项式均存在于 $V_0$ 空间（练习 5.18）。也就是说，$V_0$ 与平滑样条惩罚项的零空间是等价的。哈尔（Haar）小波函数有一个零矩，$V_0$ 包含了所有常数函数。

symlet-p 的尺度函数是很多小波函数生成函数族中的一个。对哈尔（Haar）基函数也有类似的操作：

• 若 $V_0$ 为 $\phi(x-k)$ 生成的空间，则 $V_1\supset V_0$ 为 $\phi_{1,k}(x)=\sqrt{2}\phi(2x-k)$ 生成的空间；并且 $\phi(x)=\sum_{k\in\mathcal{Z}}h(k)\phi_{1,k}(x)$，其中 $h(k)$ 为某种滤波系数。
• $W_0$ 为 $\psi(x)=\sum_{k\in\mathcal{Z}}g(k)\phi_{1,k}(x)$，其中滤波系数为 $g(k)=(-1)^{1-k}h(1-k)$。

5.9.2 自适应小波滤波 😱

小波函数非常适合用于均匀栅格上测量的数据，例如离散化的信号、图像、或时间序列。方便起见，以下考虑一维的在 $N=2^J$ 个格栅点上的数据。假设 $\mathbf{y}$ 为输出向量，$\mathbf{W}$ 为 $N\times N$ 的标准正交小波基函数向量，其取值在 $N$ 个均匀分布的观测点上。则 $\mathbf{y}^*=\mathbb{W}^T\mathbf{y}$ 被称为 $\mathbf{y}$ 的小波转换（同时也是完整最小二程回归系数）。一个常用的自适应小波拟合方法被称为 SURE（Stein Unbiased Risk Estimation） 收缩（Donoho and Johnstone (1994)）。其最小化准则为：

这与第三章中的套索回归准则一样。由于 $\mathbf{W}$ 为标准正交，可得到简单的解：

最小二乘的系数被向零点移动，同时在零点截断。拟合函数（向量）可通过拟小波转换得到 $\hat{\mathbf{f}}=\mathbf{W} \hat{\mathbf{\theta}}$。

一个简单的 $\lambda$ 选择为 $\lambda=\sigma\sqrt{2\log N}$，其中 $\sigma$ 为噪声项标准差的估计。这个选择的原因如下。由于 $\mathbf{W}$ 为标准正交转换，若 $\mathbf{y}$ 中的元素为白噪声（独立的高斯分布变量，均值为 0 方差为 $\sigma^2$），则 $\mathbf{y}^*$ 也同样为白噪声。已知当随机变量 $Z_1$，$Z_2$，……，$Z_N$ 为白噪声时，$|Z_j|,j=1,\dots,N$ 的最大值的期望近似为 $\sigma\sqrt{2\log N}$。因此，所有水平低于 $\sigma\sqrt{2\log N}$ 的系数很可能是噪声，可被设置为零。

空间 $\mathbf{W}$ 可以为任意的标准正交基函数：多项式、自然样条、或余弦曲线。小波函数的特殊之处在于所用基函数的特定形式，使函数表达式可在时间和频率上局部化。

在图 5.17 的核磁共振信号中，使用了 symlet-8 基函数计算小波转换。注意到系数并没有下降到 $V_0$ 上，而是停止在了含有 16 个基函数的 $V_4$。随着细节每提高一个层级，系数随之变小，只有在函数的凸起尖峰位置才会有留存的系数。小波系数体现了信号在时间上（不同层级的基函数为彼此的平移）和频率上的局部特征。而每一次收缩都以两倍因子增加细节，可与传统的傅立叶表达式中的两倍频率相对应。实际上从数学的角度理解，某个特定尺度的小波函数对应着一个限制在有限区域或八开（octave）频率上的傅立叶转换。

图右侧的收缩/截断是通过之前提到过的 SURE 方法得到的。标准正交的 $N\times N$ 基函数矩阵 $\mathbf{W}$，其每列为小波基函数在 $N$ 个时间点上的取值。在这个图中，矩阵有 16 个列，对应着 $\phi_{4,k}(x)$，残差（remainder）对应着 $\psi_{j,k}(x)$，$j=4,\dots,11$。在实践中，$\lambda$ 取决于噪声项的方差，需要从数据中估计得出（比如在最高等级的系数的方差）。

注意 SURE 最小化准则（5.68）与平滑样条准则（5.21）的相似之处：

• 两者均为从粗糙到精细细节的层级结构，但小波函数在每个解析度（resolution）等级内同时在时间上局部化。
• 样条函数通过差异化的收缩常数 $d_k$ 达到偏向于平滑函数的结果。早期版本的 SURE 收缩在所有的尺度下相同。S+wavelets 函数 waveshrink() 有很多选项，其中有些允许差异化的收缩。
• 样条中的 $L_2$ 惩罚项只产生收缩，而 SURE 中的 $L_1$ 惩罚项同时起收缩和变量选择的作用。

更一般地理解，平滑样条通过要求平滑程度来压缩原始信号，而小波函数是通过要求稀疏程度。图 5.19 中在两组性质不同的样本中，对比了小波函数拟合（使用 SURE 收缩）和平滑样条拟合（使用交叉验证）。在上图的核磁共振数据中，平滑样条为了捕捉到孤立的尖峰，需要在全局上进行调整；小波函数则很好地在局部拟合尖峰。在下图中，真实的关系函数是平滑的，噪声相对更高。小波函数拟合引入了多余不必要的摆动，这是其为更好适应性而在方差上付出的代价。

实际计算中，小波转换并不是通过矩阵乘法 $\mathbf{y}^*=\mathbf{W}^T\mathbf{y}$ 得出的。利用合理的金字塔状机制，$y^*$ 的计算成本为 $O(N)$，这甚至要快于快速傅立叶转换（fast Fourier transform, FFT）的 $N\log(N)$。普遍的算法构建方法超出了本书的范围，但对哈尔（Haar）基函数的构建比较简单（练习 5.19）。同样地，拟小波转换 $\mathbf{W}\hat{\mathbf{\theta}}$ 的计算量也是 $O(N)$。

这个领域的内容很多发展很快，本节只是一个简要的概览。在数学和计算机领域都有建立在小波函数上庞大的理论基础。现代的图像压缩常常是通过二维的小波表达式实现的。

本节练习

练习 5.18

The wavelet function ψ(x) of the symmlet-p wavelet basis has vanishing moments up to order p. Show that this implies that polynomials of order p are represented exactly in V0 , defined on page 176.

练习 5.19

Show that the Haar wavelet transform of a signal of length N = 2J can be computed in O(N ) computations.