正交匹配追踪算法（OMP）

---

前言

记录OMP算法的学习过程。

---

一、信号模型和逆问题

对于非齐次线性方程组$$Ax=b$$式中$b\in R^m,A\in R^{m\ast 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的位置。

---

三、伪代码

需要注意的是有个残差的更新过程，实际上原始残差就是$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')

实验结果

可见恢复了原始信号

总结

记录学习过程~