[转]广义正交匹配追踪(gOMP)

广义正交匹配追踪(Generalized OMP, gOMP)算法可以看作为OMP算法的一种推广，由文献[1]提出，第1作者本硕为哈工大毕业，发表此论文时在Korea University攻读博士学位。OMP每次只选择与残差相关最大的一个，而gOMP则是简单地选择最大的S个。之所以这里表述为“简单地选择”是相比于ROMP之类算法的，不进行任何其它处理，只是选择最大的S个而已。

下图为论文中提出的算法伪代码流程：

1 gOMP重构算法流程

参考文献的代码如下所示：

（1）构造K稀疏信号

% Generate K-sparse vector
%
% N         : original signal size.
% K         : sparsity level
%
%   Output parameters
% x_omp     : estimated signal
% iter_count: iteration count during estimating
%
% Written by Suhyuk (Seokbeop) Kwon
% Information System Lab., Korea Univ.
% http://isl.korea.ac.kr
function [x x_pos] = islsp_GenSparseVec(N, K)
    KPos    = K;
    if N/2 < K
        KPos    = N-K;
    end
    randPos     = ceil(N*rand( KPos, 1 ));
    randPos     = union(randPos,randPos);
    leftPOsLen  = KPos-length(randPos);
    while leftPOsLen > 0
        tmpPos  = ceil(N*rand( leftPOsLen, 1 ));
        randPos = union(tmpPos,randPos);
        leftPOsLen = KPos-length(randPos);
    end
    if KPos < K
        randPos = setxor((1:N), randPos);
    end
    x               = zeros( N, 1 );
    x(randPos)      = randn( K, 1 );
    x_pos           = randPos;
end

（2）gOMP函数

% Estimate the sparse signal x using generalized OMP
%
% y     : observation
% Phi   : sensing matrix
% K     : sparsity
% S     : selection length
%
%   Output parameters
% x_omp     : estimated signal
% iter_count: iteration count during estimating
%
% Written by Suhyuk (Seokbeop) Kwon
% Information System Lab., Korea Univ.
% http://isl.korea.ac.kr
function [x_ommp iter_count] = islsp_EstgOMP(y, Phi, K, S, zero_threshold)
% Check the parameters
    if nargin < 3
        error('islsp_EstgOMP : Input arguments y ,Phi and K must be specified.');
    end
        
    if nargin < 4
        S = max(K/4, 1);
    end
    
    if nargin < 5
        zero_threshold  = 1e-6;
    end
% Initialize the variables
    [nRows nCols]   = size(Phi);
    x_ommp          = zeros(size(Phi,2), 1);
    residual_prev   = y;
    supp            = [];
    iter_count      = 0;
    while (norm(residual_prev) > zero_threshold && iter_count < K)
        iter_count  = iter_count+1;
        [supp_mag supp_idx] = sort(abs(Phi'*residual_prev), 'descend');
        supp_n              = union(supp, supp_idx(1:S));
        if (length(supp_n) ~= length(supp)) && (length(supp_n) < nRows )
            x_hat           = Phi(:,supp_n)y;
            residual_prev   = y - Phi(:,supp_n)*x_hat;
            supp    = supp_n;
        else
            break;
        end
    end
    x_ommp(supp)    = Phi(:,supp)y;
    
    if nargout < 2
        clear('iter_count');
    end
end

（3）测试主函数

% Measurements size
m       = 50;
% Signal size
N       = 100;
% Sparsity level
K       = 20;
% Generate sensing matrix
Phi = randn(m,N)/sqrt(m);
% Generate sparse vector
[x x_pos]   = islsp_GenSparseVec(N, K);
y   = Phi*x;
% Using default parameters
[x1 itr1]   = islsp_EstgOMP(y, Phi, K);
% Find the sparse vector via selecting 4 indices
[x2 itr2]   = islsp_EstgOMP(y, Phi, K, 4);
% Find the sparse vector via selecting 4 indices until the residual becomes 1e-12
[x3 itr3]   = islsp_EstgOMP(y, Phi, K, 4, 1e-12);
disp('Mean square error');
[mse(x-x1) mse(x-x2) mse(x-x3)]
disp('Iteration number');
[itr1 itr2 itr3]

参考文献：

[1] Jian Wang, Seokbeop Kwon,Byonghyo Shim.  Generalized orthogonal matching pursuit, IEEE Transactions on Signal Processing, vol. 60, no. 12, pp.6202-6216, Dec. 2012.

Available at: http://islab.snu.ac.kr/paper/tsp_gOMP.pdf

[2] http://islab.snu.ac.kr/paper/gOMP.zip

[3] 彬彬有礼. 压缩感知重构算法之广义正交匹配追踪(gOMP)，http://blog.csdn.net/jbb0523/article/details/45693027