块正交匹配追踪(BOMP)算法的实现

信息处理 matlab 压缩传感 优化 线性代数 稀疏模型
2022-02-07 17:07:42

如何在 MATLAB 中实现块正交匹配追踪 (BOMP) 算法?

2个回答

块正交匹配追踪 (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


MATLAB 代码可在我的StackExchange Signal Processing Q60197 GitHub 存储库中找到(查看SignalProcessing\Q60197文件夹)。
在完整代码中,我将 Block 实现与 OMP 进行比较以验证实现。

我有使用图像莉娜

变换图像和恢复值的系数绘制为![恢复和原始系数的图] 2

恢复的图像 恢复图像

毫秒是 612896.781119