4.3. 正交匹配追踪

4.3.1. 什么是正交匹配追踪

正交匹配追踪 (Orthogonal Matching Pursuit , OMP) 是一种信号稀疏分解方法, 由 Pati 等人 [10] 于1993年提出, OMP被广泛用于稀疏信号恢复 [11]. 稀疏分解旨在将信号分解为过完备字典中的尽可能少(稀疏)的原子的线性组合, 从而获得信号的简洁有效的表示.

good

4.3.2. 正交匹配追踪算法

正交匹配追踪算法步骤见算法1 [11].

注解

算法1: 正交匹配追踪算法步骤

输入: 列归一化过完备字典矩阵 ARm×n{\bm A} \in {\mathbb R}^{m \times n}, 观测向量 yRm×1{\bm y}\in {\mathbb R}^{m \times 1}, 待恢复信号稀疏度 KK, 容忍残差模值 ϵ\epsilon

输出: 待估计稀疏信号 x^Rn×1\hat{\bm x} \in {\mathbb R}^{n\times 1}, 索引集合 IU={1,2,,n}{\mathbb I} \subset {\mathbb U} = \{1,2,\cdots,n\}, 观测向量近似 y^\hat{\bm y} 和残差向量 r=yy^{\bm r} = {\bm y} - \hat{\bm y}

Step1: 初始化稀疏信号 x0=0{\bm x}_{0} = {\bm 0}, 残差 r0=yAx0=y{\bm r}_{0} = {\bm y} - \bm{A}{\bm x}_0 = {\bm y}, 索引集合(支撑) I0={\mathbb I}_{0} = \empty, 迭代计数器 k=1k = 1.

Step2: 选择与残差最相关的原子, 即求解索引 iki_k 满足优化问题

ik=argmaxiUAiTrk1.i_k = \arg \max _{i\in{\mathbb U}} |{\bm A}_i^T{\bm r}_{k-1}|.

Step3: 更新索引集 Ik=Ik1ik{\mathbb I}_{k} = {\mathbb I}_{k-1} \cup {i_k}.

Step4: 求解新的最小二乘优化问题

xk=argminxyAIkx22{\bm x}_k = \arg \min _{\bm x} \|{\bm y} - {\bm A}_{{\mathbb I}_k}{\bm x}\|_2^2

求得的稀疏系数解为 xk=(AIkTAIk)1AIkTy{\bm x}_k = ({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k})^{-1}{\bm A}_{{\mathbb I}_k}^T{\bm y}

Step5: 计算新的投影近似, 残差

y^k=AIkxkrk=yy^k.\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+1k = k + 1k<Kk<Kr2<ϵ||{\bm r}||_2 < \epsilon, 重复 Step2 至 Step5, 否则停止迭代, 转 Step7.

Step7: 输出 x^=xk\hat{\bm x} = \bm{x}_{k}, y^=y^k\hat{\bm y} = \hat{\bm y}_{k}, r=rk{\bm r} = {\bm r}_k. 其中, x^\hat{\bm x} 在位置 Ik{\mathbb I}_k 处非零.

警告

待查证是否有相关文献, 没有的话, 是否可以发表 实际迭代过程中, 在求解逆矩阵 (AIkTAIk)1({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k})^{-1} 时容易出现其奇异的情况, 可以通过求解 (AIkTAIk+αI)1({\bm A}_{{\mathbb I}_k}^T{\bm A}_{{\mathbb I}_k} + \alpha {\bm I})^{-1} 来代替, 其中 α>0\alpha > 0.

提示

Pk=AIk(AIkTAIk)1AIkT{\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, 易知 Pk{\bm P}_k 为正交投影矩阵(线性投影), 将向量投影到由 AIk{\bm A}_{{\mathbb I}_k} 的元素(类比坐标系)张成的线性空间中, 从而 y^k=AIkxk=Pky{\hat{\bm y}}_k = {\bm A}_{{\mathbb I}_k} {\bm x}_k = {\bm P}_k {\bm y}, 则投影后的误差可以表示为

rk=yy^k=(IPk)y=(IPk)Ax+(IPk)n,\begin{aligned} {\bm r}_k &= {\bm y} - {\hat{\bm y}}_k \\ &= ({\bm I} - {\bm P}_k){\bm y} \\ &= ({\bm I} - {\bm P}_k){\bm {Ax}} + ({\bm I} - {\bm P}_k){\bm n}, \end{aligned}

其中, (IPk)Ax({\bm I} - {\bm P}_k){\bm {Ax}} 为残差中的信号成份, (IPk)n({\bm I} - {\bm P}_k){\bm n} 为残差中的噪声成份.

4.3.3. 实验与分析

主要分析OMP算法在信号稀疏分解上的性能, 有关OMP在压缩感知中的应用参见 稀疏信号实验 小节.

仿真数据实验

实验中通过仿真生成含有三个频率成分的信号 y=sin2πf1t+sin2πf2t+sin2πf3t{\bm y} = {\rm sin}2πf_1t + {\rm sin}2πf_2t +{\rm sin}2πf_3t, 仿真参数设置如下:

  • 频率成分: f1=10Hz,f2=20Hz,f3=70Hzf_1 = 10Hz, f_2 = 20Hz,f_3 = 70Hz

  • 采样频率: Fs=256HzF_s = 256Hz

  • 采样时间: Ts=1sT_s = 1s

  • 采样点数: Ns=FsTs=256N_s = F_s T_s = 256

构造过完备DCT字典(参见 过完备DCT字典), 字典大小为 M×NM\times N, 其中 M=Ns,N=RMM=N_s, N=R*M, RR 为一可调变量, 通过OMP求解稀疏系数 x{\bm x}, 并利用 式.4.1 恢复信号 y{\bm y}.

实现代码, 参见文件 demo_omp_sin3_decomposition.py .

仿真生成的信号 y{\bm y} 如下图(左)所示, 其DCT变换后的系数如下图(右)所示

Signal in Time and Frequency domains.

图 4.20 Signal in Time (left, y{\bm y} ) and Frequency (right, x{\bm x}) domains.

如图中所示, 三种频率成分对应六个峰, 峰的位置为频率的二倍(这是因为DCT变换中分母为 2N2N, 而不是 NN), 峰的上下两个方向对应正频和负频.

R=2R=2N=2MN = 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

求解的频域稀疏系数为

Recovered signal in frequency domain with different sparse degree :math:`k`.

图 4.21 Recovered signal in frequency domain with different sparse degree kk.

恢复的信号为

Recovered signal in time domain with different sparse degree :math:`k`.

图 4.22 Recovered signal in time domain with different sparse degree kk.

R=1R=1N=MN = 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

求解的频域稀疏系数为

Recovered signal in frequency domain with different sparse degree :math:`k`.

图 4.23 Recovered signal in frequency domain with different sparse degree kk.

恢复的信号为

Recovered signal in time domain with different sparse degree :math:`k`.

图 4.24 Recovered signal in time domain with different sparse degree kk.

R=0.5R=0.5N=M/2N = 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

求解的频域稀疏系数为

Recovered signal in frequency domain with different sparse degree :math:`k`.

图 4.25 Recovered signal in frequency domain with different sparse degree kk.

恢复的信号为

Recovered signal in time domain with different sparse degree :math:`k`.

图 4.26 Recovered signal in time domain with different sparse degree kk.

真实数据实验