4.3. 正交匹配追踪¶
4.3.1. 什么是正交匹配追踪¶
正交匹配追踪 (Orthogonal Matching Pursuit , OMP) 是一种信号稀疏分解方法, 由 Pati 等人 [10] 于1993年提出, OMP被广泛用于稀疏信号恢复 [11]. 稀疏分解旨在将信号分解为过完备字典中的尽可能少(稀疏)的原子的线性组合, 从而获得信号的简洁有效的表示.
4.3.2. 正交匹配追踪算法¶
正交匹配追踪算法步骤见算法1 [11].
注解
算法1: 正交匹配追踪算法步骤
输入: 列归一化过完备字典矩阵 \({\bm A} \in {\mathbb R}^{m \times n}\), 观测向量 \({\bm y}\in {\mathbb R}^{m \times 1}\), 待恢复信号稀疏度 \(K\), 容忍残差模值 \(\epsilon\)
输出: 待估计稀疏信号 \(\hat{\bm x} \in {\mathbb R}^{n\times 1}\), 索引集合 \({\mathbb I} \subset {\mathbb U} = \{1,2,\cdots,n\}\), 观测向量近似 \(\hat{\bm y}\) 和残差向量 \({\bm r} = {\bm y} - \hat{\bm y}\)
Step1: 初始化稀疏信号 \({\bm x}_{0} = {\bm 0}\), 残差 \({\bm r}_{0} = {\bm y} - \bm{A}{\bm x}_0 = {\bm y}\), 索引集合(支撑) \({\mathbb I}_{0} = \empty\), 迭代计数器 \(k = 1\).
Step2: 选择与残差最相关的原子, 即求解索引 \(i_k\) 满足优化问题
Step3: 更新索引集 \({\mathbb I}_{k} = {\mathbb I}_{k-1} \cup {i_k}\).
Step4: 求解新的最小二乘优化问题
\[{\bm x}_k = \arg \min _{\bm x} \|{\bm y} - {\bm A}_{{\mathbb I}_k}{\bm x}\|_2^2 \]求得的稀疏系数解为 \({\bm x}_k = ({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k})^{-1}{\bm A}_{{\mathbb I}_k}^T{\bm y}\)
Step5: 计算新的投影近似, 残差
\[\begin{aligned} \hat{\bm y}_{k} &=\bm{A}_{{\mathbb I}_k} \bm{x}_{k} \\ \bm{r}_{k} &=\bm{y}-\hat{\bm y}_{k} \end{aligned}. \]
Step6: 更新迭代计数器 \(k = k + 1\) 若 \(k<K\) 且 \(||{\bm r}||_2 < \epsilon\), 重复 Step2 至 Step5, 否则停止迭代, 转 Step7.
Step7: 输出 \(\hat{\bm x} = \bm{x}_{k}\), \(\hat{\bm y} = \hat{\bm y}_{k}\), \({\bm r} = {\bm r}_k\). 其中, \(\hat{\bm x}\) 在位置 \({\mathbb I}_k\) 处非零.
警告
待查证是否有相关文献, 没有的话, 是否可以发表 实际迭代过程中, 在求解逆矩阵 \(({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k})^{-1}\) 时容易出现其奇异的情况, 可以通过求解 \(({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k} + \alpha {\bm I})^{-1}\) 来代替, 其中 \(\alpha > 0\).
提示
记 \({\bm P}_k = {\bm A}_{{\mathbb I}_k}({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k})^{-1}{\bm A}_{{\mathbb I}_k}^T\), 易知 \({\bm P}_k\) 为正交投影矩阵(线性投影), 将向量投影到由 \({\bm A}_{{\mathbb I}_k}\) 的元素(类比坐标系)张成的线性空间中, 从而 \({\hat{\bm y}}_k = {\bm A}_{{\mathbb I}_k} {\bm x}_k = {\bm P}_k {\bm y}\), 则投影后的误差可以表示为
其中, \(({\bm I} - {\bm P}_k){\bm {Ax}}\) 为残差中的信号成份, \(({\bm I} - {\bm P}_k){\bm n}\) 为残差中的噪声成份.
4.3.3. 实验与分析¶
主要分析OMP算法在信号稀疏分解上的性能, 有关OMP在压缩感知中的应用参见 稀疏信号实验 小节.
仿真数据实验¶
实验中通过仿真生成含有三个频率成分的信号 \({\bm y} = {\rm sin}2πf_1t + {\rm sin}2πf_2t +{\rm sin}2πf_3t\), 仿真参数设置如下:
频率成分: \(f_1 = 10Hz, f_2 = 20Hz,f_3 = 70Hz\)
采样频率: \(F_s = 256Hz\)
采样时间: \(T_s = 1s\)
采样点数: \(N_s = F_s T_s = 256\)
构造过完备DCT字典(参见 过完备DCT字典), 字典大小为 \(M\times N\), 其中 \(M=N_s, N=R*M\), \(R\) 为一可调变量, 通过OMP求解稀疏系数 \({\bm x}\), 并利用 式.4.1 恢复信号 \({\bm y}\).
实现代码, 参见文件 demo_omp_sin3_decomposition.py .
仿真生成的信号 \({\bm y}\) 如下图(左)所示, 其DCT变换后的系数如下图(右)所示
图 4.20 Signal in Time (left, \({\bm y}\) ) and Frequency (right, \({\bm x}\)) domains.¶
如图中所示, 三种频率成分对应六个峰, 峰的位置为频率的二倍(这是因为DCT变换中分母为 \(2N\), 而不是 \(N\)), 峰的上下两个方向对应正频和负频.
\(R=2\) 即 \(N = 2M\) 时, 设置不同稀疏度下的信号重构均方误差结果如下:
---MSE(y, y1) with k = 2: 0.937048086151075
---MSE(y, y2) with k = 4: 0.5420001800721357
---MSE(y, y3) with k = 6: 0.25283463530586237
---MSE(y, y4) with k = 100: 0.002307660406098402
求解的频域稀疏系数为
图 4.21 Recovered signal in frequency domain with different sparse degree \(k\).¶
恢复的信号为
图 4.22 Recovered signal in time domain with different sparse degree \(k\).¶
\(R=1\) 即 \(N = M\) 时, 设置不同稀疏度下的信号重构均方误差结果如下:
---MSE(y, y1) with k = 2: 0.9865788586325661
---MSE(y, y2) with k = 4: 0.5829159447924545
---MSE(y, y3) with k = 6: 0.3169739132218109
---MSE(y, y4) with k = 100: 0.0023321884115169826
求解的频域稀疏系数为
图 4.23 Recovered signal in frequency domain with different sparse degree \(k\).¶
恢复的信号为
图 4.24 Recovered signal in time domain with different sparse degree \(k\).¶
\(R=0.5\) 即 \(N = M/2\) 时, 设置不同稀疏度下的信号重构均方误差结果如下:
---MSE(y, y1) with k = 2: 1.031848207375802
---MSE(y, y2) with k = 4: 0.6615399428007099
---MSE(y, y3) with k = 6: 0.4982313612363407
---MSE(y, y4) with k = 100: 0.32437796224425747
求解的频域稀疏系数为
图 4.25 Recovered signal in frequency domain with different sparse degree \(k\).¶
恢复的信号为
图 4.26 Recovered signal in time domain with different sparse degree \(k\).¶