正交匹配追踪

这篇博文是在对Koredianto Usman《Introduction to Orthogonal Matching Pursuit》文章的翻译，后面附带了一些总结.

这篇文章是前面《Matching Pursuit (MP)》文章的继续. （注：原文中还是有一些小细节错误，请大家睁大眼睛阅读）

简介

考虑下面的问题：给定x=[−1.210]和A=[−0.7070.7070.80.60−1]，计算y=A⋅x.

这相当简单！

y = A \cdot x = [- 0.707 0.707 0.8 0.6 0 - 1] \cdot [- 1.2 10] = [1.65 - 0.25]

现在是比较困难的部分！给定y=[1.65−0.25]和A=[−0.7070.7070.80.60−1]，如何找到原始的x（或者尽可能接近原始的x）？

在压缩感知（Compressive Sensing）术语中，从x和A中得到y叫做压缩，反之，从y和A得到x叫重建（个人感觉原文中好像说反了，还是我理解错了？）. 重建问题不是一个很简单的问题。

一般地，x被称为original signal或者original vector，A被称为compression matrix或者sensing matrix，y称为compressed signal或者compressed vector.

基本概念

和前面在MP算法中讨论过的一样，感知矩阵A可以看做列向量的集合：

A = [- 0.707 0.707 0.8 0.6 0 - 1] = [b 1 b 2 b 3]

其中，

b 1 = [- 0.707 0.707], b 2 = [- 0.8 0.6], b 1 = [0 - 1]

这些列被称为基（basis）或者原子（atom）.

现在，我们令x=⎡⎣⎢x1x2x3⎤⎦⎥，则

A \cdot x = [- 0.707 0.707] \cdot x 1 + [0.8 0.6] \cdot x 2 + [0 - 1] \cdot x 3 = y

从该方程可以看出y是使用x中的稀疏对原子的线性组合. 实际上，我们知道x1=−1.2，x2=1，x3=0. 换言之，对于得到y值，b1的贡献值为-1.2，b2为1，而b2为0.

OMP算法和MP算法类似，都是从字典中找出哪一个原子对y值的贡献最大，接下来是哪个原子的贡献值大，以此类推. 我们现在知道这个过程需要N次迭代，N是字典中原子的个数. 在该例子中，N等于3.

计算原子的贡献值

字典中有三个原子，

b 1 = [- 0.707 0.707], b 2 = [- 0.8 0.6], b 1 = [0 - 1]

且

y = [1.65 - 0.25]

因此，各原子对y的贡献值为

< b 1, y > = [- 0.707 0.707] T \cdot [1.65 - 0.25] = - 0.707 \cdot - 1.65 + 0.707 \cdot - 0.25 = - 1.34

< b 2, y > = 1.17

<b3,y>=0.25

显然，原子b1对y值的贡献最大为-1.34（按绝对值计算）.

在第一次迭代过程中，我们选择b1作为基，基的系数为-1.34.

当然，我们也可以使用矩阵的形式计算点积：

w = A T \cdot y = ⎡ ⎣ ⎢ - 0.707 0.8 0 0.707 0.6 - 1 ⎤ ⎦ ⎥ \cdot [1.65 - 0.25] = ⎡ ⎣ ⎢ - 1.34 1.17 0.25 ⎤ ⎦ ⎥

原子和y的几何意义如图1：

正交匹配追踪

计算残差

现在已经选定了第一个基b1=[−0.7070.707]，贡献值系数为λ1=−1.34. 将该部分贡献值从y中减去，残差为

r = y - λ 1 \cdot b 1 = [1.65 - 0.25] - (- 1.34) \cdot [- 0.707 0.707] = [0.7 0.7]

残差的几何意义如图2：

正交匹配追踪

重复迭代

经过第一次迭代，我们选择出了基b1. 将选择出的基加入新的矩阵Anew，即：

A n e w = [b 1] = [- 0.707 0.707]

贡献系数矩阵可以写成xrec=⎡⎣⎢−1.3400⎤⎦⎥

第一个元素是第一个贡献值-1.34，对应着第一个第一个基b1.

残差值为r=[0.70.7]

现在需要从剩下的基b2和b3中选择对残差r贡献最大的.

w = [b 2 b 3] T \cdot r = [0.98 - 0.7]

所以，b2的贡献值更大一些是0.98.

我们将选择的基b2添加到Anew得到

A n e w = [b 1 b 2] = [- 0.707 0.707 0.8 0.6]

接下来的步骤和MP算法稍微有些不同. 这里，我们需要计算Anew对y的贡献（而不是b2对r的贡献）

我们使用最小二乘法（Least Square Problem）解决该该贡献值问题：

如何求得λ1和λ2的值使得λ1⋅[−0.7070.707]+λ2⋅[0.80.6]最接近y=[1.65−0.25]？

表示成数学语言为：

min ∥ A n e w \cdot λ - y ∥ 2

在这个例子中，即

min ∥ [b 1 b 2] = [- 0.707 0.707 0.8 0.6] \cdot [λ 1 λ 2] - [1.65 - 0.25] ∥ 2

对于min∥Anew⋅λ−y∥2的解可以直接使用公式：

λ = A + n e w \cdot y

其中，A+new是Anew的谓逆，A+new=(ATnew⋅Anew)−1⋅ATnew

（注：可参看我的博文《Moore-Penrose广义逆矩阵》和《矛盾方程的最小二乘解法》）

在我们的例子中A+new=[−0.7070.7070.80.6]+=[−0.60620.71430.80820.7143]

可以使用MATLAB中的内置函数pinv()进行伪逆的计算.

得到A+new后，可以计算λ

λ=A+new⋅y=[−0.60620.71430.80820.7143]⋅[1.65−0.25]=[−1.21]

得到λ以后，更新系数矩阵xrec为

xrec=⎡⎣⎢−1.210⎤⎦⎥

同样的，λ在xrec第一个元素和第二个元素对应着选择出来的第一个和第二个基.

得到Anew和λ以后，我们可以计算残差

r=y−Anew⋅λ=[1.65−0.25]−[−0.7070.7070.80.6]⋅[—1.21]=[00]

此时，残差值为[00]，我们可以停止计算或者继续下一步计算（这里停止，可以节省一些计算工作）.

最后一次迭代

这一步不是必须的，因为残差已经完全消除了（很多实现OMP的软件都需要输入稀疏度K参数，这样经过K次迭代以后，无论残差大小都会停止迭代）.

重新组织信号系数xrec=⎡⎣⎢−1.210⎤⎦⎥，这刚好是原始的系数.

其他例子

给定x=⎡⎣⎢⎢⎢0312⎤⎦⎥⎥⎥和A=⎡⎣⎢−0.8−0.20.20.30.411−0.3−0.10.4−0.40.8⎤⎦⎥，

求得y=A⋅x=⎡⎣⎢2.70.14.5⎤⎦⎥

现在给定A和y，使用OMP算法求解x.

我们有4个基（原子）：b1=⎡⎣⎢−0.8−0.20.2⎤⎦⎥b2=⎡⎣⎢0.30.41⎤⎦⎥b3=⎡⎣⎢1−0.3−0.1⎤⎦⎥b4=⎡⎣⎢0.4−0.40.8⎤⎦⎥
由于基向量的长度不是1，所以我们首先进行标准化.

$b 1^= b 1 / ∥ b 1 ∥ = ⎡ ⎣ ⎢ - 0.8 - 0.2 0.2 ⎤ ⎦ ⎥ / (- 0.8) 2 + (- 0.4) 2 + (0.2) 2 - - - - - - - - - - - - - - - - - - - - - - \sqrt = ⎡ ⎣ ⎢ - 0.9428 - 0.2357 0.2357 ⎤ ⎦ ⎥$

$b 2^= b 2 / ∥ b 2 ∥ = ⎡ ⎣ ⎢ 0.2680 0.3578 0.8940 ⎤ ⎦ ⎥$

$b 3^= b 3 / ∥ b 3 ∥ = ⎡ ⎣ ⎢ 0.9535 - 0.2860 0.0953 ⎤ ⎦ ⎥$

$b 2^= b 2 / ∥ b 2 ∥ = ⎡ ⎣ ⎢ 0.4082 - 0.4082 - 0.8165 ⎤ ⎦ ⎥$
标准化的基向量对y的贡献w

$w = A^T \cdot y = ⎡ ⎣ ⎢ - 0.9428 - 0.2357 0.2357 0.2680 0.3578 0.9840 0.9535 - 0.2860 - 0.0953 0.4082 - 0.4082 - 0.8165 ⎤ ⎦ ⎥ T \cdot ⎡ ⎣ ⎢ 2.7 0.1 4.5 ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ - 1.5085 4.7852 2.1167 4.7357 ⎤ ⎦ ⎥ ⎥ ⎥$
第二个基向量b2贡献值最大，所以将b2加入到Anew

Anew=b2=⎡⎣⎢0.30.41⎤⎦⎥
利用最小二乘法计算xrec

Lp=A+new⋅y=[4.28]

因为Lp对应着第二个基向量，所以xrec=⎡⎣⎢⎢⎢04.2800⎤⎦⎥⎥⎥
接下来计算残差

$r = y - A n e w \cdot L p = ⎡ ⎣ ⎢ 2.7 0.1 4.5 ⎤ ⎦ ⎥ - ⎡ ⎣ ⎢ 0.3 0.4 1 ⎤ ⎦ ⎥ \cdot 4.28 = ⎡ ⎣ ⎢ 1.416 - 1.612 0.22 ⎤ ⎦ ⎥$
接下来重复第3步，计算b1^，b1^和b1^对r的贡献.

$w = A^T \cdot r = ⎡ ⎣ ⎢ - 0.9428 - 0.2357 0.2357 0.9535 - 0.2860 - 0.0953 0.4082 - 0.4082 - 0.8165 ⎤ ⎦ ⎥ T \cdot ⎡ ⎣ ⎢ 1.416 - 1.612 0.22 ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ - 0.9032 1.7902 1.4158 ⎤ ⎦ ⎥$

选择第二个贡献最大的基b3，其贡献值为1.7902.
将选择的b3加入到Anew中

Anew=[b1b2]=⎡⎣⎢0.30.411−0.3−0.1⎤⎦⎥
用最小二乘法计算Lp=A+new⋅y=[4.17021.7149]

这个Lp对应着b2和b3，因此xrec=⎡⎣⎢⎢⎢04.1721.71490⎤⎦⎥⎥⎥
接着计算残差
$r = y - A n e w \cdot L p = ⎡ ⎣ ⎢ 2.7 0.1 4.5 ⎤ ⎦ ⎥ - ⎡ ⎣ ⎢ 0.3 0.4 1 1 - 0.3 - 0.1 ⎤ ⎦ ⎥ \cdot [4.172 1.7149] = ⎡ ⎣ ⎢ - 0.266 - 1.0536 0.5012 ⎤ ⎦ ⎥$
重复步骤2，进行最后一次迭代
计算b1和b4的贡献值w=[b1^b2^]⋅r=⎡⎣⎢−0.9428−0.23570.23570.4082−0.4082−0.8165⎤⎦⎥T⋅⎡⎣⎢−0.266−1.05360.5012⎤⎦⎥=[0.61720.7308]

这里b4提供了最大贡献值0.7308.
将b4加入Anew=⎡⎣⎢0.30.411−0.3−0.10.4−0.40.8⎤⎦⎥
利用最小二乘计算Lp=A+new⋅y=⎡⎣⎢321⎤⎦⎥，对应的xrec=⎡⎣⎢⎢⎢0312⎤⎦⎥⎥⎥
接着计算残差
$r = y - A n e w \cdot L p = ⎡ ⎣ ⎢ 2.7 0.1 4.5 ⎤ ⎦ ⎥ - ⎡ ⎣ ⎢ 0.3 0.4 1 1 - 0.3 - 0.1 0.4 - 0.4 0.8 ⎤ ⎦ ⎥ \cdot ⎡ ⎣ ⎢ 312 ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ 000 ⎤ ⎦ ⎥$
我们的迭代到此为止，因为此时残差已经为0，重建的x为xrec=⎡⎣⎢⎢⎢0312⎤⎦⎥⎥⎥，和原来的信号相同.

需要注意的问题

通过上面的迭代计算过程，我们应该注意如下几点：

OMP中最大贡献值的计算需要对基向量进行标准化处理，不是由原始基得到的.
如果给定的基向量已经是单位向量，则不需要进行标准化.
基向量对y值的贡献的计算是通过标准化以后的基向量和残差的点积进行计算的.
在MP算法中，重建系数xrec的计算是通过基向量和残差的点积进行计算的，在OMP算法中，xrec的计算是通过最小二乘法得到Anew相对于y的系数得到的，即Lp的计算. Lp向量中的值直接对应于xrec的相应位置. Anew通过每次对基向量的选择得到. 这个过程是需要时间的，因此OMP比MP慢. （不是应该OMP收敛的快吗？）
残差r的计算通过原始的y值和Anew，Lp得到.
迭代的次数最多等于A矩阵的行数M，或者如果给定了稀疏度K，则迭代K次. 如果K<M，则已知的K可以加快计算结束，如果K未知，则迭代M次.

说明总结

在正交匹配追踪OMP中，残差总是与已经选择过的原子正交的。这意味着一个原子不会被选择两次，结果会在有限的几步收敛。

OMP算法步骤描述：

输入：字典矩阵A，采样向量y，稀疏度k.

输出：x的K稀疏逼近x^.

初始化：残差r0=y，索引集Λ0=∅，t=1.

循环执行步骤1-5：

找出残差r和字典矩阵的列Ai积中最大值所对应的脚标λ，即λt=argmaxi=1,⋯,N|<rt−1,Ai>|.
更新索引集Λt=Λt−1∪{λt}，记录找到的字典矩阵中的重建原子集合At=[At−1,Aλt].
由最小二乘得到x^t=argmin∣∣|y−Atx^|∣∣.
更新残差rt=y−Atx^t，t=t+1.
判断是否满足t>K，若满足，则迭代停止；若不满足，则继续执行步骤1.