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]可能容易理解:在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的值是由求法很简单,通过对(7)左右两边添加作内积消减得到:


后边的第二项因为它们正交,所以为0,所以可以得出ak的第一部分。对于,在(4)左右两边中与作内积,可以得到ak的第二部分。

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

3.2 收敛性证明

通过(7),由于正交,将两个残值移到右边后求二范的平方,并将ak的值代入可以得到:


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

3.3 算法步骤

整个OMP算法的步骤如下:

由于有了上面的来龙去脉,这个算法就相当好理解了。

到这里还不算完,后来OMP的迭代运算用另外一种方法可以计算得知,有位同学的论文[2]描述就非常好,我就直接引用进来:

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

再贴另外一个洋牛paper[3]中关于OMP的描述,之所以引入,是因为它描述的非常严谨,但是也有点苦涩难懂,不过有了上面的基础,就容易多了。


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


它的算法步骤Update Provisional Solution中求很简单,就是在 b = Ax 已知 A和b求x, 在x的最小二范就是A的伪逆与b相乘,即:


看上去头疼,其实用matlab非常简单,看看上面的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


本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/452251.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【牛客刷题记录】2021-03-10

牛客代码刷题记录1 问题一 给定一个数组序列, 需要求选出一个区间, 使得该区间是所有区间中经过如下计算的值最大的一个&#xff1a; 区间中的最小数 * 区间所有数的和最后程序输出经过计算后的最大值即可&#xff0c;不需要输出具体的区间。如给定序列 [6 2 1]则根据上述公式…

SpringMVC详细示例实战教程

一、SpringMVC基础入门&#xff0c;创建一个HelloWorld程序 1.首先&#xff0c;导入SpringMVC需要的jar包。 2.添加Web.xml配置文件中关于SpringMVC的配置 123456789101112131415<!--configure the setting of springmvcDispatcherServlet and configure the mapping-->&…

学成在线--8.Freemarker入门教程

文章目录1.FreeMarker介绍1&#xff09;常用的java模板引擎还有哪些&#xff1f;2&#xff09;freemarker是一个用Java开发的模板引擎3&#xff09;模板数据模型输出2.FreeMarker快速入门1&#xff09;创建测试工程2&#xff09;配置文件3&#xff09;创建模型类4&#xff09;创…

面向对象建模

用面向对象方法开发软件&#xff0c;通常需要建立3种形式的模型&#xff0c;它们分别是 描述系统数据结构的对象模型 &#xff0c; 描述系统控制结构的动态模型 和 描述系统功能的功能模型 。一个典型的软件系统使用数据结构(对象模型)&#xff0c;执行操作(动态模型)&#xff…

if...elif...else...fi和case...esac的脚本条件判断式

注意1&#xff1a; if 表达式和case 表达式的区别及什么时候使用哪个要有明确的区分&#xff1f; 『 if .... then .... fi 』对于变量的判断是以『比对』的方式来分辨的&#xff0c; 如果符合状态就进行某些行为&#xff0c;并且透过较多层次 (就是elif ) 的方式来进行多个变量…

安卓自定义View进阶-分类与流程

自定义View绘制流程函数调用链(简化版) 一.自定义View分类 我将自定义View分为了两类(sloop个人分类法&#xff0c;非官方)&#xff1a; 1.自定义ViewGroup 自定义ViewGroup一般是利用现有的组件根据特定的布局方式来组成新的组件&#xff0c;大多继承自ViewGroup或各种Layout&…

【笔试记录】2021/3/13美团

2021/3/13美团笔试 1.旋转二维数组 m, n map(int, input().split()) a [] for _ in range(m):a.append(list(map(int, input().split()))) for j in range(n):for k in range(m):print(a[k][j], end )print(end\n)输入输出 2.找出输入字符串中的所有数字 s input() n …

Ubuntu品牌机批量涌入世界市场,中国不会例外

根据种种信息来源&#xff0c;我们可以预见&#xff1a; Ubuntu 品牌机不久即将涌入世界市场&#xff0c;中国也不会例外。对此现象&#xff0c;我们必须有一定的思想准备。 2004年4月&#xff0c;在澳大利亚一次学术会议上萌发了Ubuntu发行版的念头&#xff0c;至今已经有8个年…

学成在线--9.页面静态化

文章目录一.页面静态化流程二.数据模型1.轮播图DataUrl接口1&#xff09;需求分析2&#xff09;接口定义3&#xff09;Dao4&#xff09;Service5&#xff09;Controller6&#xff09;测试2.远程请求接口1&#xff09;添加依赖2&#xff09;配置RestTemplate3&#xff09;测试Re…

数据库的完整性和安全性

数据的完整性防止数据库中存在不符合语义的数据&#xff0c;也就是防止数据库中存在不正确的数据防范对象&#xff1a;不合语义的、不正确的数据数据的安全性保护数据库 防止恶意的破坏和非法的存取防范对象&#xff1a;非法用户和非法操作

数据库MySQL/mariadb知识点——日志记录(2)二进制日志

二进制日志 记录已提交事务导致数据改变或潜在导致数据改变的SQL语句&#xff0c;通过“重放”日志文件中的事件来生成数据副本&#xff0c;不依赖存储引擎类型。 开启二进制日志&#xff0c;默认是关闭的&#xff0c;二进制日志和数据分开存放 开启记录二进制文件的功能&#…

【面试记录】Python常见面试200题及答案总结

Python常见面试200题及答案总结 /待完善/ 1. 列出5个常用python标准库&#xff1f; os&#xff1a;提供了不少与操作系统相关联的函数&#xff0c;提供了一种可移植的使用操作系统功能的方法。使用os模块中提供的接口&#xff0c;可实现跨平台访问。但是&#xff0c;并不是所…

Linux负载均衡软件LVS之一(概念篇)

2019独角兽企业重金招聘Python工程师标准>>> 一、 LVS简介 LVS是Linux Virtual Server的简称&#xff0c;也就是Linux虚拟服务器, 是一个由章文嵩博士发起的自由软件项目&#xff0c;它的官方站点是www.linuxvirtualserver.org。现在LVS已经是 Linux标准内核的一部分…

C语言目录操作 (Linux/Unix)

获得工作目录&#xff1a;#include <unistd.h> char *getcwd(char *buf,size_t size);char *getwd(char *buf);/*this is for FreeBSD*/ 改变当前目录&#xff1a;#include <unistd.h> int chdir(const char *path); 保存当前目录&#xff1a;#include <unist…

学成在线--10.页面预览

文章目录一.需求分析二. 搭建环境1.在cms服务中集成freemarker&#xff1a;2.在application.yml配置freemarker三.Service四.Controller五.页面预览测试-- 配置Nginx代理六.页面预览测试-- 添加“页面预览”链接一.需求分析 页面在发布前增加页面预览的步骤&#xff0c;方便用…

参照完整性违约处理

参照完整性违约处理 &#xff08;1&#xff09; 拒绝&#xff08;NO ACTION&#xff09;执行 不允许该操作执行。该策略一般设置为默认策略 &#xff08;2&#xff09; 级联&#xff08;CASCADE&#xff09;操作 当删除或修改被参照表&#xff08;Student&#xff09;的一个元…

BSON和JSON的区别

BSON目前主要用于MongoDB中&#xff0c;是mongoDB的数据存储格式&#xff0c;BSON基于JSON格式&#xff0c;选择JSON进行改造的原因主要是JSON的通用性及JSON的schemaless的特性。 BSON主要实现下面三点&#xff1a; 1、更快的遍历速度 对json格式来说&#xff0c;太大的json结…

【代码刷题】排序算法总结(python实现)

排序算法总结&#xff08;Python实现&#xff09;算法介绍算法分类相关概念1. 冒泡排序&#xff08;Bubble Sort&#xff09;1.1 思想1.2 python实现1.3 复杂度1.4 稳定性2. 快速排序&#xff08;Quick Sort&#xff09;2.1 思想&#xff08;伪代码&#xff09;2.2 python实现2…

C语言遍历目录

C语言遍历目录&#xff0c;可以循环的遍历子目录#include <stdio.h>#include <string.h>#include <stdlib.h>#include <dirent.h>#include <sys/stat.h>#include <unistd.h>#include <sys/types.h>void listDir(char *path){struct …

学成在线--11.RabbitMQ快速入门

文章目录一.RabbitMQ简介二.相关知识1.AMQP2.JMS是什么 &#xff1f;三.RabbitMQ的工作原理四.Hello World1.创建Maven工程2.生产者3.消费者五.总结一.RabbitMQ简介 MQ全称为Message Queue&#xff0c;即消息队列&#xff0c; RabbitMQ是由erlang语言开发&#xff0c;基于AMQP…