MP算法和OMP算法及其思想

主要介绍MP(Matching Pursuits)算法和OMP(Orthogonal Matching Pursuit)算法[1]，这两个算法尽管在90年代初就提出来了，但作为经典的算法，国内文献(可能有我没有搜索到)都仅描写叙述了算法步骤和简单的应用，并未对其进行详尽的分析，国外的文献还是分析的非常透彻，所以我结合自己的理解，来分析一下写到博客里，算作笔记。

1. 信号的稀疏表示(sparse representation of signals)

给定一个过完备字典矩阵，当中它的每列表示一种原型信号的原子。给定一个信号y，它能够被表示成这些原子的稀疏线性组合。信号 y 能够被表达为 y = Dx ，或者。 字典矩阵中所谓过完备性，指的是原子的个数远远大于信号y的长度(其长度非常显然是n)，即n<<k。

2.MP算法(匹配追踪算法)

2.1 算法描写叙述

作为对信号进行稀疏分解的方法之中的一个，将信号在完备字典库上进行分解。

假定被表示的信号为y，其长度为n。假定H表示Hilbert空间，在这个空间H里，由一组向量构成字典矩阵D，当中每一个向量能够称为原子(atom)，其长度与被表示信号 y 的长度n同样，并且这些向量已作为归一化处理，即|，也就是单位向量长度为1。MP算法的基本思想：从字典矩阵D（也称为过完备原子库中），选择一个与信号 y 最匹配的原子(也就是某列)，构建一个稀疏逼近，并求出信号残差，然后继续选择与信号残差最匹配的原子，重复迭代，信号y能够由这些原子来线性和，再加上最后的残差值来表示。非常显然，假设残差值在能够忽略的范围内，则信号y就是这些原子的线性组合。假设选择与信号y最匹配的原子？怎样构建稀疏逼近并求残差？怎样进行迭代？我们来具体介绍使用MP进行信号分解的步骤：[1] 计算信号 y 与字典矩阵中每列(原子)的内积，选择绝对值最大的一个原子，它就是与信号 y 在本次迭代运算中最匹配的。用专业术语来描写叙述：令信号，从字典矩阵中选择一个最为匹配的原子，满足，r0 表示一个字典矩阵的列索引。这样，信号 y 就被分解为在最匹配原子的垂直投影分量和残值两部分，即：。[2]对残值R1f进行步骤[1]相同的分解，那么第K步能够得到：

， 当中 满足。可见，经过K步分解后，信号 y 被分解为：，当中。

2.2 继续讨论

(1)为什么要假定在Hilbert空间中？Hilbert空间就是定义了完备的内积空。非常显然，MP中的计算使用向量的内积运算，所以在在Hilbert空间中进行信号分解理所当然了。什么是完备的内积空间？篇幅有限就请自己搜索一下吧。

(2)为什么原子要事先被归一化处理了，即上面的描写叙述。内积经常使用于计算一个矢量在一个方向上的投影长度，这时方向的矢量必须是单位矢量。MP中选择最匹配的原子是，是选择内积最大的一个，也就是信号(或是残值)在原子(单位的)垂直投影长度最长的一个，比方第一次分解过程中，投影长度就是。，三个向量，构成一个三角形，且和正交（不能说垂直，可是能够想象二维空间这两个矢量是垂直的）。

(3)MP算法是收敛的，由于，和正交，由这两个能够得出，得出每个残值比上一次的小，故而收敛。

2.3 MP算法的缺点

如上所述，假设信号(残值)在已选择的原子进行垂直投影是非正交性的，这会使得每次迭代的结果并不少最优的而是次最优的，收敛须要非常多次迭代。举个样例说明一下：在二维空间上，有一个信号 y 被 D=[x1, x2]来表达，MP算法迭代会发现总是在x1和x2上重复迭代，即，这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。再用严谨的方式描写叙述[1]可能easy理解:在Hilbert空间H中，，，定义，就是它是这些向量的张成中的一个，MP构造一种表达形式：;这里的Pvf表示 f在V上的一个正交投影操作，那么MP算法的第 k 次迭代的结果能够表演示样例如以下(前面描写叙述时信号为y，这里变成f了，请注意)：

假设  是最优的k项近似值，当且仅当。因为MP仅能保证，所以普通情况下是次优的。这是什么意思呢？是k个项的线性表示，这个组合的值作为近似值，仅仅有在第k个残差和正交，才是最优的。假设第k个残值与正交，意味这个残值与fk的随意一项都线性无关，那么第k个残值在后面的分解过程中，不可能出现fk中已经出现的项，这才是最优的。而普通情况下，不能满足这个条件，MP一般仅仅能满足第k个残差和xk正交，这也就是前面为什么提到“信号(残值)在已选择的原子进行垂直投影是非正交性的”的原因。假设第k个残差和fk不正交，那么后面的迭代还会出现fk中已经出现的项，非常显然fk就不是最优的，这也就是为什么说MP收敛就须要很多其它次迭代的原因。不是说MP一定得到不到最优解，并且其前面描写叙述的特性导致一般得到不到最优解而是次优解。那么，有没有办法让第k个残差与正交，方法是有的，这就是以下要谈到的OMP算法。

3.OMP算法

3.1 算法描写叙述

OMP算法的改进之处在于：在分解的每一步对所选择的所有原子进行正交化处理，这使得在精度要求同样的情况下，OMP算法的收敛速度更快。

那么在每一步中怎样对所选择的所有原子进行正交化处理呢？在正式描写叙述OMP算法前，先看一点基础思想。

先看一个 k  阶模型，表示信号 f 经过 k 步分解后的情况，似乎非常眼熟，但要注意它与MP算法不同之处，它的残值与前面每一个分量正交，这就是为什么这个算法多了一个正交的原因，MP中仅与近期选出的的那一项正交。

(1)

 k + 1 阶模型例如以下：

(2)

应用 k + 1阶模型减去k 阶模型，得到例如以下：

(3)

我们知道，字典矩阵D的原子是非正交的，引入一个辅助模型，它是表示对前k个项的依赖，描写叙述例如以下：

(4)

和前面描写叙述类似，在span(x1, …xk)之中的一个上的正交投影操作，后面的项是残值。这个关系用数学符号描写叙述：

请注意，这里的 a 和 b 的上标表示第 k 步时的取值。

将(4)带入(3)中，有：

（5）

假设一下两个式子成立，(5)必定成立。

(6)

(7)

令，有

当中。

ak的值是由求法非常easy，通过对(7)左右两边加入�作内积消减得到：

后边的第二项由于它们正交，所以为0，所以能够得出ak的第一部分。对于，在（4）左右两边中与作内积，能够得到ak的第二部分。

对于(4)，能够求出，求的步骤请參见參考文件的计算细节部分。为什么这里不提，由于后面会介绍更简单的方法来计算。

3.2 收敛性证明

通过(7)，因为与正交，将两个残值移到右边后求二范的平方，并将ak的值代入能够得到：

可见每一次残差比上一次残差小，可见是收敛的。

3.3 算法步骤

整个OMP算法的过程例如以下：

因为有了上面的来龙去脉，这个算法就相当好理解了。

到这里还不算完，后来OMP的迭代运算用第二种方法能够计算得知，有位同学的论文[2]描写叙述就很好，我就直接引用进来：

对照中英文描写叙述，本质都是一样，仅仅是有细微的区别。这里顺便贴出网一哥们写的OMP算法的代码，源出处不得而知，共享给大家。

再贴另外一个洋牛paper[3]中关于OMP的描写叙述，之所以引入，是由于它描写叙述的很严谨，可是也有点苦涩难懂，只是有了上面的基础，就easy多了。

它的描写叙述中的Sweep步骤就是寻找与当前残差最大的内积时列在字典矩阵D中的索引，它的这个步骤描写叙述说明为什么要选择内积最大的以及怎样选择。见下图，说的很清晰。

它的算法步骤Update Provisional Solution中求非常easy，就是在 b = Ax 已知 A和b求x， 在x的最小二范就是A的伪逆与b相乘，即：

看上去头疼，其有用matlab很easy，看看上面的matlab的代码就明确了。

我们能够看得出来，算法流程清晰明了，还是非常好理解的。这正是OMP算法的魅力，作为工具使用简单，背后却隐藏着非常有趣的思想。

写这篇博客的目的，是由于搜索了一下，MP和OMP没有人非常具体的介绍。文献[1]讲的非常清楚的，大家有兴趣能够找来看看。不要被老板发现我竟然在搜中文文献还写中文博客。

參考文献：

[1] Orthogonal Matching Pursuit:Recursive Function Approximat ion with Applications to Wavelet Decomposition
[2]http://wenku.baidu.com/view/22f3171614791711cc7917e4.html

[3] From Sparse Solutions of Systems of Equations to Sparse Modeling of Signals and Images