之前我们介绍过了什么是稀疏表示,用一种浅显的方式介绍了稀疏表示和它的应用。得出了一个简单的稀疏表示公式:
那么自然我们就会得出一个问题,现在我们已知信号y,以及一个过完备字典D,如何才能求解出这个信号的稀疏表示α呢?本文要介绍的匹配追踪算法(Matching Pursuit, MP)就是最经典的稀疏表示求解算法之一。网络上有很多关于MP算法的资料,笔者开始学的时候也一脸懵逼,今天就介绍我对于MP算法的理解,希望对刚入门的同学有所帮助。
(〇)前言
首先我们要知道,对于(1)式存在的不等式约束的优化问题,我们可以将其优化式写作其增广拉格朗日(Augmented Lagrangian method)形式:
式(2)为我们具体要求解的优化式,但是式中存在非凸项
(一)内积
对于一段信号
[y1,y2] | [y1,y3] | [y1,y4] | |
---|---|---|---|
相关性 | 1 | 0.8 | 0.9 |
残差 | 0 | 0 | 0 |
内积 | 55 | 37 | 54 |
可以看出内积是更加灵敏的相关性分析工具。
实际上在MP算法内部,就是利用信号y与字典D的内积结果最大位置处的原子作为当前循环寻找的最优解,不难看出MP算法每次只寻找当前最优解,是一个典型的贪婪算法。具体来说:
其中
(二)支撑集
在第k次循环中,理论上我们已经得到了k个原子
(三)迭代终止条件
现在是时候考虑如何退出循环了。我们的目标是找到字典D中最匹配信号y的k个原子,因此最朴素的思想就是当支撑集中原子的数目大于常数c时退出循环,即固定稀疏度的MP算法,这样的迭代终止条件很常用。但实际上一段信号中信息量相差很大,这样一刀切的方式很难实现所有信号信息重构,因此我们通常用残差的二范数大小来充当迭代终止条件,即
(四)正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法
MP算法有一个很致命的问题,就是支撑集的重复选择。即在第k次循环中找到的最佳原子
OMP算法的实现形式为在每次循环的残差计算时采用最小二乘法,即舍弃式(4),首先利用最小二乘公式求解在当前支撑集下与原信号y的最小二乘解。
再利用最小二乘解β计算残差r。
除此之外与OMP算法与MP算法的流程完全一致。此外随着学者的不断研究近年来也提出了多种OMP算法的改进,包括StOMP,ROMP等。
(五)补充
实际上针对式(2),已经有多种近似求解算法,除了本文介绍的MP算法,还包括将式(2)松弛为1范数,整体转化为LASSO问题求解,可运用包括基于梯度的ADMM算法,以及阈值算法ISTA等。可以证明,松弛为1范数的式(2)同样能保证稀疏系数的稀疏度。
(六)Python实现。
在这里贴上本人作为初学者期间实现的OMP算法代码,值得一提的是在sklearn中以已经集成了OMP算法,但手撕代码的快乐只有我们能懂。
import numpy as np
def omp(y, D, S=20):
M, N=D.shape[0],D.shape[1]
result = np.zeros((N, 1))
indices = []
R = y
rr=0
for i in range(S):
projection = np.dot(D.T, R)
pos = projection.index(max(projection))
indices.append(pos)
At = D[:,indices[:]]
Ls = np.dot(np.linalg.pinv(At),y)
R = y - np.dot(At, Ls)
rr=R
if (R**2).sum()<1e-6:
break
for t, s in zip(indices, Ls):
result[t] = s
return result, rr