压缩感知重构算法之基追踪(Basis Pursuit, BP)_基追踪bp_Aruen24的博客-程序员秘密

题目:压缩感知重构算法之基追踪(Basis Pursuit, BP)

        除匹配追踪类贪婪迭代算法之外,压缩感知重构算法另一大类就是凸优化算法或最优化逼近方法,这类方法通过将非凸问题转化为凸问题求解找到信号的逼近,其中最常用的方法就是基追踪(Basis Pursuit, BP),该方法提出使用l1范数替代l0范数来解决最优化问题,以便使用线性规划方法来求解[1]。本篇我们就来讲解基追踪方法。

        理解基追踪方法需要一定的最优化知识基础,可参见最优化方法分类中的内容。

1、l1范数和l0范数最小化的等价问题

        在文献【2】的第4部分,较为详细的证明了l1范数与l0范数最小化在某条件下等价。证明过程是一个比较复杂的数学推导,这里尽量引用文献中的原文来说明。

        首先,在文献【2】的4.1节,给出了(P1)问题,并给出了(P1)的线性规划等价形式(LP),这个等价关系后面再详叙。

        然后在文献【2】的4.2节直接谈到l1l0最小化的关系,先是定义了压缩感知要解决的(P0)问题,然后指出“当(P0)有一个稀疏解,(P1)会找到这个解”,若并在Theorem 8中以定理形式指出“(P0)和(P1)都有相同的惟一解”。

        接下来是一段为了证明Theorem 8过渡性的描述,里面提到l1l0最小化的等价问题已经有很多文献了。


        为了证明Theorem 8,引入了一个引理Lemma4.1 :

        证明完Lemma 4.1后,开始证明Theorem 8 :


        证明过程还是比较复杂的,有兴趣的好好学习研究一下吧。

2、基追踪实现工具箱l1-MAGIC

        若要谈基追踪方法的实现,就必须提到l1-MAGIC工具箱(工具箱主页:http://users.ece.gatech.edu/~justin/l1magic/),在工具箱主页有介绍:L1-MAGIC is a collection of MATLAB routines for solving the convexoptimization programs central to compressive sampling. The algorithms are basedon standard interior-point methods, and are suitable for large-scale problems.

        另外,该工具箱专门有一个说明文档《l1-magic: Recovery of Sparse Signals via Convex Programming》,可以在工具箱主页下载。

        该工具箱一共解决了七个问题,其中第一个问题即是Basis Pursuit :

        工具箱中给出了专门针对(P1)的代码l1eq_pd.m,使用方法可以参见l1eq_example.m,说明文档的3.1节也进行了介绍。

        在附录A中,给出了将(P1)问题转化为线性规划问题的过程,但这个似乎并不怎么容易看明白:


3、如何将(P1)转化为线性规划问题?

        尽管在l1-MAGIC给出了一种基追踪的实现,但需要基于它的l1eq_pd.m文件,既然基追踪是用线性规划求解,那么就应该可以用MATLAB自带的linprog函数求解,究竟该如何将(P1)转化为标准的线性规划问题呢?我们来看文献【3】的介绍:


        这里,文献【3】的转化说明跟文献【2】中4.1节的说明差不多,但对初学者来说仍然会有一定的困难,下面我们就以文献【3】中的符号为准来解读一下。

        首先,式(3.1)中的变量a没有非负约束,所以要将a变为两个非负变量u和v的差a=u-v,由于u可以大于也可以小于v,所以a可以是正的也可以是负的[4]。也就是说,约束条件Φa=s要变为Φ(u-v)=s,而这个还可以写为[Φ,-Φ][u;v]=s,更清晰的写法如下:

        然后,根据范数的定义,目标函数可进一点写为:

        目标函数中有绝对值,怎么去掉呢?这里得看一下文献【5】:

        到现在一切应该都清晰明白了,总结如下:

        问题可以转化为线性规划问题,其中:

求得最优化解x0后可得变量a的最优化解a0=x0(1:p)-x0(p+1:2p) 。

4、基于linprog的基追踪MATLAB代码(BP_linprog.m)

[plain]  view plain   copy
  在CODE上查看代码片 派生到我的代码片
  1. function [ alpha ] = BP_linprog( s,Phi )  
  2. %BP_linprog(Basis Pursuit with linprog) Summary of this function goes here  
  3. %Version: 1.0 written by jbb0523 @2016-07-21   
  4. %Reference:Chen S S, Donoho D L, Saunders M A. Atomic decomposition by  
  5. %basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at:   
  6. %http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)  
  7. %   Detailed explanation goes here  
  8. %   s = Phi * alpha (alpha is a sparse vector)    
  9. %   Given s & Phi, try to derive alpha  
  10.     [s_rows,s_columns] = size(s);    
  11.     if s_rows<s_columns    
  12.         s = s';%s should be a column vector    
  13.     end   
  14.     p = size(Phi,2);  
  15.     %according to section 3.1 of the reference  
  16.     c = ones(2*p,1);  
  17.     A = [Phi,-Phi];  
  18.     b = s;  
  19.     lb = zeros(2*p,1);  
  20.     x0 = linprog(c,[],[],A,b,lb);  
  21.     alpha = x0(1:p) - x0(p+1:2*p);  
  22. end  

5、基追踪单次重构测试代码(CS_Reconstuction_Test.m)

        测试代码与OMP测试单码相同,仅仅是修改了重构函数。

[plain]  view plain   copy
  在CODE上查看代码片 派生到我的代码片
  1. %压缩感知重构算法测试    
  2. clear all;close all;clc;    
  3. M = 64;%观测值个数    
  4. N = 256;%信号x的长度    
  5. K = 10;%信号x的稀疏度    
  6. Index_K = randperm(N);    
  7. x = zeros(N,1);    
  8. x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的    
  9. Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta    
  10. Phi = randn(M,N);%测量矩阵为高斯矩阵    
  11. A = Phi * Psi;%传感矩阵    
  12. y = Phi * x;%得到观测向量y    
  13. %% 恢复重构信号x    
  14. tic    
  15. theta = BP_linprog(y,A);    
  16. x_r = Psi * theta;% x=Psi * theta    
  17. toc    
  18. %% 绘图    
  19. figure;    
  20. plot(x_r,'k.-');%绘出x的恢复信号    
  21. hold on;    
  22. plot(x,'r');%绘出原信号x    
  23. hold off;    
  24. legend('Recovery','Original')    
  25. fprintf('\n恢复残差:');    
  26. norm(x_r-x)%恢复残差   

        运行结果如下:(信号为随机生成,所以每次结果均不一样)

         1)图:

2)Command Windows

Optimization terminated.

Elapsed time is 0.304111 seconds.

恢复残差:

ans =

  6.5455e-010

6、结束语

        值得一提的是,基追踪并不能称为一个具体的算法,而是一种最优化准则,文献【3】对此进行了明确的说明,基追踪实现方法可以使用单纯形法(simplex algorithm),也可以使用内点法(interior-pointmethods), 因此,有些文献里说凸松弛算法包括基追踪、内点法等,个人感觉这是不恰当的,因为内点法只是基追踪的一种实现形式而己,再说了,内点法也有很多种实现方法……

        本文实现方法基于MATLAB自带的线性规划函数linprog,当然也可以采用l1-magic中的l1eq_pd.m,有兴趣的可以做一下对比。

7、参考文献

【1】李珅, 马彩文, 李艳, 等. 压缩感知重构算法综述[J]. 红外与激光工程, 2013, 42(S01): 225-232.

【2】Donoho D L. Compressedsensing[J]. IEEE Transactions on information theory, 2006, 52(4):1289-1306. (Available at: http://www.signallake.com/innovation/CompressedSensing091604.pdf)

【3】Chen S S, Donoho D L,Saunders M A.Atomicdecomposition by basis pursuit[J]. SIAM review, 2001, 43(1): 129-159. (Availableat:http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)

【4】孙文瑜, 徐成贤, 朱德通.最优化方法(第二版)[M]. 北京:高等教育出版社, 2010:49-51.

【5】L1范数优化的线性化方法如何证明? 链接:http://www.zhihu.com/question/21427075

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/wyw921027/article/details/52101780

智能推荐

hibernate在控制台打印 SQL 语句_hibernate 打印sql_春风化作秋雨的博客-程序员秘密

在 Config 里面把 hibernate. show_SQL 设置为 true 即可。不建议开启,开启后会降低程序的运行效率。1、spring boot 之ymlspring: jpa: properties: hibernate: format_sql: true //格式化sql语句 show_sql: true //控制台是否打印 use_sql_comments: true /

Model Degradation Hinders Deep Graph Neural Networks(AIR)论文阅读笔记_刘大彪的博客-程序员秘密

Model Degradation Hinders Deep Graph Neural Networks(AIR)论文阅读笔记

影刀RPA宣布完成Coatue领投5000万美金B轮融资,高瓴、金沙江等老股东跟投投_数据猿的博客-程序员秘密

数据智能产业创新服务媒体——聚焦数智· 改变商业“影刀RPA在等待IT精英们,如同蒲公英一般的伞兵,在黑夜里从天而降,长驱直入,用最智慧的产品、最优质的服务一起干一家伟大的科技公司,让生...

隐藏的模块中编译错误:ThisWorkbook._程序员正茂的博客-程序员秘密

打开Excel时总提示隐藏的模块中编译错误:ThisWorkbook.解决办法:关掉所有可用的加载宏。

Leetcode算法——54、螺旋矩阵(spiral matrix)_从左上角开始,以顺时针螺旋的形式依次排列的结果就是_HappyRocking的博客-程序员秘密

给定一个矩阵 m*n,返回所有元素的螺旋排列顺序(从左上角开始,顺时针旋转,由外向内)。示例1:Input:[ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ]]Output: [1,2,3,6,9,8,7,4,5]示例2:Input:[ [1, 2, 3, 4], [5, 6, 7, 8], [9,10,11,12]]Output...

FFT快速傅里叶算法——基2FFT时域抽取基本原理解析_基2频分傅里叶算法与基2时分傅里叶算法_ysc_11的博客-程序员秘密

学习数字信号处理时学到如果直接对有限长的序列进行N点DFT,那么其需要进行N2N^2N2次的复数乘以及复数加运算,其算法复杂度十分高,将会浪费大量时间在计算过程中。为此研究了FFT(快速傅里叶算法)这里给大家介绍两种:时域抽取法(DIT-DFT)和频域抽取法(DIF-FFT)。时域抽取法:设有一序列x(n)x(n)x(n)的长度为N,且满足N=2MN=2^MN=2M,M为自然数。按nnn的奇偶性分成两个N/2N/2N/2的子序列。x1(r)=x(2r)r=0,1,2…N/2−1x_1(r)=x(2r)

随便推点

分布式架构的十四次演进之路_公众号:ITIL之家的博客-程序员秘密

公众号回复&#39;架构&#39;获取架构师电子书及视频课程导读:本文介绍从一百个到千万级并发情况下服务端的架构的演进过程,同时列举出每个演进阶段会遇到的相关技术,让大家对架构的演进有一个...

Centos ping不通百度彻底解决_我是潇洒哥的弟弟黑大帅的博客-程序员秘密

1.简介在职场过程中,有一个重要的内容就是熟练使用Linux操作系统,操作命令,而安装虚拟机,配置虚拟机环境则是开发的首要步骤。本文主要是对在VMWare中配置网络,使得通过虚拟机可以访问外部网络,ping通www.baidu.com而写。2.安装步骤2.1 安装VMware VMWare下载2.2 安装Centos2.2.1 安装Centos其他的参数不再赘述。如果在安装时提示Intel-VTx未启用,则重启,并在BIOS中开启Intel-VTx.在安装Centos,使用1708版

【原创】ubuntu 挂起后唤醒解决办法_ubuntu挂起后怎么唤醒_GENGLUT的博客-程序员秘密

待机计算机将目前的运行状态等数据存放在内存,关闭硬盘、外设等设备,进入等待状态。此时内存仍然需要电力维持其数据,但整机耗电很少。恢复时计算机从内存读 出数据,回到挂起前的状态,恢复速度较快。一般笔记本在电池无故障且充满的情况下可以支持这种挂起数小时甚至数天(依具体机型有差别)。其他名称:Suspend, STR(Suspend To RAM), 挂起, 挂起到内存休眠计算机将目前的

【C#】一次性清空textbox、combobox中所有的内容_c#窗体清空整个页面数据_光哥_帅的博客-程序员秘密

在做机房重构的时候,总是会遇到清空所有的代码,比如注册的窗体,如果你每个窗体的清空都写成,像这样:txtcard.text=”“;这样就会出现大量的冗余的代码!这时请看下面的代码,批量清除了所有的控件的内容,如果此时用到这个方法的窗体过多,就可以将它封装成一个类! //清空所有控件里边的内容 private void btnClear_Click(object s...

Intellij idea设置调用方法提示方法注释(2018版)_idea 2018显示方法的引用_GuawaUp的博客-程序员秘密

如何在Intellij idea中像Eclipse一样再调用方法时出现方法的注释?效果图如下设置步骤如下:file-settings-Editor-General-CodeCompletion-&gt;勾选上Auto-displaydocumentation选项即可此方法适用于Intellij idea2018版 其他版本可以参考https://blog.csdn.net...

iqooneo5隐藏应用方法分享(2021)_iqooneo5怎么分享已安装软件_爱学习的小兔子的博客-程序员秘密

在我们的手机中难免会有一些隐私应用不想被别人看到。这时我们可以对该软件进行应用隐藏。那具体该怎么操作呢?不要着急换换早就为大家准备好了iqooneo5应用隐藏的详细方法。快来接着往下看看吧!iqooneo5隐藏应用方法分享1、打开手机设置。点击【指纹、面部与密码】。2、点击【隐私与应用加密】。设置好隐私密码。3、点击【应用隐藏】。4、开启需要隐藏应用的开关即可。...

推荐文章

热门文章

相关标签