MP、OMP与施密特正交化(转载)

浅谈压缩感知（十九）：MP、OMP与施密特正交化

        </h1><div class="clear"></div><div class="postBody">

关于MP、OMP的相关算法与收敛证明，可以参考：http://www.cnblogs.com/AndyJee/p/5047174.html，这里仅简单陈述算法流程及二者的不同之处。

主要内容：

MP的算法流程及其MATLAB实现
OMP的算法流程以及MATLAB实现
MP与OMP的区别
施密特正交化与OMP的关系

一、MP（匹配追踪）的算法流程：

二、MP的MATLAB实现：

% MP：匹配追踪算法
% dictionary: 超完备字典
% x: 待表示信号
% M = 4; N = 10;
% Phi = randn(M,N); % 字典
% for nn = 1:N
%     Phi(:,nn) = Phi(:,nn)/norm(Phi(:,nn));
% end
% b = randn(M,1); % 信号
function x = MP(dictionary,x,iter)
[M,N] = size(dictionary);
residual = zeros(M,iter); %残差矩阵，保存每次迭代后的残差
residual(:,1) = x; %初始化残差为x
L = size(residual,2); %得到残差矩阵的列
pos_num = zeros(1,L); %用来保存每次选择的列序号
resi_norm = zeros(1,L); %用来保存每次迭代后的残差的2范数
resi_norm(1) = norm(x); %因为前面已初始化残差为x
iter_out = 1e-3;
iter_count = 0;

for mm = 1:iter
%迭代退出条件
if resi_norm(mm) < iter_out
break;
end
%求出dictionary每列与上次残差的内积
scalarproducts = dictionary’residual(:,mm);
%找到内积中最大的列及其内积值
[val,pos] = max(abs(scalarproducts));
%更新残差
residual(:,mm+1) = residual(:,mm) - scalarproducts(pos)dictionary(:,pos);
%计算残差的2范数（平方和再开根号）
resi_norm(mm+1) = norm(residual(:,mm+1));
%保存选择的列序号
pos_num(mm) = pos;
iter_count = iter_count + 1;
end
%绘出残差的2范数曲线
resi_norm = resi_norm(1:iter_count+1);
plot(resi_norm);grid;
%显示选择的字典原子
pos_num = pos_num(1:iter_count);
disp(pos_num);
%稀疏系数(稀疏表示)
dict = dictionary(:,pos_num);
y_vec = (dict’*dict)^(-1)dict’x;
disp(y_vec);
figure;plot(y_vec);

三、OMP（正交匹配追踪）的算法流程：

四、OMP的MATLAB实现：

% MP：匹配追踪算法
% dictionary: 超完备字典
% x: 待表示信号
% M = 4; N = 10;
% Phi = randn(M,N); % 字典
% for nn = 1:N
%     Phi(:,nn) = Phi(:,nn)/norm(Phi(:,nn));
% end
% b = randn(M,1); % 信号
function x = OMP(dictionary,x,iter)
[M,N] = size(dictionary);
residual = zeros(M,iter); %残差矩阵，保存每次迭代后的残差
residual(:,1) = x; %初始化残差为x
L = size(residual,2); %得到残差矩阵的列
pos_num = zeros(1,L); %用来保存每次选择的列序号
resi_norm = zeros(1,L); %用来保存每次迭代后的残差的2范数
resi_norm(1) = norm(x); %因为前面已初始化残差为x
iter_out = 1e-3;
iter_count = 0;
aug_mat = [];

for mm = 1:iter
%迭代退出条件
if resi_norm(mm) < iter_out
break;
end
%求出dictionary每列与上次残差的内积
scalarproducts = dictionary’residual(:,mm);
%找到内积中最大的列及其内积值
[val,pos] = max(abs(scalarproducts));
%最小二乘的增广矩阵
aug_mat = [aug_mat dictionary(:,pos)];
%最小二乘投影
proj_y = aug_mat(aug_mat’*aug_mat)^(-1)aug_mat’x;
%更新残差
residual(:,mm+1) = x - proj_y;
%计算残差的2范数（平方和再开根号）
resi_norm(mm+1) = norm(residual(:,mm+1));
%保存选择的列序号
pos_num(mm) = pos;
iter_count = iter_count + 1;
end
%绘出残差的2范数曲线
resi_norm = resi_norm(1:iter_count+1);
plot(resi_norm);grid;
%显示选择的字典原子
pos_num = pos_num(1:iter_count);
disp(pos_num);
%稀疏系数
dict = dictionary(:,pos_num);
y_vec = (dict’*dict)^(-1)dict’x;
disp(y_vec);
figure;plot(y_vec);

五、MP与OMP的区别：

OMP与MP的不同根本在于残差更新过程：OMP减去的Pe_m是e_m在所有被选择过的原子组成的矩阵Φ_t所张成空间上的正交投影，而MP减去的Pe_m是e_m在本次被选择的原子φ_m所张成空间上的正交投影。基于此，OMP可以保证已经选择过的原子不会再被选择。

六、施密特（Schimidt）正交化与OMP

1、施密特(Schimidt)正交化的过程:

上面的的[x,y]表示向量内积，[x,y]=x^Ty=y^Tx=[x,y]。施密特正交化公式中的b_r实际上可写为：

分子之所以可以这么变化是由于[x,y]实际上为一个数，因此[x,y]x=x[x,y]= xx^Ty。

2、OMP与施密特(Schimidt)正交化的关系：

结论：OMP分解过程，实际上是将所选原子依次进行Schimidt正交化，然后将待分解信号减去在正交化后的原子上各自的分量即可得残差。其实(式3)求残差的过程也是在进行施密特正交化。

3、验证OMP残差求解过程与Schmidt正交化的关系

% 验证OMP残差求解过程与Schmidt正交化的关系
%
clc;clear;close all;
M = 4; N = 10;
Phi = randn(M,N); % 字典
for nn = 1:NPhi(:,nn) = Phi(:,nn)/norm(Phi(:,nn));
end
b = randn(M,1); % 信号
res0 = b; % 初始化残差为待稀疏信号b
% OMP
% 选择字典第一个原子
c1 = Phi'* res0; % 求矩阵Phi各列与b的内积
[val1,pos1] = max(abs(c1)); % 找到内积中最大的列及其内积值
phit = [Phi(:,pos1)]; % 由所有选出的列组合的矩阵
Pphi = phit*(phit'*phit)^(-1)*phit'; % 正交投影变换矩阵
omp_res1 = res0 - Pphi*res0; % OMP用上一次残差减去残差在phit列空间的正交投影
omp_resb = b - Pphi*b; % OMP用待稀疏信号b减去b在phit列空间的正交投影
% Schimidt
x = Phi(:,pos1); % Schimidt正交化第一个向量
Px = x*(x'*x)^(-1)*x';
smt_res1 = res0 - Px*b; % 实际上是b - Px*b
% test
norm(omp_res1-omp_resb)
norm(omp_resb-smt_res1)

% OMP
% 选择字典第二列
c2 = Phi’ * omp_res1;
[val2,pos2] = max(abs(c2));
phit = [Phi(:,pos1) Phi(:,pos2)];
Pphi = phit*(phit’phit)^(-1)phit’;
omp_res2 = omp_res1 - Pphiomp_res1;
omp_resb = b - Pphib;
% Schimidt
y = Phi(:,pos2) - PxPhi(:,pos2); % Schimidt正交化第二个向量
Py = y(y’y)^(-1)y’;
smt_res2 = smt_res1 - Pyb; % 实际上是b - Pxb - Py*b，上一次残差减去b在第2列正交化所得z上的投影
% test
norm(omp_res2-omp_resb)
norm(omp_resb-smt_res2)

% OMP
% 选择字典第三列
c3 = Phi’ * omp_res2;
[val3,pos3] = max(abs(c3));
phit = [Phi(:,pos1) Phi(:,pos2) Phi(:,pos3)];
Pphi = phit*(phit’phit)^(-1)phit’;
omp_res3 = omp_res2 - Pphiomp_res2;
omp_resb = b - Pphib;
% Schimidt
z = Phi(:,pos3) - PxPhi(:,pos3) - PyPhi(:,pos3); % Schimidt正交化第三个向量
Pz = z*(z’z)^(-1)z’;
smt_res3 = smt_res2 - Pzb; % 实际上是b - Pxb - Pyb - Pzb，上一次残差减去b在第3列正交化所得z上的投影
% test
norm(omp_res3-omp_resb)
norm(omp_resb-smt_res3)

参考文章：

http://blog.csdn.net/jbb0523/article/details/45099655

http://blog.csdn.net/jbb0523/article/details/45100351

分类: 压缩感知compressive sensing, 数学Mathematics

<div id="blog_post_info">

好文要顶关注我收藏该文

AndyJee
关注 - 0
粉丝 - 324

+加关注

<div class="clear"></div>
<div id="post_next_prev"><a href="https://www.cnblogs.com/AndyJee/p/5091932.html" class="p_n_p_prefix">« </a> 上一篇：    <a href="https://www.cnblogs.com/AndyJee/p/5091932.html" title="发布于 2015-12-31 15:41">浅谈压缩感知（十八）：常见测量矩阵及其实现</a>
<br>
<a href="https://www.cnblogs.com/AndyJee/p/5112251.html" class="p_n_p_prefix">» </a> 下一篇：    <a href="https://www.cnblogs.com/AndyJee/p/5112251.html" title="发布于 2016-01-08 10:21">浅谈压缩感知（二十）：OMP与压缩感知</a>

posted @ 2016-01-07 17:44 AndyJee 阅读(3158) 评论(0) 编辑收藏

</div><!--end: topics 文章、评论容器-->

刷新评论刷新页面返回顶部

注册用户登录后才能发表评论，请登录或注册，访问网站首页。

【推荐】超50万C++/C#源码: 大型实时仿真组态图形源码
【推荐】零基础轻松玩转华为云产品，获壕礼加返百元大礼
【推荐】新手上天翼云，数十款云产品、新一代主机0元体验
【推荐】华为IoT平台开发者套餐9.9元起，购买即送免费课程

相关博文：
· 浅谈压缩感知（九）：正交匹配追踪算法OMP
· 浅谈压缩感知（二十二）：压缩感知重构算法之正则化正交匹配追踪（ROMP）
· 浅谈压缩感知（二十五）：压缩感知重构算法之分段正交匹配追踪（StOMP）
· 浅谈压缩感知（二十一）：压缩感知重构算法之正交匹配追踪（OMP）
· 浅谈压缩感知（二十）：OMP与压缩感知

    <div id="google_ads_iframe_/1090369/C2_0__container__" style="border: 0pt none;"><iframe id="google_ads_iframe_/1090369/C2_0" title="3rd party ad content" name="google_ads_iframe_/1090369/C2_0" width="468" height="60" scrolling="no" marginwidth="0" marginheight="0" frameborder="0" srcdoc="" style="border: 0px; vertical-align: bottom;" data-google-container-id="2" data-load-complete="true"></iframe></div></div>
</div>
<div id="under_post_kb">

最新 IT 新闻:
· Windows 10终于拿下操作系统市场半壁江山
· 藏在1.85亿人体内的隐形致癌病毒，有人确诊即是晚期
· 嘀嗒出行：用户突破1.3亿，车主突破1500万
· 无人机泡沫破裂：创企纷纷关门风投遭重创
· 基因编辑明星动物遇挫：无角牛体内发现细菌基因污染 | 专访
» 更多新闻...