精华内容
下载资源
问答
  • 在希尔伯特空间,信号可以由这个空间的表征正交基线性表示,正交基稀疏表达一个信号也有很多缺点,因为一组基表达信号的能力取决于信号的特征是否基向量的特征相吻合;例如,光滑连续信号可以被傅里叶基稀疏的表达...

     

    在希尔伯特空间,信号可以由这个空间的表征正交基线性表示,正交基稀疏表达一个信号也有很多缺点,因为一组基表达信号的能力取决于信号的特征是否与基向量的特征相吻合;例如,光滑连续信号可以被傅里叶基稀疏的表达,但脉冲信号就不行;再如,带有孤立不连续点的平滑信号可以被小波基稀疏表达,但小波基在表达傅里叶频谱中有窄带高频支撑的信号时却是无效的。现实世界中的信号经常包含有用单一基所不能表达的特征,对于这些信号,你或许希望可以选择来自不同基的向量(如用小波基和傅里叶基来联合表达一个信号)。因为你想保证你可以表达一个信号空间的所有信号向量,所以由所有可选向量组成的字典应该能够张成这个信号空间。然而由于这组字典中的向量来自不同的基,它们可能不是线性独立的。由于这组字典中的向量不是线性独立的,会造成用这组字典做信号表达时系数不唯一。然而如果创建一组冗余字典,你就可以把你的信号展开在一组可以适应各种时频或时间-尺度特性的向量上。这样构造的字典可以极大地增加你稀疏表达各种特性信号的能力。

    通过构建字典D,继而求解信号在这个空间的稀疏表征,可以用来进行去等去噪,超分辨,信号复原(压缩感知)等操作,稀疏求解模型构建如下:

    0范数约束并非是凸问题,无法通过凸优化获得最优解。经典的方式是匹配追踪MP和正交匹配追踪OMP。

    1.匹配追踪MP

    资料较多此处不赘述,每一次迭代时,寻找残差与当前字典中最为匹配(内积绝对值最大)的原子。引用他人的博客如下:

    从向量分解的角度,k次残差到k+1次残差进行了向量的三角分解,且依赖x_{r_{k+1}}与第k+1残差正交。假设当前

                                                                               y=\sum_{n=1}^{k}a_{k}x_{k}+R}

    R项在每次三角分解时,最为匹配的原子可能是之前已经选中的原子,换而言之,每次的迭代并非是最优迭代,解也并非是最有解,因此需要更多次的迭代计算代价大。理论上,最优解要求残差R与当前拟合项\sum_{n=1}^{k}a_{k}x_{k}}正交,这一思想即促生正交匹配追中(OMP)。

    2.正交匹配追踪OMP

    先给出计算方式:

    现在证明如下:

    注意,公式(4)本质上是施密特正交化,残留的正交项gamma本事是对已经选中的原子进行施密特正交化的结果。继而有下列推导,其中(7)

    从上列的推导可以看出,每一次迭代中,为了满足假设残差与已迭代项正交,已经选中的原子的系数也需要更新。更加广泛的使用是如下算法:

     

    给出的代码来自港大的沙威,一维信号重建:

    %  1-D信号压缩传感的实现(正交匹配追踪法Orthogonal Matching Pursuit)
    %  测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
    %  编程人--香港大学电子工程系 沙威  Email: wsha@eee.hku.hk
    %  编程时间:2008年11月18日
    %  文档下载: http://www.eee.hku.hk/~wsha/Freecode/freecode.htm 
    %  参考文献:Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    
    clc;clear
    
    %%  1. 时域测试信号生成
    K=7;      %  稀疏度(做FFT可以看出来)
    N=256;    %  信号长度
    M=64;     %  测量数(M>=K*log(N/K),至少40,但有出错的概率)
    f1=50;    %  信号频率1
    f2=100;   %  信号频率2
    f3=200;   %  信号频率3
    f4=400;   %  信号频率4
    fs=800;   %  采样频率
    ts=1/fs;  %  采样间隔
    Ts=1:N;   %  采样序列
    x=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)+0.1*cos(2*pi*f3*Ts*ts)+0.9*cos(2*pi*f4*Ts*ts);  %  完整信号,由4个信号叠加而来
    
    %%  2.  时域信号压缩传感
    Phi=randn(M,N);                                   %  测量矩阵(高斯分布白噪声)64*256的扁矩阵,Phi也就是文中说的D矩阵
    s=Phi*x.';                                        %  获得线性测量 ,s相当于文中的y矩阵
    
    %%  3.  正交匹配追踪法重构信号(本质上是L_1范数最优化问题)
    %匹配追踪:找到一个其标记看上去与收集到的数据相关的小波;在数据中去除这个标记的所有印迹;不断重复直到我们能用小波标记“解释”收集到的所有数据。
    
    m=2*K;                                            %  算法迭代次数(m>=K),设x是K-sparse的
    Psi=fft(eye(N,N))/sqrt(N);                        %  傅里叶正变换矩阵
    T=Phi*Psi';                                       %  恢复矩阵(测量矩阵*正交反变换矩阵)
    
    hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r_n=s;                                            %  残差值
    
    for times=1:m;                                    %  迭代次数(有噪声的情况下,该迭代次数为K)
        for col=1:N;                                  %  恢复矩阵的所有列向量
            product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置,即找到一个其标记看上去与收集到的数据相关的小波
        Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充    
        
        T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零),在数据中去除这个标记的所有印迹
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
        r_n=s-Aug_t*aug_y;                            %  残差
        pos_array(times)=pos;                         %  纪录最大投影系数的位置
    end
    hat_y(pos_array)=aug_y;                           %  重构的谱域向量
    hat_x=real(Psi'*hat_y.');                         %  做逆傅里叶变换重构得到时域信号
    
    %%  4.  恢复信号和原始信号对比
    figure(1);
    hold on;
    plot(hat_x,'k.-')                                 %  重建信号
    plot(x,'r')                                       %  原始信号
    legend('Recovery','Original')
    norm(hat_x.'-x)/norm(x)                           %  重构误差

    二维图像重建:

    %  本程序实现图像LENA的压缩传感
    %  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk
    %  算法采用正交匹配法,参考文献 Joel A. Tropp and Anna C. Gilbert 
    %  Signal Recovery From Random Measurements Via Orthogonal Matching
    %  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
    %  DECEMBER 2007.
    %  该程序没有经过任何优化
    
    function Wavelet_OMP
    
    clc;clear
    
    %  读文件
    X=imread('lena256.bmp');
    X=double(X);
    [a,b]=size(X);
    
    %  小波变换矩阵生成
    ww=DWT(a);
    
    %  小波变换让图像稀疏化(注意该步骤会耗费时间,但是会增大稀疏度)
    X1=ww*sparse(X)*ww';
    % X1=X;
    X1=full(X1);
    
    %  随机矩阵生成
    M=190;
    R=randn(M,a);
    % R=mapminmax(R,0,255);
    % R=round(R);
    
    %  测量值
    Y=R*X1;
    
    %  OMP算法
    %  恢复矩阵
    X2=zeros(a,b); 
    %  按列循环
    for i=1:b 
        %  通过OMP,返回每一列信号对应的恢复值(小波域)
        rec=omp(Y(:,i),R,a);
        %  恢复值矩阵,用于反变换
        X2(:,i)=rec;
    end
    
    
    %  原始图像
    figure(1);
    imshow(uint8(X));
    title('原始图像');
    
    %  变换图像
    figure(2);
    imshow(uint8(X1));
    title('小波变换后的图像');
    
    %  压缩传感恢复的图像
    figure(3);
    %  小波反变换
    X3=ww'*sparse(X2)*ww; 
    % X3=X2;
    X3=full(X3);
    imshow(uint8(X3));
    title('恢复的图像');
    
    %  误差(PSNR)
    %  MSE误差
    errorx=sum(sum(abs(X3-X).^2));        
    %  PSNR
    psnr=10*log10(255*255/(errorx/a/b))   
    
    
    %  OMP的函数
    %  s-测量;T-观测矩阵;N-向量大小
    function hat_y=omp(s,T,N)
    Size=size(T);                                     %  观测矩阵大小
    M=Size(1);                                        %  测量
    hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
    Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
    r_n=s;                                            %  残差值
    
    for times=1:M;                                  %  迭代次数(稀疏度是测量的1/4)
        for col=1:N;                                  %  恢复矩阵的所有列向量
            product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
        end
        [val,pos]=max(product);                       %  最大投影系数对应的位置
        Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充
        T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零)
        aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
        r_n=s-Aug_t*aug_y;                            %  残差
        pos_array(times)=pos;                         %  纪录最大投影系数的位置
        
        if (norm(r_n)<0.9)                            %  残差足够小
            break;
        end
    end
    hat_y(pos_array)=aug_y;                           %  重构的向量

     

     

    参考:[1]压缩感知与稀疏表征

               [2]浅谈压缩感知(九):正交匹配追踪算法OMP

    展开全文
  • 施密特(Schimidt)正交化与正交匹配追踪

    万次阅读 多人点赞 2015-04-17 18:10:45
    上一篇《稀疏表示匹配追踪》中详细的解释了匹配追踪(Matching Pursuit,MP)的流程,在最后给出了正交匹配追踪(Orthogonal Matching Pursuit,OMP)的流程,并指出OMPMP的不同根本在于残差更新过程:OMP减去的...

    题目:施密特(Schimidt)正交化与正交匹配追踪

            文献[1]中给出了施密特(Schimidt)正交化的过程:


    上面的的[x,y]表示向量内积,[x,y]=xTy=yTx=[x,y]。施密特正交化公式中的br实际上可写为:


    分子之所以可以这么变化是由于[x,y]实际上为一个数,因此[x,y]x=x[x,y]= xxTy

            上一篇《稀疏表示与匹配追踪》中详细的解释了匹配追踪(Matching Pursuit,MP)的流程,在最后给出了正交匹配追踪(Orthogonal Matching Pursuit,OMP)的流程,并指出OMP与MP的不同根本在于残差更新过程:OMP减去的Pem是在所有被选择过的原子组成的矩阵φt所张成空间上的正交投影,而MP减去的Pem是在本次被选择的原子φp所张成空间上的正交投影。正交投影及正交投影变换矩阵的概念可以参见《压缩感知中的数学知识:投影矩阵(projectionmatrix)》。

            首先给出一个结论:

            设OMP共从冗余字典中选择了r个原子,分别是a1a2,……,ar,根据正交匹配追踪的流程可以知道待分解信号x最后剩余的残差eromp

             (式1)

            该残差也可以表示为

             (式2)

    其中矩阵A为选择的r个原子组成的矩阵,e(r-1)omp为选择(r-1)个原子时的残差。

            将选择的r个原子a1a2,……,ar进行施密特正交化:


    则残差eromp还可以写为

             (式3)

    (式1)一般出现在稀疏分解算法中,(式2)一般出现在重构算法中,(式3)是自己琢磨出来的(受到沙威的文档中提到的施密特正交化的启发,但沙威只限于向量情况下,详情可参见[3],此处相当于一个推广)。

            而实际上(式1)、(式2)和(式3)三种正交匹配追踪算法的残差计算方法是等价的,纯数学证明比较复杂,这里给一个MATLAB程序,可以验证三者的等价性:

    %% 验证OMP残差求解过程与Schmidt正交化的关系
    M = 4;N = 10;
    Phi = randn(M,N);
    for nn = 1:N
        Phi(:,nn) = Phi(:,nn)/norm(Phi(:,nn));
    end
    b = randn(M,1);
    e0 = b;%初始化残差为待稀疏信号b
    %选第1列
    c1 = Phi'*e0;%求出矩阵Phi每列与b的内积
    [val1,pos1]=max(abs(c1));%找到内积中最大的列及其内积值
    phit = [Phi(:,pos1)];%由所有选出的列组合的矩阵
    Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵
    e1ompe = e0 - Pphi*e0;%OMP用上一次残差减去残差在phit列空间的正交投影
    e1ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影
    x = Phi(:,pos1);%Schimidt正交化第一个向量
    Px = x*(x'*x)^(-1)*x';
    %实际上是b - Px*b
    e1ompsmt = e0 - Px*b;
    e1 = e1ompe;
    norm(e1ompe-e1ompb)+norm(e1ompsmt-e1ompb)
    %选第2列
    c2 = Phi'*e1;%求出矩阵Phi每列与e1的内积
    [val2,pos2]=max(abs(c2));%找到内积中最大的列及其内积值
    phit = [Phi(:,pos1) Phi(:,pos2)];%由所有选出的列组合的矩阵
    Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵
    e2ompe = e1 - Pphi*e1;%OMP用上一次残差减去残差在phit列空间的正交投影
    e2ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影
    y = Phi(:,pos2) - Px*Phi(:,pos2);%Schimidt正交化第二个向量
    Py = y*(y'*y)^(-1)*y';
    %实际上是b - Px*b - Py*b
    e2ompsmt = e1 - Py*b;%上一次残差减去b在第2列正交化所得z上的投影
    e2 = e2ompe;
    norm(e2ompe-e2ompb)+norm(e2ompsmt-e2ompb)
    %选第3列
    c3 = Phi'*e2;%求出矩阵Phi每列与e2的内积
    [val3,pos3]=max(abs(c3));%找到内积中最大的列及其内积值
    phit = [Phi(:,pos1) Phi(:,pos2) Phi(:,pos3)];%由所有选出的列组合的矩阵
    Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵
    e3ompe = e2 - Pphi*e2;%OMP用上一次残差减去残差在phit列空间的正交投影
    e3ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影
    z = Phi(:,pos3) - Px*Phi(:,pos3) - Py*Phi(:,pos3);%Schimidt正交化第三个向量
    Pz = z*(z'*z)^(-1)*z';
    %实际上是b - Px*b - Py*b - Pz*b
    e3ompsmt = e2 - Pz*b;%上一次残差减去b在第3列正交化所得z上的投影
    e3 = e3ompe;
    norm(e3ompe-e3ompb)+norm(e3ompsmt-e3ompb)

    程序最后一行的输出结果为零,说明三个残差是相等的。

            这里可以看出OMP与Schimidt正交化的关系:

            OMP分解过程,实际上是将所选原子依次进行Schimidt正交化,然后将待分解信号减去在正交化后的原子上各自的分量即可得残差。其实(式3)求残差的过程也是在进行施密特正交化

            有个关键问题还是要说的,分解后在所选择各原子上的系数是多少呢?答案其实也很简单,各个系数是(ATA)-1ATx,即最小二乘解,这个解是一个列向量,每一个元素分别是组成矩阵A的各原子的线性组合系数,这个在下一篇《正交匹配追踪(OMP)在稀疏分解与压缩感知重构中的异同》也会明确再次说明。

            同理,若设MP共从冗余字典中选择了r个原子,分别是a1a2,……,ar,根据匹配追踪的流程可以知道待分解信号x每次迭代后剩余的残差ermp

    比较式(3)的第2个等号表示的eromp与此处的ermp也可以体会出OMP与MP的区别吧。

    【参考文献】

    【1】同济大学数学系. 线性代数(第五版)[M].高等教育出版社,2007:114.

    【2】Joel A. Tropp and AnnaC. Gilbert . Signal Recovery From Random Measurements Via Orthogonal MatchingPursuit[J]. IEEE Transactions on Information Theory, VOL. 53, NO. 12, DECEMBER2007.

    【3】沙威. “压缩传感”引论.http://www.eee.hku.hk/~wsha/Freecode/Files/Compressive_Sensing.pdf

    展开全文
  • 浅析匹配追踪算法(Matching Pursuit,MP)与正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)前言第一个例子基本概念OMP算法求解过程MP算法求解过程第二个例子MP和OMP的缺点后记参考文献 前言 考虑下列数学...

    浅析匹配追踪算法(Matching Pursuit,MP)与正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)

    前言

    考虑下列数学模型:
    Alt
    在压缩感知(Compressed sensing)术语中,从X和A中找到Y称为压缩(compression)。 相反,从Y和A中找到原始X称为重构(reconstruction)
    x被称为原始信号(original signal)原始向量(original vector),A被称为压缩矩阵(Compression matrix)或感知矩阵(sensing matrix),而Y被称为压缩信号(compressed signal)或压缩向量(compressed vector)

    第一个例子

    给定压缩矩阵A和原始向量X:
    Alt
    很容易就可以得到压缩向量Y:
    Alt
    然而,反过来给定压缩向量Y和原始向量X:
    [Alt]
    很难得到找到原始的(或尽可能接近的) X。

    基本概念

    首先介绍一些基本概念。

    我们将压缩矩阵A看作是列向量的集合,即:
    Alt
    其中:
    在这里插入图片描述
    这些列向量(b1,b2,b3)被称为基(basis)或原子(atoms)
    令:
    在这里插入图片描述
    那么:
    在这里插入图片描述
    上述方程告诉我们,Y是用X中的系数(x1, x2, x3)表示的基(b1,b2,b3)的线性组合,且已知x1=-1.2,x2=1,x3=0。换句话说,在上述方程中,基b1对生成Y的贡献分别是-1.2,基b2对生成Y的贡献分别是1,基b3对生成Y的贡献分别是0。

    OMP算法求解过程

    OMP的工作原理与MP相似,但MP是找出哪一个基对Y贡献最大,然后哪一个基贡献第二,哪一个基贡献第三,依此类推。OMP计算过程需要进行N次迭代,其中N是A中基的个数。在上述的示例中,N为3。

    1. 计算基的贡献度
    上述例子中有三个基,分别是:
    在这里插入图片描述
    以及
    在这里插入图片描述
    所以Y中每个基的贡献度分别是:
    在这里插入图片描述
    显然,基b1的贡献最大,其贡献度为-1.34(使用绝对值)。
    在第一轮迭代中,选择基b1,该基的贡献度为-1.34。
    当然,可以在一个单独的步骤中计算每一个基的点积,即:
    在这里插入图片描述
    2. 计算残差r1
    首先,选择第一个基,即基b1和贡献度λ1=-1.34。残差为:
    在这里插入图片描述
    如图2所示,残差r1与第一个基b1是垂直的。
    在这里插入图片描述
    3. 重复迭代
    第一次迭代后,我们得到所选基b1。把这些基重新组合成一个新的基Anew。即:
    在这里插入图片描述
    重写贡献度值,即:
    在这里插入图片描述
    值-1.34在第一个元素因为它是第一个基b1的贡献。
    第一次迭代后残差为:
    在这里插入图片描述
    现在我们要从剩下的基b2和b3中选择对残差r贡献最大的基,在一步中:
    在这里插入图片描述

    因此,基b2的贡献更大,即0.98。
    我们把选定的基b2添加到矩阵Anew中,
    在这里插入图片描述

    现在,这一步不同于MP。这里,我们需要重新计算基b1和b2共同对Y的贡献(而不是基b2对r1的贡献)

    为了解决这个贡献问题,我们必须求出最小值平方问题,很容易表述为:
    在这里插入图片描述
    在数学术语中,这个公式可以写成:
    在这里插入图片描述
    根据最小二乘原理,上述公式的解为:
    在这里插入图片描述
    其中,Anew+是Anew的广义逆矩阵。
    在这里插入图片描述
    在上述例子中
    在这里插入图片描述
    在Matlab中计算广义逆矩阵很简单,只需使用pinv()命令即可。
    在这里插入图片描述
    计算λ之后,更新xrec
    在这里插入图片描述
    在这里我们把λ放在xrec在第一和第二行,因为它对应于第一和第二被选择的基b1和b2
    重新得到Anew和λ之后,现在我们可以更新新的残差r2
    在这里插入图片描述
    到此,残差为0。我们可以停在这里,或者继续下一步。(如果在这里停止,我们可以节省一些计算工作) 。

    4. 最后迭代

    由于残差已经消失,因此不需要此步骤。(许多OMP实现软件需要信号稀疏度K的参数。这是该软件将迭代K次的线索。此后它将停止,无论剩下多少残差)
    因此,我们的重构信号为:
    在这里插入图片描述
    它和原始信号是一样的.

    以上便是OMP的整个计算流程。

    MP算法求解过程

    1. 计算贡献度
    Y中每个基的贡献度分别是:
    在这里插入图片描述
    显然,基b1的贡献最大,其贡献度为-1.34(使用绝对值)。

    ==================================================================================

    第一轮迭代

    选择基b1,该基的贡献度为-1.34
    重写贡献度值,即:
    在这里插入图片描述

    计算残差
    在这里插入图片描述
    接下来我们计算这个残差和基b2和基b3的点积。(不需要用基b1计算,因为这个残差必定垂直于基b1).
    在这里插入图片描述
    取绝对值,我们得到基b2是第二强的影响。
    贡献度值,即:
    在这里插入图片描述
    第一轮迭代完成

    ==================================================================================

    第二轮迭代,计算残差:
    在这里插入图片描述
    接下来我们计算这个残差r2和基b1基b3的点积。
    在这里插入图片描述
    在第二轮迭代中选择基b1

    贡献度值,即:
    在这里插入图片描述
    第二轮迭代完成

    ==================================================================================

    第三轮迭代,计算残差:
    在这里插入图片描述
    接下来我们计算这个残差r3和基b2基b3的点积。
    在这里插入图片描述
    在第三轮迭代中选择基b2

    贡献度值,即:
    在这里插入图片描述
    第三轮迭代完成

    ==================================================================================

    第四轮迭代,计算残差:
    在这里插入图片描述
    重复上述过程,直到残差达到指定阈值。

    ==================================================================================
    以上便是MP算法的计算流程。

    第二个例子

    现在我们继续看另一个例子,给定:
    在这里插入图片描述
    很容易就可以得到:

    在这里插入图片描述
    现在给定A和Y用OMP计算X。

    步骤如下:

    1.A中有四个基
    在这里插入图片描述
    2. 因为基的长度不是1,我们需要把它们都归一化。
    在这里插入图片描述
    3. 计算这些归一化基的贡献度:
    在这里插入图片描述
    4.第二个标准化的基b ̂2的贡献最高。因此,我们将第二个基b2放入矩阵Anew中。
    在这里插入图片描述
    5.计算最小二乘问题的解λ1
    在这里插入图片描述
    因为这个λ1是基b2的系数,那么:
    在这里插入图片描述
    6.计算残差
    在这里插入图片描述
    7.重复第3步。计算基b ̂1,b ̂3,b ̂4中哪一个基对残差r1的贡献最大。
    在这里插入图片描述
    在这里,在基b ̂3的贡献最大的,是1.7902。
    8.将基b3而不是(b ̂3)添加到Anew中。
    在这里插入图片描述
    4. 计算最小二乘问题的解λ2
    在这里插入图片描述
    由于λ2对应基b2和基b3,因此:
    在这里插入图片描述
    10.计算残差在这里插入图片描述
    11.重复第3步,直到最后一次迭代。
    12.计基b1和基b4的贡献度
    在这里插入图片描述
    因此,基b4的贡献度最大,为0.7308.
    5. 将基b4而不是(b ̂4)添加到Anew中。
    在这里插入图片描述
    6. 计算最小二乘问题的解λ3.在这里插入图片描述
    因为λ2对应基b2、基b3和基b4,因此,
    在这里插入图片描述
    15.计算残差
    在这里插入图片描述
    16.迭代到此停止,因为残差已经全部是0了。重构的原始向量xrec为:
    在这里插入图片描述

    MP和OMP的缺点

    OMP很快。与凸优化相比,速度非常快。这是贪婪算法的一个优点。但是OMP也遭受矩阵A高度相干性的困扰。矩阵中的相干性是什么?简单来说,矩阵A中的相干性表示矩阵A中两列之间的相似性(similarity)。
    请看下面两个矩阵:
    在这里插入图片描述
    矩阵A2具有非常高的相干性,因为第二类和第一列非常相似。矩阵的相干性较低,因为第二列与第一列和第三列不太相似。相干值用μ表示:
    在这里插入图片描述
    μ的值在0和1之间。如果μ的值很大,则通常OMP给出一个错误的重构结果。

    试试下面的例子:
    已知
    在这里插入图片描述
    求出Y。
    然后,给定A和Y,使用MP和OMP求X。
    下面直接给出结果:
    在这里插入图片描述
    可以看出,MP和OMP有类似特性,且都不能计算出正确结果

    后记

    正如我们在上面的过程中看到的那样,我们现在注意到以下事实。
    (1)我们从OMP计算得出的最高贡献来自标准化基础。并非来自最初的非标准化基础。
    (2)如果测量或重构矩阵A中的所有基础都已被归一化,则无需进行归一化步骤即可继续进行计算。
    (3)每个基的贡献是由残基和归一化基的点积完成的。
    (4)在MP的情况下,重建系数xrec是根据基和残基的点积计算得出的,而在OMP上,则是根据Anew到Y的最小二乘解计算得出的系数xrec。最小二乘解λ表示基对应位置的xrec。Anew本身是从非标准化的基得来的,这个过程需要时间。因此,OMP比MP更慢。
    (5)残差r是由Y、Anew和λ计算得来的
    (6)迭代将持续与A中的行数(即M)一样多。或者,如果信号k的稀疏度已知,则迭代重复k次。如果k <M,则知道k将有助于我们更快地进行计算。但是,如果k未知,则迭代次数最多应为M。

    参考文献

    [1]: Usman, Koredianto, (2017), Introduction to Orthogonal Matching Pursuit, Telkom University, Online : http://korediantousman.staff.telkomuniversity.ac.id.
    [2]: Usman, Koredianto, 2017, Introduction to Matching Pursuit, Telkom University, Online : http://korediantousman.staff.telkomuniversity.ac.id/tutorial.

    展开全文
  • 设计了一种基于FPGA的正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法的硬件优化结构,对OMP算法进行了改进,大大减少了乘法运算次数;在矩阵分解部分采用了交替柯列斯基分解(Alternative Cholesky ...
  • 稀疏编码中的正交匹配追踪(OMP)代码
                   


    分类: AI and Computer Vision908人阅读评论(11)收藏举报

    最近在看有关匹配追踪与相关优化的文章,发现了这篇http://blog.csdn.net/scucj/article/details/7467955,感觉作者写得很不错,这里也再写写自己的理解。文中有Matlab的代码,为了方便以后的使用,我顺便写了一个C++版本的,方便与OpenCV配合。


    为了方便理解,我将所有向量都表示为平面二维向量,待用原子表征的目标向量y,用红色表示,原子向量用蓝色表示,残差向量用绿色表示。于是匹配追踪算法(MP)实际上可以用下图表示。


    注意原子向量和目标向量都已归一化到单位长度,MP算法首先在所有原子向量中找到向OA投影最大的向量,即OB,然后计算OA - <OA, OB>OB,其中的<OA, OB>OB也就是图中的OD了,被OA减掉后,剩下的就是残差DA,根据初中几何知识就可以知道,DA是一定垂直于OB的,也就是说MP的残差始终与最近选出来的那个原子向量正交。


    而OMP要做的,就是让残差与已经选出来的所有原子向量都正交,这一点在图上不好画出来,但上面的那篇博文写的已经很详尽了,这里不再敖述。

    下面是用C++实现的OMP算法,具体流程参考上面博文中的一张图:



    其中的最小二乘可以直接通过矩阵运算得到,也可以使用OpenCV的solve方法,该方法专门用于求解线性方程组或最小二乘问题。代码如下:

    void OrthMatchPursuitvector<Mat>& dic,//字典 Mat& target,  float min_residual,  int sparsity,  Mat& x,  //返回每个原子对应的系数; vector<int>& patch_indices  //返回选出的原子序号 ){ Mat residual = target.clone();  Mat phi; //保存已选出的原子向量 x.create(0, 1, CV_32FC1);  float max_coefficient; unsigned int patch_index; for(;;) {  max_coefficient = 0;  for (int i = 0; i < dic.size(); i++)  {   float coefficient = (float)dic[i].dot(residual);    if (abs(coefficient) > abs(max_coefficient))   {   max_coefficient = coefficient;   patch_index = i;   }  }  patch_indices.push_back(patch_index); //添加选出的原子序号    Mat& matched_patch = dic[patch_index];  if (phi.cols == 0)   phi = matched_patch;  else   hconcat(phi,matched_patch,phi); //将新原子合并到原子集合中(都是列向量)    x.push_back(0.0f); //对系数矩阵新加一项  solve(phi, target, x, DECOMP_SVD); //求解最小二乘问题  residual = target - phi*x;  //更新残差  float res_norm = (float)norm(residual);  if (x.rows >= sparsity || res_norm <= min_residual) //如果残差小于阈值或达到要求的稀疏度,就返回   return; }}

    代码写得有点乱,基本上完全按照算法步骤来的,应该还有很大的性能提升空间。


    ===================================无耻的分割线========================================

    之前说上面的代码还有很大的优化空间,这几天捣鼓了一下,发现优化还是很有成效的,下面是具体方法。

    为方便起见,这里用A代表从字典当中选出的原子的集合,对于上面求解最小二乘的一步,可以表示为下式:


    其中的是一个对称正定矩阵,现在假如经过一次搜索后,又找到了一个原子向量v,那么新的原子集合可以表示为:


    那么用这个新的原子集合计算x时,可以得到:


    可以看到,新的集合乘积,有一部分是上次的结果(上式最右边矩阵的第一个元素),因此没有必要每次都从新计算,而只需对原来的矩阵更新一列和一行就行了。同时,上式的第二和第三个元素互为转置,也只需要计算其中一个,第四个元素是v的二范数的平方,直接调用norm()函数求得。

    在对矩阵进行行和列的添加时,我放弃了使用vconcat和hconcat方法,这两个方法效率较低,每次添加都会把原来的部分复制一遍。我现在一次分配好需要的大小,然后通过Mat的括号操作符取需要的子矩阵进行更新和计算。

    求解x时,还有一步是求的逆,既然是对称正定的,可以进行Cholesky分解,那么它的逆也可以很快求出。刚好OpenCV中有这样的方法,即调用inv()方法时,用DECOMP_CHOLESKY作为参数,根据官方文档,这样的速度是普通矩阵求逆的两倍!

    我做了一个测试,使用一个有600多个原子的字典,每个原子的维度为200,稀疏度设定为10,匹配一个信号,原来的方法需要200ms左右,而用上面的方法优化后,只需10ms!!快了一个数量级!!


    下面是优化后的代码:

    void DictionaryLearning::OrthMatchPursuit( Mat& target,  float min_residual,  int sparsity, //Store matched patches' coefficient vector<float>& coefficients,  //Store matched patches vector<DicPatch>& matched_patches, //Store indices of matched patches vector<int>& matched_indices ){ Mat residual = target.clone();  //the atoms' set; Mat ori_phi = Mat::zeros(m_vec_dims,sparsity,CV_32FC1); Mat phi; //phi.t()*phi which is a SPD matrix Mat ori_spd = Mat::ones(sparsity,sparsity,CV_32FC1); Mat spd = ori_spd(Rect(0,0,1,1));  //reserve enough memory. matched_patches.reserve(sparsity); matched_indices.reserve(sparsity); float max_coefficient; int matched_index; deque<DicPatch>::iterator matched_patch_it;  for(int spars = 1;;spars++) {  max_coefficient = 0;  matched_index = 0;  int current_index = 0;  for (deque<DicPatch>::iterator patch_it = m_patches.begin();    patch_it != m_patches.end();    ++patch_it   )  {   Mat& cur_vec = (*patch_it).vector;   float coefficient = (float)cur_vec.dot(residual);    //Find the maxmum coefficient   if (abs(coefficient) > abs(max_coefficient))   {    max_coefficient = coefficient;    matched_patch_it = patch_it;    matched_index = current_index;   }   current_index++;  }  matched_patches.push_back((*matched_patch_it));  matched_indices.push_back(matched_index);    Mat& matched_vec = (*matched_patch_it).vector;    //update the spd matrix via symply appending a single row and column to it.  if (spars > 1)  {   Mat v = matched_vec.t()*phi;   float c = (float)norm(matched_vec);   Mat new_row = ori_spd(Rect(0, spars - 1, spars - 1, 1));   Mat new_col = ori_spd(Rect(spars - 1, 0, 1, spars - 1));   v.copyTo(new_row);   ((Mat)v.t()).copyTo(new_col);   *ori_spd.ptr<float>(spars - 1, spars - 1) = c*c;   spd = ori_spd(Rect(0, 0, spars, spars));  }  //Add the new matched patch to the vectors' set.  phi = ori_phi(Rect(0, 0, spars, m_vec_dims));  matched_vec.copyTo(phi.col(spars - 1));     //A SPD matrix! Use Cholesky process to speed up.  Mat x = spd.inv(DECOMP_CHOLESKY)*phi.t()*target;  residual = target - phi*x;    float res_norm = (float)norm(residual);  if (spars >= sparsity || res_norm <= min_residual)  {   coefficients.clear();   coefficients.reserve(x.cols);   x.copyTo(coefficients);   return;  } }}
               

    再分享一下我老师大神的人工智能教程吧。零基础!通俗易懂!风趣幽默!还带黄段子!希望你也加入到我们人工智能的队伍中来!https://blog.csdn.net/jiangjunshow

    展开全文
  • 受压缩感知理论的启发, 基于参数向量所具有的稀疏特性, 提出一种新的阈值正交匹配追踪算法辨识系统的参数和时滞. 仿真实验表明, 所提出的算法能在少量采样时有效地辨识系统参数、估计未知时滞, 同时验证了算法的...
  • 转载自:...  文献[1]中给出了施密特(Schimidt)正交化的过程: 上面的的[x,y]表示向量内积,[x,y]=xTy=yTx=[x,y]。施密特正交化公式中的br实际上可写为: 分子之所以可以这么变化是由于[x,y]
  • 利用信号噪声比和均方误差对改进后的正交匹配追踪算法进行去噪和识别的检验,以此判断输电信号的识别效果。通过对改进后的算法进行实验仿真表明:在信号的去噪过程中,MPA算法最接近原始输电信号,而改进后的正交...
  • 正交匹配追踪(OMP)C++代码

    千次阅读 2016-12-01 19:19:21
    一、前言 本文参考自两篇博文... OMP(正交匹配追踪)理论早在90年代就提出来了,其为将信号分解为超完备字典上的稀疏表示的经典方法之一,这两篇博客分析得很透出,原理上这里不再重复。
  • 现有的正交匹配追踪(OMP)算法都是在给定迭代次数(待重建图像的稀疏度)的条件下重建,这使其需要通过非常多的线性测量来保证精确重建.为此,文中提出一种改进的后退型最优OMP方法:首先利用最优正交匹配追踪(OOMP)算法在...
  • 压缩感知重构算法之正交匹配追踪(OMP)

    万次阅读 多人点赞 2015-04-19 17:26:12
     前面经过几篇的基础铺垫,本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码,并且给出单次测试例程代码、测量数M重构成功概率关系曲线绘制例程代码、信号稀疏度K重构成功概率关系曲线绘制例程代码。
  • 针对含有未知时滞的多输入输出误差系统的时滞参数辨识问题,提出一种基于辅助模型的正交匹配追踪迭代算法.首先,由于各输入通道的时滞未知,通过设定输入回归长度,对系统模型进行过参数化,得到一个高维的辨识模型,且...
  • 1.1.9.正交匹配追踪OrthogonalMatchingPursuit和orthogonal_mp实现OMP算法,用于近似线性模型的拟合,其中非线性系数的数量(即L 0伪范数)受到约束。作为最小角度回归等前向特征选择方法,正交匹配追踪可以用固定...
  • 最近在看有关匹配追踪与相关优化的文章,发现了这篇http://blog.csdn.net/scucj/article/details/7467955,感觉作者写得很不错,这里也再写写自己的理解。文中有Matlab的代码,为了方便以后的使用,我顺便写了一个...
  •  前面经过几篇的基础铺垫,本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码,并且给出单次测试例程代码、测量数M重构成功概率关系曲线绘制例程代码、信号稀疏度K重构成功概率关系曲线绘制例程代码。 0、符号...
  • 压缩感知重构算法之正则化正交匹配追踪(ROMP)

    万次阅读 多人点赞 2015-04-25 10:32:56
     正交匹配追踪算法每次迭代均只选择残差最相关的一列,自然人们会想:“每次迭代是否可以多选几列呢?”,正则化正交匹配追踪(RegularizedOMP)就是其中一种改进方法。本篇将在上一篇《压缩感知重构算法之正交匹配...
  •  正交匹配追踪算法每次迭代均只选择残差最相关的一列,自然人们会想:“每次迭代是否可以多选几列呢?”,正则化正交匹配追踪(RegularizedOMP)就是其中一种改进方法。本篇将在上一篇《压缩感知重构算法之
  • 正交匹配追踪算法【Orthogonal Matching Pursuit (OMP)】0.简介OMP算法沿用了MP算法中的原子选择准则,即通过求余量r感知矩阵Φ中原子之间内积的绝对值来选择最佳原子组合,不同的是OMP在分解的每一步中对所选择的...
  • 稀疏表示与匹配追踪

    2017-07-23 22:45:46
     如果研究压缩感知重构算法,那么匹配追踪绝对是一个绕不过去的概念,虽然最原始的匹配追踪算法应用于重构的文献并不多,但基于匹配追踪改进后的正交匹配追踪算法在研究压缩感知的人群中可以说是家喻户晓、人尽皆知...
  • 广义正交匹配追踪(Generalized OMP, gOMP)算法可以看作为OMP算法的一种推广。OMP每次只选择残差相关最大的一个,而gOMP则是简单地选择最大的S个。之所以这里表述为"简单地选择"是相比于ROMP之类算法的,不进行任何...
  • 题目:正交匹配追踪(OMP)在稀疏分解压缩感知重构中的异同 &nbsp; &nbsp; &nbsp; &nbsp; 如果研究了稀疏分解再来研究压缩感知可能会有一个疑惑:在稀疏分解中有一个OMP算法,在压缩感知的重构算.....

空空如也

空空如也

1 2 3 4 5 ... 7
收藏数 133
精华内容 53
关键字:

匹配追踪与正交匹配追踪