OMP 算法
问题描述
给定一个过完备字典矩阵,其中它的每列表示一种原型信号的原子。给定一个信号y,它可以被表示成这些原子的稀疏线性组合。信号 y 可以被表达为 y = Dx ,或者。 字典矩阵中所谓过完备性,指的是原子的个数远远大于信号y的长度(其长度很显然是n),即n<<k。
主要问题
数学模型
输入输出
输入:测量值y,字典D;
输出:稀疏稀疏向量x;
算法求解
代码实现
function [ theta ] = CS_OMP( y,A,t )
%CS_OMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-18
% Detailed explanation goes here
% y = Phi * x
% x = Psi * theta
% y = Phi*Psi * theta
% 令 A = Phi*Psi, 则y=A*theta
% 现在已知y和A,求theta
[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
[M,N] = size(A);%传感矩阵A为M*N矩阵
theta = zeros(N,1);%用来存储恢复的theta(列向量)
At = zeros(M,t);%用来迭代过程中存储A被选择的列
Pos_theta = zeros(1,t);%用来迭代过程中存储A被选择的列序号
r_n = y;%初始化残差(residual)为y
for ii=1:t%迭代t次,t为输入参数
product = A'*r_n;%传感矩阵A各列与残差的内积
[val,pos] = max(abs(product));%找到最大内积绝对值,即与残差最相关的列
At(:,ii) = A(:,pos);%存储这一列
Pos_theta(ii) = pos;%存储这一列的序号
A(:,pos) = zeros(M,1);%清零A的这一列,其实此行可以不要,因为它与残差正交
%y=At(:,1:ii)*theta,以下求theta的最小二乘解(Least Square)
theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解
%At(:,1:ii)*theta_ls是y在At(:,1:ii)列空间上的正交投影
r_n = y - At(:,1:ii)*theta_ls;%更新残差
end
theta(Pos_theta)=theta_ls;%恢复出的theta
end
总结
参考文献
MP算法和OMP算法及其思想:http://blog.csdn.net/scucj/article/details/7467955
Orthogonal Matching Pursuit(OMP)正交匹配追踪算法学习笔记:http://blog.sciencenet.cn/blog-810210-653094.html