算法代码中的循环矩阵在哪体现_机器学习基础（二）——匹配追踪算法（Matching Pursuit）...

最新推荐文章于 2021-12-10 15:58:55 发布

股票期货软件开发

最新推荐文章于 2021-12-10 15:58:55 发布

阅读量281

点赞数

文章标签：算法代码中的循环矩阵在哪体现

本文链接：https://blog.csdn.net/weixin_35680189/article/details/112730057

版权

之前我们介绍过了什么是稀疏表示，用一种浅显的方式介绍了稀疏表示和它的应用。得出了一个简单的稀疏表示公式：

(1)

那么自然我们就会得出一个问题，现在我们已知信号y，以及一个过完备字典D，如何才能求解出这个信号的稀疏表示α呢？本文要介绍的匹配追踪算法（Matching Pursuit， MP）就是最经典的稀疏表示求解算法之一。网络上有很多关于MP算法的资料，笔者开始学的时候也一脸懵逼，今天就介绍我对于MP算法的理解，希望对刚入门的同学有所帮助。

（〇）前言

首先我们要知道，对于（1）式存在的不等式约束的优化问题，我们可以将其优化式写作其增广拉格朗日（Augmented Lagrangian method）形式：

（2）

式（2）为我们具体要求解的优化式，但是式中存在非凸项

让整体优化式变为了

NP-hard问题，没有解析解，只能用穷举的办法求解。简单来说就是不能直接利用公式推导求梯度的方式求解，因此只能用一些歪门邪道求解式（2）的近似解。 MP算法就是这样一种求解近似解的方法。

（一）内积

对于一段信号

，字典

，我们首先要解决的问题是如何才能找到字典

D中与 y最为匹配的原子。为此我们有多种解决方案，如相关性分析（ correlation analysis），残差（ residual）等，MP算法选择用内积（ inner product）来计算两个原子的相关性，相较于残差内积具有更好的数学性质，相较于相关函数内积的计算结果更加简单稳定。假设存在四个长度为5的信号， y1=[1,2,3,4,5]， y2=[1,2,3,4,5]， y3=[5,3,4,1,2]， y4=[1,2,3,5,4]，下表可以看出以上三个指标作为评判标准下 y1与其他三个信号的匹配程度。

	[y1,y2]	[y1,y3]	[y1,y4]
相关性	1	0.8	0.9
残差	0	0	0
内积	55	37	54

可以看出内积是更加灵敏的相关性分析工具。

实际上在MP算法内部，就是利用信号y与字典D的内积结果最大位置处的原子作为当前循环寻找的最优解，不难看出MP算法每次只寻找当前最优解，是一个典型的贪婪算法。具体来说：

（3）

其中

为内积计算结果，取内积向量

h中最大值处索引 i处所对应字典 D中的原子即为当前循环中与信号 y最为匹配的原子，记该原子为

。

（二）支撑集

在第k次循环中，理论上我们已经得到了k个原子

，

i=1.2…k。由于字典是固定不变的，因此将 k个原子的在字典中的索引放在矩阵

中，矩阵t就被称为支撑集。同时每次循环都能得出一个最大内积

h*，同样存入一个矩阵

中。因此在第

k次循环中利用支撑集 t找到字典 D中对应的原子，并将其构成一个矩阵

，并计算信号

y与

和矩阵

s内积的残差得到第 k次循环的残差

(4)

将作为第

k+1次循环的输入信号，寻找字典 D中与

最匹配的原子。

（三）迭代终止条件

现在是时候考虑如何退出循环了。我们的目标是找到字典D中最匹配信号y的k个原子，因此最朴素的思想就是当支撑集中原子的数目大于常数c时退出循环，即固定稀疏度的MP算法，这样的迭代终止条件很常用。但实际上一段信号中信息量相差很大，这样一刀切的方式很难实现所有信号信息重构，因此我们通常用残差的二范数大小来充当迭代终止条件，即

。实际上从利用残差二范数的方式的迭代结果来看一段信号的最终支撑集长度通常也仅为

2-3。

（四）正交匹配追踪（Orthogonal Matching Pursuit，OMP）算法

MP算法有一个很致命的问题，就是支撑集的重复选择。即在第k次循环中找到的最佳原子

与第

L次循环找到的最佳原子

可能是同一个原子。这使得

MP算法的收敛速度大幅降低，从而学者提出了正交匹配追踪算法（ OMP），他能保证在每次计算得出的残差与当前支撑集中原子为正交关系，从而在下次循环时支撑集中的原子与残差的内积为 0，因此不会出现原子复选的情况。在实际的运用中 OMP算法由于其更快的收敛效率受到更多学者的青睐，但事实上 OMP算法与 MP算法理论上会收敛于同一解。

OMP算法的实现形式为在每次循环的残差计算时采用最小二乘法，即舍弃式（4），首先利用最小二乘公式求解在当前支撑集下与原信号y的最小二乘解。

(5)

再利用最小二乘解β计算残差r。

(6)

除此之外与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