块正交匹配追踪(BOMP)算法的实现 - 修复给定代码

信息处理 matlab 优化 压缩传感 线性代数 稀疏模型
2022-02-10 10:18:52

这是我的实现,它不起作用:

function [result,gamma,T]=block_OMP4(phi,signal,N,d,maxIter)
[m n]=size(phi);
r=signal;
gamma=[];
T=[];
i=1;
tol=1e-4;
err=1;
x=zeros(N*d,1);
val=[];
figure,
while  (i<=maxIter && err>tol)

    maxv=[];
    for j=1:N
        startx=(j-1)*d+1;
        endx=startx+d-1;
        block=phi(:,startx:endx);
        t=block'*r;
        maxv(j)=norm(t,2);
%        
    end
    [maxvl,maxvindx]=max(maxv);
    gamma=union(gamma,maxvindx);
    x2=lsqr(phi,r);
    xb=zeros(N*d,1);
    for i=1:length(gamma)
        stratx=(gamma(i)-1)*d;
        endx=stratx+d-1;
        xb(stratx:endx)=x2(stratx:endx);
        x(stratx:endx)=x(stratx:endx)+xb(stratx:endx);
    end

    y=phi*x;
    r=signal-y;

    val(i)=norm(r);
    err=norm(r)
    i=i+1;

end
result=x;
end

任何人都有一个有效的实施?

1个回答

答案取自实施块正交匹配追踪 (BOMP) 算法

块正交匹配追踪 (BOMP) 算法基本上是具有单一主要区别的正交匹配追踪 (OMP) 算法 - 我们没有选择最大化相关性的单一索引,而是选择了一组索引、矩阵的列子集和解决方案向量。

该算法的一个很好的参考在:

  1. 块正交匹配追踪算法的一个最优条件
  2. 块稀疏性:一致性和高效恢复

代码由以下给出:

function [ vX ] = SolveLsL0Bomp( mA, vB, numBlocks, paramK, tolVal )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsL0Omp( mA, vB, paramK, tolVal )
% Minimizes Least Squares of Linear System with L0 Constraint Using
% Block Orthogonal Matching Pursuit (OMP) Method.
% \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} subject to {\left\| x
% \right\|}_{2, 0} \leq K
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix (Fat Matrix). Assumed to be
%                           normlaized. Namely norm(mA(:, ii)) = 1 for any
%                           ii.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'. 
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numBlocks         -   Number of Blocks.
%                           The number of blocks in the problem structure.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - paramK            -   Parameter K.
%                           The L0 constraint parameter. Basically the
%                           maximal number of active blocks in the
%                           solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - tolVal            -   Tolerance Value.
%                           Tolerance value for equality of the Linear
%                           System.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  An Optimal Condition for the Block Orthogonal Matching Pursuit
%       Algorithm - https://ieeexplore.ieee.org/document/8404118.
%   2.  Block Sparsity: Coherence and Efficient Recovery - https://ieeexplore.ieee.org/document/4960226.
% Remarks:
%   1.  The algorithm assumes 'mA' is normalized (Each column).
%   2.  The number of columns in matrix 'mA' must be an integer
%       multiplication of the number of blocks.
%   3.  For 'numBlocks = numColumns' (Equivalent of 'numElmBlock = 1') the
%       algorithm becomes the classic OMP.
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre Process 'mA' by normalizing its columns.
% Release Notes:
%   -   1.0.000     19/08/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

numElmBlock = numCols / numBlocks;
if(round(numElmBlock) ~= numElmBlock)
    error('Number of Blocks Doesn''t Match Size of Arrays');
end

vActiveIdx      = false([numCols, 1]);
vR              = vB;
vX              = zeros([numCols, 1]);
activeBlckIdx   = [];

for ii = 1:paramK

    maxCorr         = 0;

    for jj = 1:numBlocks
        vBlockIdx = (((jj - 1) * numElmBlock) + 1):(jj * numElmBlock);

        currCorr = abs(mA(:, vBlockIdx).' * vR);
        if(currCorr > maxCorr)
            activeBlckIdx = jj;
            maxCorr = currCorr;
        end
    end

    vBlockIdx = (((activeBlckIdx - 1) * numElmBlock) + 1):(activeBlckIdx * numElmBlock);
    vActiveIdx(vBlockIdx) = true();

    vX(vActiveIdx) = mA(:, vActiveIdx) \ vB;
    vR = vB - (mA(:, vActiveIdx) * vX(vActiveIdx));

    resNorm = norm(vR);

    if(resNorm < tolVal)
        break;
    end

end


end