正交匹配追踪算法（OMP）-Toy模板网

这篇具有很好参考价值的文章主要介绍了正交匹配追踪算法（OMP）。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

前言

记录OMP算法的学习过程。

一、信号模型和逆问题

对于非齐次线性方程组 $A x = b$ 式中 $\in R^m,A \in R^{m*n},x \in R^m$ 。
一般如果我们考虑 $A, x$ 已知，那么求 $b$ 是一个很简单的问题。
这个问题的逆问题为， $b, A$ 已知，去求 $x$ 。
当 $n > > m$ 时，该方程有无穷多解，如果我们想得到唯一解，就需要限定 $x$ ,实际上在压缩感知领域，就是限定 $x$ 是稀疏的，也就是 $x$ 中有很多0，在这种情况下去求解 $x$ ，其中 $k$ 为 $x$ 的稀疏度，实际上就是 $x$ 中非0的个数。

二、OMP原理

现在假设我们已知 $A, b$ ，想从中恢复出 $x$ ，显然我们需要充分利用 $x$ 是稀疏的这个事实，由线性代数知识显然可以知道 $b$ 其实是矩阵 $A$ 的列向量的线性组合，也就是 $x$ 作为权重与 $A$ 的列向量加权求和后得到的结果，由于 $x$ 是稀疏的，那么显然可以发现一个事实， $A$ 中仅仅有很少的列向量对 $b$ 做出了贡献，我们的目的就是找出这些对 $b$ 贡献较大的列向量，与此同时，根据列向量在 $A$ 中的位置，可以判断出 $x$ 中非零元素的位置，这就是OMP算法的基本思想。
那么关键是如何刻画 $A$ 中列向量对 $b$ 的贡献，在欧几里和空间中，我们常常用内积去定义两个向量的距离，实际上我们可以将 $b$ 往 $A$ 中列向量方向去投影，由此判断某个列向量对其的贡献。公式化表达 $contribution(b,a_i) =| \frac{<b,a_i>}{|a_i|}|$ 式中 $a_i$ 为 $A$ 的第 $i$ 个列向量，实际上在考虑投影时回出现正负号问题，我们只用考虑贡献的大小，而不去考虑方向，所以加上了绝对值。(<,>为内积运算)
实际上如果我们事先将 $A$ 矩阵的列向量单位化，那么公式可以简化为 $contribution(b,a_i) =|<b,a_i>|$

实际上给定稀疏度 $k$ ，我们只需要迭代 $k$ 次算法就可以求出 $A$ 中 $k$ 个贡献最大的列向量，以及 $x$ 中 $k$ 个不为0的位置。

三、伪代码

omp算法,压缩感知,matlab,矩阵,线性代数
需要注意的是有个残差的更新过程，实际上原始残差就是 $y$ ，每一次找到一个和他相关的列向量后，就得到了这部分的信息，所以要减去这部分信息，剩余的信息再去和列向量相关。

四、MATLAB代码

clear;
clc;
m = 64;
n = 256; % n>>m;
CN  = [];

A = randn(m,n);
x = zeros(n,1);
x(1) =0.4;
x(40) =0.6;
x(32) = 0.8;  % 构造稀疏向量x

k = 3;  % 代表稀疏度
b = A * x;
%% initialization 
r = b;   % 初始残差
Cn = []; % 用于记录存放的列的序号
An = []; % 用于存放列的列向量

%% Normalization
abs_colmn = sqrt(sum(A.^2)); % 每列的模长
abs_matrix = repmat(abs_colmn,m,1);
norm_A  = A./abs_matrix;
% disp(sum(norm_A.^2))
%% 迭代求解
for ii = 1:k
    product = norm_A' * r; % 实际上就是A的每个列向量与r相乘
    [val,index] = max(abs(product));
    Cn = [Cn index];
    An = [An norm_A(:,index)];
    xk   = inv(An' * An) * An' * b; % 最小二乘解
    r    = b - An * xk;   
end
x_recovery = zeros(n,1);
x_recovery(Cn) = xk;

figure;
subplot(2,1,1)
stem(x);
title('origin signal')
subplot(2,1,2)
stem(x_recovery ./ (abs_colmn')); % 反归一化
title('recovery signal')

实验结果
omp算法,压缩感知,matlab,矩阵,线性代数
可见恢复了原始信号