精华内容
下载资源
问答
  • 基追踪算法

    2017-12-14 11:23:07
    基追踪算法可用于信号的去噪,以及其他处理。最后的结果可用于故障分析等。
  • 基追踪算法,含去噪matlab,可以直接作为库使用
  • 该方案利用小波变换实现图像稀疏化,利用标准伪随机数均匀分布和二维中心傅里叶变换生成随机测量矩阵,并对小波变换后的高频子带进行加权采样,用改进的基追踪算法实现二维图像压缩感知重建。仿真实验结果表明,该...
  • 压缩感知-CS-BP算法(基追踪算法

    热门讨论 2012-12-31 13:50:36
    压缩传感 压缩感知 CS BP算法 基追踪算法 经测试可以实现数据压缩感知重构。
  • 这是一个关于利用基追踪RP算法解稀疏信号重构的代码包,在无线传感器网络中起到良好的定位作用,非常详细
  • 适合初学者学习
  • 基于MATLAB算法,基追踪匹配追踪算法,稀疏分解或者压缩感知算法,应用广泛,优化求解算法等,可利于初学者有效学习
  • 该算法首先求出并存储正交向量之间的内积,然后根据向量正交展开系数为其与正交向量内积的性质将内积运算转化为代数运算,得到一种快速匹配追踪算法。实验结果表明,基于Dirac和DCT构成的完备字典对信号...
  • 基于基追踪的压缩感知算法

    热门讨论 2013-04-04 20:34:32
    基于基追踪的压缩感知算法,非常实用,研究热点,可直接使用
  • 该代码是压缩感知重构算法基追踪(BP),注释很详细,可以直接运行
  • 实验结果表明,在一定条件下,该算法重构精度高于正交匹配追踪算法(orthogonal matching pursuit,OMP)、子空间匹配算法(subspace pursuit,SP)、基追踪算法(basis pursuit,BP)和前向后向追踪算法(forward-backward ...
  • 基于基追踪—Moore-Penrose 逆矩阵算法的稀疏信号重构
  • 而现有的基于正交匹配追踪的改进LLE算法(LLE-OMP),由于0范数的限制,实质上近邻点的选取还是不会大于高维流形维数。 本文提出了一种新的LLE改进算法,允许扩大近邻点集合,从而更好地提取高维流形的全局性几何...
  • 压缩感知重构算法基追踪(Basis Pursuit, BP)

    万次阅读 多人点赞 2016-07-21 21:03:19
    题目:压缩感知重构算法基追踪(Basis Pursuit, BP)  除匹配追踪类贪婪迭代算法之外,压缩感知重构算法另一大类就是凸优化算法或最优化逼近方法,这类方法通过将非凸问题转化为凸问题求解找到信号的逼近,其中最...

    题目:压缩感知重构算法之基追踪(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)

    function [ alpha ] = BP_linprog( s,Phi )
    %BP_linprog(Basis Pursuit with linprog) Summary of this function goes here
    %Version: 1.0 written by jbb0523 @2016-07-21 
    %Reference:Chen S S, Donoho D L, Saunders M A. Atomic decomposition by
    %basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at: 
    %http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)
    %   Detailed explanation goes here
    %   s = Phi * alpha (alpha is a sparse vector)  
    %   Given s & Phi, try to derive alpha
        [s_rows,s_columns] = size(s);  
        if s_rows<s_columns  
            s = s';%s should be a column vector  
        end 
        p = size(Phi,2);
        %according to section 3.1 of the reference
        c = ones(2*p,1);
        A = [Phi,-Phi];
        b = s;
        lb = zeros(2*p,1);
        x0 = linprog(c,[],[],A,b,lb);
        alpha = x0(1:p) - x0(p+1:2*p);
    end

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

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

    %压缩感知重构算法测试  
    clear all;close all;clc;  
    M = 64;%观测值个数  
    N = 256;%信号x的长度  
    K = 10;%信号x的稀疏度  
    Index_K = randperm(N);  
    x = zeros(N,1);  
    x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的  
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta  
    Phi = randn(M,N);%测量矩阵为高斯矩阵  
    A = Phi * Psi;%传感矩阵  
    y = Phi * x;%得到观测向量y  
    %% 恢复重构信号x  
    tic  
    theta = BP_linprog(y,A);  
    x_r = Psi * theta;% x=Psi * theta  
    toc  
    %% 绘图  
    figure;  
    plot(x_r,'k.-');%绘出x的恢复信号  
    hold on;  
    plot(x,'r');%绘出原信号x  
    hold off;  
    legend('Recovery','Original')  
    fprintf('\n恢复残差:');  
    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

    展开全文
  • 为提高一维信号去除噪声的稀疏分解基追踪算法的效率,提出了采用修正的拟牛顿法来解决基追踪去噪过程中的无约束优化问题。该算法在传统拟牛顿法的基础上,对BFGS(Broyden-Fletcher-Goldfarb-Shanno)公式进行修正,...
  • CS的基础追踪算法

    2018-04-19 10:49:54
    除匹配追踪类贪婪迭代算法之外,压缩感知重构算法另一大类就是凸优化算法或最优化逼近方法,这类方法通过将非凸问题转化为凸问题求解找到信号的逼近,其中最常用的方法就是基追踪(Basis Pursuit, BP),该方法提出...
  • 原创 压缩感知重构算法基追踪(Basis Pursuit, BP) ...

    基追踪(basis pursuit)算法是一种用来求解未知参量L1范数最小化的等式约束问题的算法。
    基追踪是通常在信号处理中使用的一种对已知系数稀疏化的手段。将优化问题中的L0范数转化为L1范数的求解就是基追踪的基本思想。
    比如我原先有一个优化问题:
    min ||x||_0(就是L0范数的最小值)subject to y=Ax。
    这个||x||_0,就是表示x中有多少个非零元素;那么我们要求min ||x||_0,就是想知道含有最多0元素的那个解x是什么。
    但是呢,L0范数有非凸性,不怎么好求解,这时我们就转而求解L1范数的优化问题。
    那么,基追踪算法就是转而求解
    min||x||_1(就是L1范数的最小值)subject to||y-Ax||_2=0(2范数)
    这个||x||_1,就是x的绝对值;那么我们要求min||x||_1,就是求绝对值最小的那个解x是什么。
    更通俗一点来讲,比如我要求一个线性方程组
    Ax=b
    x就是我们要求的未知量。这个A矩阵不是个方阵,是个欠定矩阵,那么就导致这个线性方程组会有若干组解。那么我们到底要哪组解好呢?
    如果在一般情况下,可以直接用最小二乘法来获得一组最小二乘解,就是x=(A’A)^(-1)A’b。但是我们现在利用基追踪,就是想要来获得一组含0元素最多的解。
    那么我们为什么希望我们获得的解里面0元素越多越好呢?这就要谈到“稀疏化”了。所谓稀疏化,就是希望我获得的这个解放眼望去全是0,非0元素稀稀疏疏的。这样在大样本或者高维数的情况下,计算速度不会太慢,也不会太占计算机的内存。当然,所谓稀疏解是有一定精度误差的,想要提高计算速度,必然会损失一点精度,这是不可避免的。

            除匹配追踪类贪婪迭代算法之外,压缩感知重构算法另一大类就是凸优化算法或最优化逼近方法,这类方法通过将非凸问题转化为凸问题求解找到信号的逼近,其中最常用的方法就是基追踪(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)

    function [ alpha ] = BP_linprog( s,Phi )
    %BP_linprog(Basis Pursuit with linprog) Summary of this function goes here
    %Version: 1.0 written by jbb0523 @2016-07-21 
    %Reference:Chen S S, Donoho D L, Saunders M A. Atomic decomposition by
    %basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at: 
    %http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)
    %   Detailed explanation goes here
    %   s = Phi * alpha (alpha is a sparse vector)  
    %   Given s & Phi, try to derive alpha
        [s_rows,s_columns] = size(s);  
        if s_rows<s_columns  
            s = s';%s should be a column vector  
        end 
        p = size(Phi,2);
        %according to section 3.1 of the reference
        c = ones(2*p,1);
        A = [Phi,-Phi];
        b = s;
        lb = zeros(2*p,1);
        x0 = linprog(c,[],[],A,b,lb);
        alpha = x0(1:p) - x0(p+1:2*p);
    end

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

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

    %压缩感知重构算法测试  
    clear all;close all;clc;  
    M = 64;%观测值个数  
    N = 256;%信号x的长度  
    K = 10;%信号x的稀疏度  
    Index_K = randperm(N);  
    x = zeros(N,1);  
    x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的,且位置是随机的  
    Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta  
    Phi = randn(M,N);%测量矩阵为高斯矩阵  
    A = Phi * Psi;%传感矩阵  
    y = Phi * x;%得到观测向量y  
    %% 恢复重构信号x  
    tic  
    theta = BP_linprog(y,A);  
    x_r = Psi * theta;% x=Psi * theta  
    toc  
    %% 绘图  
    figure;  
    plot(x_r,'k.-');%绘出x的恢复信号  
    hold on;  
    plot(x,'r');%绘出原信号x  
    hold off;  
    legend('Recovery','Original')  
    fprintf('\n恢复残差:');  
    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

    文章最后发布于: 2016-07-21 21:03:19
    展开全文
  • 题目:压缩感知重构算法基追踪降噪(Basis Pursuit De-Noising, BPDN)  本篇来探讨基追踪降噪(Basis Pursuit De-Noising, BPDN),与基追踪不同之处在于,基追踪降噪在模型中考虑到了噪声的存在,而这在实际中是...

    题目:压缩感知重构算法之基追踪降噪(Basis Pursuit De-Noising, BPDN)

            本篇来探讨基追踪降噪(Basis Pursuit De-Noising, BPDN),与基追踪不同之处在于,基追踪降噪在模型中考虑到了噪声的存在,而这在实际中是非常有意义的。因为考虑到了噪声,所以不同于BP的最优化模型可以转化为线性规划问题,BPDN的最优化模型可以转化为二次规划问题。

    1、基追踪降噪的提出

            文献【1】为了使基追踪能够适应有噪声的数据,提出了基追踪降噪:

            然后将基追踪降噪转化为扰动线性规划(perturbed linear program)问题:

            基追踪降噪最优化模型中的参数λ可按如下规划选取:

            并没有查到简单的MATLAB函数求解扰动线性规划问题。作者发布了文献【1】中的算法MATLAB源代码工具箱Atomizer(下载链接:http://sparselab.stanford.edu/atomizer/),并有一个说明文档,在工具箱\Documentation目录下(AboutAtomizer1208.ps,若想打开需要转换一下格式,不成功者可以参考链接http://download.csdn.net/detail/jbb0523/9583505)

    2、利用l1-magic工具箱实现基追踪降噪

            在l1-magic工具箱(链接:http://users.ece.gatech.edu/~justin/l1magic/)解决的七类问题中,并没有出现文献【1】式(5.1)类型,但在文献【2】中提到:


            即问题(3)与问题(1)等价。下面是l1-magic的(P2)问题:

            可以看出文献【2】的问题(3)与l1-magic中的问题(P2)等价,因此,我们完全可以使用l1-magic工具箱中求解(P2)的函数l1qc_logbarrier.m解决基追踪降噪(使用时注意,函数l1qc_logbarrier还调用了工具箱其它的函数)。

    3、将基追踪降噪转化为二次规划问题

            文献【2】给出了将基追踪降噪转化为二次规划问题的过程:


            对于线性代数和矩阵分析基础不好的人,这个过程可能并不是很容易看明白,现详细推导如下:

    Step0:

    与基追踪推导类似,将变量x变为两个非负变量uv之差:

    Step1:

     ,其中 

    Step2:

    Step3:

    因为是一个数(或者说是1×1的矩阵)

    所以

    Step4:

    因为

    所以

    Step5:

    Step6:

    综合Step3、Step4、Step5的结果,Step2可化为:

    Step7:

    与基追踪推导类似:

    其中, 表示转置

    Step8:

    综合Step6与Step7的结果,Step0可化为:

    其中

    Step9:

    因为是一个与变量无关的常数,因此不影响最优化极值点的位置

    所以Step8可写为标准的二次规划形式:

            求得最优化解z0后可得变量x的最优化解x0=z0(1:n)-z0(n+1:2n)。

    4、基于quadprog的基追踪降噪MATLAB代码(BPDN_quadprog.m)

            根据以上推导结果,可以写出基于MATLAB自带二次规划函数quadprog的代码:

    [plain] view plain copy
     在CODE上查看代码片派生到我的代码片
    1. function [ x ] = BPDN_quadprog( y,A,tao )  
    2. %BPDN_quadprog(Basis Pursuit DeNoising with quadprog) Summary of this function goes here  
    3. %Version: 1.0 written by jbb0523 @2016-07-22   
    4. %Reference:Nowak R D, Wright S J. Gradient projection for sparse reconstruction:  
    5. %Application to compressed sensing and other inverse problems[J].   
    6. %IEEE Journal of selected topics in signal processing, 2007, 1(4): 586-597.  
    7. %(Available at: http://pages.cs.wisc.edu/~swright/papers/FigNW07a.pdf)  
    8.     [y_rows,y_columns] = size(y);    
    9.     if y_rows<y_columns    
    10.         y = y';%y should be a column vector    
    11.     end   
    12.     n = size(A,2);  
    13.     %according to section II-A of the reference  
    14.     b = A'*y;  
    15.     c = tao*ones(2*n,1)+[-b;b];  
    16.     B = [A'*A,-A'*A;-A'*A,A'*A];  
    17.     lb = zeros(2*n,1);  
    18.     z0 = quadprog(B,c,[],[],[],[],lb);  
    19.     x = z0(1:n) - z0(n+1:2*n);  
    20. end  

    5、基追踪降噪单次重构测试代码

            测试代码与OMP测试代码有改动,首先是测量矩阵Phi使用函数orth进行了正交化,其次是观测值y加入了噪声e,可以用比较软件Beyond Compare具体比较不同之处: 

    [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. Phi = orth(Phi')';  
    12. A = Phi * Psi;%传感矩阵    
    13. sigma = 0.005;  
    14. e = sigma*randn(M,1);  
    15. y = Phi * x + e;%得到观测向量y    
    16. %% 恢复重构信号x    
    17. tic    
    18. lamda = sigma*sqrt(2*log(N));  
    19. theta =  BPDN_quadprog(y,A,lamda);  
    20. x_r = Psi * theta;% x=Psi * theta    
    21. toc    
    22. %% 绘图    
    23. figure;    
    24. plot(x_r,'k.-');%绘出x的恢复信号    
    25. hold on;    
    26. plot(x,'r');%绘出原信号x    
    27. hold off;    
    28. legend('Recovery','Original')    
    29. fprintf('\n恢复残差:');    
    30. norm(x_r-x)%恢复残差   

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

             1)图:

            2)Command Windows

    Optimization terminated: relative functionvalue changing by less

     thansqrt(OPTIONS.TolFun), no negative curvature detected in current

     trust region model and the rate of progress(change in f(x)) is slow.

    Elapsed time is 8.706162 seconds.

     恢复残差:

    ans =

       0.2723

    6、结束语

            从Command Windows输出可以看到,本次运行花费了8.706162秒,误差为0.2723。相比于之前的单次重构测试,这个运行时长非常长,误差非常大。为什么呢?

            有关运行时间长,是不是由于该二次规划变量只有一个非负约束,故可行域范围比较大,所以时间会很长;而误差大也是可以理解的,因为本来该方法的观测值就是含有噪声的,但问题也就来了,这不是基追踪“降噪”么?为什么没把“噪声”给“降”掉呢?

    7、参考文献:

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

     2Nowak R D, Wright SJ. Gradientprojection for sparse reconstruction: Application to compressed sensing andother inverse problems[J]. IEEE Journal of selected topics in signalprocessing, 2007, 1(4): 586-597. (Available at: http://pages.cs.wisc.edu/~swright/papers/FigNW07a.pdf)

    展开全文
  • 使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)

    万次阅读 热门讨论 2016-07-24 16:31:35
     基追踪(Basis Pursuit, BP)和基追踪降噪(Basis PursuitDe-Noising, BPDN)都不能称为是具体的算法,实际上分别是求解一个优化问题: 基追踪基追踪降噪:  BP和BPDN可以由多种优化求解方法去解。为了求解...

    题目:使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)

            基追踪(Basis Pursuit, BP)和基追踪降噪(Basis PursuitDe-Noising, BPDN)都不能称为是具体的算法,实际上分别是求解一个优化问题:

    基追踪:

    基追踪降噪:

            BP和BPDN可以由多种优化求解方法去解。为了求解这两个优化问题,基追踪可以转化为一个线性规划问题(参见《压缩感知重构算法之基追踪(Basis Pursuit, BP)》,http://blog.csdn.net/jbb0523/article/details/51986554),可以使用MATLAB自带的linprog函数求解;而基追踪降噪可以转化为一个二次规划问题(参见《压缩感知重构算法之基追踪降噪(Basis Pursuit De-Noising, BPDN)》,http://blog.csdn.net/jbb0523/article/details/52013669),可以使用MATLAB自带的quadprog函数求解。

            当然也可以用其它方法求解基追踪和基追踪降噪这两个问题。实际上同一个算法,不同的实现方式,运行效率也会有差异。

            由Emmanuel Candès and JustinRomberg等人开发的l1-magic工具箱(主页链接:http://users.ece.gatech.edu/justin/l1magic/)可以求解七类优化问题,其中(P1)问题即为基追踪,(P2)问题等价于基追踪降噪,详情可参见工具箱的说明文档[2]

    (P1)问题:

    (P2)问题:

    1、(P1)问题

            求解(P1)问题的函数是l1eq_pd.m,位于工具箱\l1magic\Optimization目录下。该函数的使用例程参见l1eq_example.m,位于工具箱根目录\l1magic下。在例程l1eq_example中,函数l1eq_pd共有两种调用方式,如下图所示:

            默认情况是第一种,第二种为large scale(大规模),需要注意的是,当使用largescale这种调用方式时,需要额外使用工具箱中的cgsolve.m函数,位于工具箱\l1magic\Optimization目录下。若使用第二种方式,只需将第50行到54行解除注释即可(参见附录代码)。

    2、(P2)问题

            求解(P2)问题的函数是l1qc_logbarrier.m,位于工具箱\l1magic\Optimization目录下。该函数的使用例程参见l1qc_example.m,位于工具箱根目录\l1magic下。需要注意的是,该函数需要调用l1qc_newton.m函数(位于工具箱\l1magic\Optimization目录下),在例程l1qc_example中,函数l1qc_logbarrier.m共有两种调用方式,如下图所示:

            默认情况是第一种,第二种为large scale(大规模),需要注意的是,当使用largescale这种调用方式时,也需要额外使用工具箱中的cgsolve.m函数。若使用第二种方式,只需将第49行到53行解除注释即可(参见附录代码)。

            先简单了解这么一点即可,若实际需要,再详细琢磨一下。

    3、总结

            总结一下,对于(P1)问题,相关文件有三个:l1eq_pd.m(求解函数)、cgsolve.m(大规模方式时l1eq_pd需调用此文件)、l1eq_example.m(使用范例);对于(P2)问题,相关文件有四个:l1qc_logbarrier.m(求解函数)、l1qc_newton.m(l1qc_logbarrier需要调用此文件)、cgsolve.m(大规模方式时l1qc_logbarrier需调用此文件)、l1qc_example.m(使用范例)。

            值得注意的是,函数l1qc_logbarrier求解基追踪降噪问题的效率和重构效果均要优于本人基于MATLAB自带的quadprog函数编写的BPDN_quadprog(参见《压缩感知重构算法之基追踪降噪(BasisPursuit De-Noising, BPDN)》,http://blog.csdn.net/jbb0523/article/details/52013669),具体原因没有深究。

    4、参考文献

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

    【2】Emmanuel Candès and Justin Romberg.l1-magic : Recovery of Sparse Signals via Convex Programming[EB/OL].http://users.ece.gatech.edu/justin/l1magic/downloads/l1magic.pdf

    5、附录

            作为一种备份,也作为一篇文档的完整性,这里将(P1)(P2)所涉及的MATLAB文件附在此,如果不愿意自己去下载,可以将这些代码保存为MATLAB文件即可。

    (1)第一个是两个问题在large scale方式下都要调用的cgsolve.m

    % cgsolve.m
    %
    % Solve a symmetric positive definite system Ax = b via conjugate gradients.
    %
    % Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
    %
    % A - Either an NxN matrix, or a function handle.
    %
    % b - N vector
    %
    % tol - Desired precision.  Algorithm terminates when 
    %    norm(Ax-b)/norm(b) < tol .
    %
    % maxiter - Maximum number of iterations.
    %
    % verbose - If 0, do not print out progress messages.
    %    If and integer greater than 0, print out progress every 'verbose' iters.
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
    
    if (nargin < 5), verbose = 1; end
    
    implicit = isa(A,'function_handle');
    
    x = zeros(length(b),1);
    r = b;
    d = r;
    delta = r'*r;
    delta0 = b'*b;
    numiter = 0;
    bestx = x;
    bestres = sqrt(delta/delta0); 
    while ((numiter < maxiter) && (delta > tol^2*delta0))
    
      % q = A*d
      if (implicit), q = A(d);  else  q = A*d;  end
     
      alpha = delta/(d'*q);
      x = x + alpha*d;
      
      if (mod(numiter+1,50) == 0)
        % r = b - Aux*x
        if (implicit), r = b - A(x);  else  r = b - A*x;  end
      else
        r = r - alpha*q;
      end
      
      deltaold = delta;
      delta = r'*r;
      beta = delta/deltaold;
      d = r + beta*d;
      numiter = numiter + 1;
      if (sqrt(delta/delta0) < bestres)
        bestx = x;
        bestres = sqrt(delta/delta0);
      end    
      
      if ((verbose) && (mod(numiter,verbose)==0))
        disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...
          numiter, bestres, sqrt(delta/delta0)));
      end
      
    end
    
    if (verbose)
      disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
    end
    x = bestx;
    res = bestres;
    iter = numiter;
    

    (2)第二个是(P1)问题求解函数l1eq_pd.m

    % l1eq_pd.m
    %
    % Solve
    % min_x ||x||_1  s.t.  Ax = b
    %
    % Recast as linear program
    % min_{x,u} sum(u)  s.t.  -u <= x <= u,  Ax=b
    % and use primal-dual interior point method
    %
    % Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
    %
    % x0 - Nx1 vector, initial point.
    %
    % A - Either a handle to a function that takes a N vector and returns a K 
    %     vector , or a KxN matrix.  If A is a function handle, the algorithm
    %     operates in "largescale" mode, solving the Newton systems via the
    %     Conjugate Gradients algorithm.
    %
    % At - Handle to a function that takes a K vector and returns an N vector.
    %      If A is a KxN matrix, At is ignored.
    %
    % b - Kx1 vector of observations.
    %
    % pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
    %     the duality gap is less than pdtol).  
    %     Default = 1e-3.
    %
    % pdmaxiter - Maximum number of primal-dual iterations.  
    %     Default = 50.
    %
    % cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
    %     Default = 1e-8.
    %
    % cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
    %     if A is a matrix.
    %     Default = 200.
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
    
    largescale = isa(A,'function_handle');
    
    if (nargin < 5), pdtol = 1e-3;  end
    if (nargin < 6), pdmaxiter = 50;  end
    if (nargin < 7), cgtol = 1e-8;  end
    if (nargin < 8), cgmaxiter = 200;  end
    
    N = length(x0);
    
    alpha = 0.01;
    beta = 0.5;
    mu = 10;
    
    gradf0 = [zeros(N,1); ones(N,1)];
    
    % starting point --- make sure that it is feasible
    if (largescale)
      if (norm(A(x0)-b)/norm(b) > cgtol)
        disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
        AAt = @(z) A(At(z));
        [w, cgres, cgiter] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
        if (cgres > 1/2)
          disp('A*At is ill-conditioned: cannot find starting point');
          xp = x0;
          return;
        end
        x0 = At(w);
      end
    else
      if (norm(A*x0-b)/norm(b) > cgtol)
        disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
        opts.POSDEF = true; opts.SYM = true;
        [w, hcond] = linsolve(A*A', b, opts);
        if (hcond < 1e-14)
          disp('A*At is ill-conditioned: cannot find starting point');
          xp = x0;
          return;
        end
        x0 = A'*w;
      end  
    end
    x = x0;
    u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
    
    % set up for the first iteration
    fu1 = x - u;
    fu2 = -x - u;
    lamu1 = -1./fu1;
    lamu2 = -1./fu2;
    if (largescale)
      v = -A(lamu1-lamu2);
      Atv = At(v);
      rpri = A(x) - b;
    else
      v = -A*(lamu1-lamu2);
      Atv = A'*v;
      rpri = A*x - b;
    end
    
    sdg = -(fu1'*lamu1 + fu2'*lamu2);
    tau = mu*2*N/sdg;
    
    rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
    rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
    resnorm = norm([rdual; rcent; rpri]);
    
    pditer = 0;
    done = (sdg < pdtol) | (pditer >= pdmaxiter);
    while (~done)
      
      pditer = pditer + 1;
      
      w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;
      w2 = -1 - 1/tau*(1./fu1 + 1./fu2);
      w3 = -rpri;
      
      sig1 = -lamu1./fu1 - lamu2./fu2;
      sig2 = lamu1./fu1 - lamu2./fu2;
      sigx = sig1 - sig2.^2./sig1;
      
      if (largescale)
        w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));
        h11pfun = @(z) -A(1./sigx.*At(z));
        [dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
        if (cgres > 1/2)
          disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
          xp = x;
          return
        end
        dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;
        Adx = A(dx);
        Atdv = At(dv);
      else
        w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));
        H11p = A*(sparse(diag(1./sigx))*A');
        opts.POSDEF = true; opts.SYM = true;
        [dv,hcond] = linsolve(H11p, w1p, opts);
        if (hcond < 1e-14)
          disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
          xp = x;
          return
        end
        dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;
        Adx = A*dx;
        Atdv = A'*dv;
      end
      
      du = (w2 - sig2.*dx)./sig1;
      
      dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;
      dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2;
      
      % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
      indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);
      s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);
      indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);
      s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);
      
      % backtracking line search
      suffdec = 0;
      backiter = 0;
      while (~suffdec)
        xp = x + s*dx;  up = u + s*du; 
        vp = v + s*dv;  Atvp = Atv + s*Atdv; 
        lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
        fu1p = xp - up;  fu2p = -xp - up;  
        rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];
        rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
        rpp = rpri + s*Adx;
        suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);
        s = beta*s;
        backiter = backiter + 1;
        if (backiter > 32)
          disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
          xp = x;
          return
        end
      end
      
      
      % next iteration
      x = xp;  u = up;
      v = vp;  Atv = Atvp; 
      lamu1 = lamu1p;  lamu2 = lamu2p;
      fu1 = fu1p;  fu2 = fu2p;
      
      % surrogate duality gap
      sdg = -(fu1'*lamu1 + fu2'*lamu2);
      tau = mu*2*N/sdg;
      rpri = rpp;
      rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
      rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
      resnorm = norm([rdual; rcent; rpri]);
      
      done = (sdg < pdtol) | (pditer >= pdmaxiter);
      
      disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',...
        pditer, tau, sum(u), sdg, norm(rdual), norm(rpri)));
      if (largescale)
        disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
      else
        disp(sprintf('                  H11p condition number = %8.3e', hcond));
      end
      
    end
    

    (3)第三个是(P1)问题使用范例l1eq_example.m(处于large scale方式下)

    % l1eq_example.m
    %
    % Test out l1eq code (l1 minimization with equality constraints).
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    % put key subdirectories in path if not already there
    % path(path, './Optimization');
    % path(path, './Data');
    
    % To reproduce the example in the documentation, uncomment the 
    % two lines below
    %load RandomStates
    %rand('state', rand_state);
    %randn('state', randn_state);
    
    % signal length
    N = 512;
    % number of spikes in the signal
    T = 20;
    % number of observations to make
    K = 120;
    
    % random +/- 1 signal
    x = zeros(N,1);
    q = randperm(N);
    x(q(1:T)) = sign(randn(T,1));
    
    % measurement matrix
    disp('Creating measurment matrix...');
    A = randn(K,N);
    A = orth(A')';
    disp('Done.');
    	
    % observations
    y = A*x;
    
    % initial guess = min energy
    x0 = A'*y;
    
    % solve the LP
    % tic
    % xp = l1eq_pd(x0, A, [], y, 1e-3);
    % toc
    
    % large scale
    Afun = @(z) A*z;
    Atfun = @(z) A'*z;
    tic
    xp = l1eq_pd(x0, Afun, Atfun, y, 1e-3, 30, 1e-8, 200);
    toc
    

    (4)第四个是(P2)问题求解函数l1qc_logbarrier.m

    % l1qc_logbarrier.m
    %
    % Solve quadratically constrained l1 minimization:
    % min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
    %
    % Reformulate as the second-order cone program
    % min_{x,u}  sum(u)   s.t.    x - u <= 0,
    %                            -x - u <= 0,
    %      1/2(||Ax-b||^2 - \epsilon^2) <= 0
    % and use a log barrier algorithm.
    %
    % Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
    %
    % x0 - Nx1 vector, initial point.
    %
    % A - Either a handle to a function that takes a N vector and returns a K 
    %     vector , or a KxN matrix.  If A is a function handle, the algorithm
    %     operates in "largescale" mode, solving the Newton systems via the
    %     Conjugate Gradients algorithm.
    %
    % At - Handle to a function that takes a K vector and returns an N vector.
    %      If A is a KxN matrix, At is ignored.
    %
    % b - Kx1 vector of observations.
    %
    % epsilon - scalar, constraint relaxation parameter
    %
    % lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
    %         Also, the number of log barrier iterations is completely
    %         determined by lbtol.
    %         Default = 1e-3.
    %
    % mu - Factor by which to increase the barrier constant at each iteration.
    %      Default = 10.
    %
    % cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
    %     Default = 1e-8.
    %
    % cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
    %     if A is a matrix.
    %     Default = 200.
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  
    
    largescale = isa(A,'function_handle');
    
    if (nargin < 6), lbtol = 1e-3; end
    if (nargin < 7), mu = 10; end
    if (nargin < 8), cgtol = 1e-8; end
    if (nargin < 9), cgmaxiter = 200; end
    
    newtontol = lbtol;
    newtonmaxiter = 50;
    
    N = length(x0);
    
    % starting point --- make sure that it is feasible
    if (largescale)
      if (norm(A(x0)-b) > epsilon)
        disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
        AAt = @(z) A(At(z));
        [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
        if (cgres > 1/2)
          disp('A*At is ill-conditioned: cannot find starting point');
          xp = x0;
          return;
        end
        x0 = At(w);
      end
    else
      if (norm(A*x0-b) > epsilon)
        disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
        opts.POSDEF = true; opts.SYM = true;
        [w, hcond] = linsolve(A*A', b, opts);
        if (hcond < 1e-14)
          disp('A*At is ill-conditioned: cannot find starting point');
          xp = x0;
          return;
        end
        x0 = A'*w;
      end  
    end
    x = x0;
    u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
    
    disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u)));
    
    % choose initial value of tau so that the duality gap after the first
    % step will be about the origial norm
    tau = max((2*N+1)/sum(abs(x0)), 1);
                                                                                                                              
    lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu));
    disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
    
    totaliter = 0;
    
    for ii = 1:lbiter
    
      [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
      totaliter = totaliter + ntiter;
      
      disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
        ii, sum(abs(xp)), sum(up), tau, totaliter));
      
      x = xp;
      u = up;
     
      tau = mu*tau;
      
    end
    

    (5)第五个是(P2)问题求解函数须调用的l1qc_newton.m

    % l1qc_newton.m
    %
    % Newton algorithm for log-barrier subproblems for l1 minimization
    % with quadratic constraints.
    %
    % Usage: 
    % [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, 
    %                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
    %
    % x0,u0 - starting points
    %
    % A - Either a handle to a function that takes a N vector and returns a K 
    %     vector , or a KxN matrix.  If A is a function handle, the algorithm
    %     operates in "largescale" mode, solving the Newton systems via the
    %     Conjugate Gradients algorithm.
    %
    % At - Handle to a function that takes a K vector and returns an N vector.
    %      If A is a KxN matrix, At is ignored.
    %
    % b - Kx1 vector of observations.
    %
    % epsilon - scalar, constraint relaxation parameter
    %
    % tau - Log barrier parameter.
    %
    % newtontol - Terminate when the Newton decrement is <= newtontol.
    %         Default = 1e-3.
    %
    % newtonmaxiter - Maximum number of iterations.
    %         Default = 50.
    %
    % cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
    %     Default = 1e-8.
    %
    % cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
    %     if A is a matrix.
    %     Default = 200.
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    
    function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 
    
    % check if the matrix A is implicit or explicit
    largescale = isa(A,'function_handle');
    
    % line search parameters
    alpha = 0.01;
    beta = 0.5;  
    
    if (~largescale), AtA = A'*A; end
    
    % initial point
    x = x0;
    u = u0;
    if (largescale), r = A(x) - b; else  r = A*x - b; end
    fu1 = x - u;
    fu2 = -x - u;
    fe = 1/2*(r'*r - epsilon^2);
    f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe));
    
    niter = 0;
    done = 0;
    while (~done)
      
      if (largescale), atr = At(r); else  atr = A'*r; end
      
      ntgz = 1./fu1 - 1./fu2 + 1/fe*atr;
      ntgu = -tau - 1./fu1 - 1./fu2;
      gradf = -(1/tau)*[ntgz; ntgu];
      
      sig11 = 1./fu1.^2 + 1./fu2.^2;
      sig12 = -1./fu1.^2 + 1./fu2.^2;
      sigx = sig11 - sig12.^2./sig11;
        
      w1p = ntgz - sig12./sig11.*ntgu;
      if (largescale)
        h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr;
        [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
        if (cgres > 1/2)
          disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
          xp = x;  up = u;
          return
        end
        Adx = A(dx);
      else
        H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr';
        opts.POSDEF = true; opts.SYM = true;
        [dx,hcond] = linsolve(H11p, w1p, opts);
        if (hcond < 1e-14)
          disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
          xp = x;  up = u;
          return
        end
        Adx = A*dx;
      end
      du = (1./sig11).*ntgu - (sig12./sig11).*dx;  
     
      % minimum step size that stays in the interior
      ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0);
      aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;
      smax = min(1,min([...
        -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...
        (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)
        ]));
      s = (0.99)*smax;
      
      % backtracking line search
      suffdec = 0;
      backiter = 0;
      while (~suffdec)
        xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
        fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);
        fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));
        flin = f + alpha*s*(gradf'*[dx; du]);
        suffdec = (fp <= flin);
        s = beta*s;
        backiter = backiter + 1;
        if (backiter > 32)
          disp('Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)');
          xp = x;  up = u;
          return
        end
      end
      
      % set up for next iteration
      x = xp; u = up;  r = rp;
      fu1 = fu1p;  fu2 = fu2p;  fe = fep;  f = fp;
      
      lambda2 = -(gradf'*[dx; du]);
      stepsize = s*norm([dx; du]);
      niter = niter + 1;
      done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
      
      disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
        niter, f, lambda2/2, stepsize));
      if (largescale)
        disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
      else
        disp(sprintf('                  H11p condition number = %8.3e', hcond));
      end
          
    end
    

    (6)第六个是(P2)问题使用范例l1qc_example.m(处于large scale方式下)

    % l1qc_example.m
    %
    % Test out l1qc code (l1 minimization with quadratic constraint).
    %
    % Written by: Justin Romberg, Caltech
    % Email: jrom@acm.caltech.edu
    % Created: October 2005
    %
    
    % put optimization code in path if not already there
    path(path, './Optimization');
    
    % signal length
    N = 512;
    
    % number of spikes to put down
    T = 20;
    
    % number of observations to make
    K = 120;
    
    % random +/- 1 signal
    x = zeros(N,1);
    q = randperm(N);
    x(q(1:T)) = sign(randn(T,1));
    
    % measurement matrix
    disp('Creating measurment matrix...');
    A = randn(K,N);
    A = orth(A')';
    disp('Done.');
    	
    % noisy observations
    sigma = 0.005;
    e = sigma*randn(K,1);
    y = A*x + e;
    
    % initial guess = min energy
    x0 = A'*y;
    
    % take epsilon a little bigger than sigma*sqrt(K)
    epsilon =  sigma*sqrt(K)*sqrt(1 + 2*sqrt(2)/sqrt(K));
                                                                                                                               
    % tic
    % xp = l1qc_logbarrier(x0, A, [], y, epsilon, 1e-3);
    % toc
    
    % large scale
    Afun = @(z) A*z;
    Atfun = @(z) A'*z;
    tic
    xp = l1qc_logbarrier(x0, Afun, Atfun, y, epsilon, 1e-3, 50, 1e-8, 500);
    toc
    


    展开全文
  • 摘要:提出了一种新的基追踪求解算法。依据信号特性自适应地选取字典;通过l1范数的近似表示,将有约束的极值问题转化为无约束问题,并利用一种新的迭代算法进行快速求解;几类典型信号实验结果验证了本方法具有良好...
  • 凸优化之基追踪

    2018-01-23 09:32:00
    关于基追踪(BP),在压缩感知重构中我们所待求解的问题是L0范数问题,因为L1范数与L0范数等价,所以将L0范数转换为L1范数问题来求解,基追踪是将L1范数问题转为成为线性规划问题来进行求解,博主还提到了基追踪降噪...
  • 并基于新的模型, 提出了一种改进的压缩采样匹配追踪算法。该算法通过构造一个感知测量矩阵, 在信号替代阶段中取代随机测量矩阵来减少相关性对支撑集筛选的影响, 最后可在乘性噪声存在的情况下实现了信号的精确重建...
  • 压缩感知理论是一种充分利用信号稀疏性或者可压缩性的伞新的信号采样理论。该理论表明,通过采集少 ...于正则化的自适应匹配追踪算法(Regularized Adaptive Matching Pursuit,RAMP)用于压缩感知信号的重建。
  • 大规模MIMO系统中联合扰动正交匹配追踪算法的稀疏信道估计,张晋鹏,徐文波,压缩感知技术可以从已知的中重构稀疏信号。目前这项技术已经广泛应用于信号处理和图像压缩等领域中。在多用户大规模MIMO通信系��
  • 针对传统的l2-范数信道估计精度低的问题, 提出了一种基于基追踪去噪(BPDN)的水声正交频分复用稀疏信道估计方法, 该方法针对水声信道的稀疏特性, 利用少量的观测值即可以很高的精度估计出信道冲激响应. 与贪婪追踪类...

空空如也

空空如也

1 2 3 4 5
收藏数 92
精华内容 36
关键字:

基追踪算法