淺談壓縮感知（十九）：MP、OMP與施密特正交化

淺談壓縮感知（十九）：MP、OMP與施密特正交化

關(guān)于MP、OMP的相關(guān)算法與收斂證明，這里僅簡單陳述算法流程及二者的不同之處。

主要內(nèi)容：

MP的算法流程及其MATLAB實現(xiàn)OMP的算法流程以及MATLAB實現(xiàn)MP與OMP的區(qū)別施密特正交化與OMP的關(guān)系一、MP（匹配追蹤）的算法流程：

二、MP的MATLAB實現(xiàn)：

%?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范數(shù)
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每列與上次殘差的內(nèi)積
????scalarproducts?=?dictionary'*residual(:,mm);
?????%找到內(nèi)積中最大的列及其內(nèi)積值
????[val,pos]?=?max(abs(scalarproducts));
????%更新殘差
????residual(:,mm+1)?=?residual(:,mm)?-?scalarproducts(pos)*dictionary(:,pos);
????%計算殘差的2范數(shù)（平方和再開根號）
????resi_norm(mm+1)?=?norm(residual(:,mm+1));
?????%保存選擇的列序號
????pos_num(mm)?=?pos;
????iter_count?=?iter_count?+?1;
end
%繪出殘差的2范數(shù)曲線
resi_norm?=?resi_norm(1:iter_count+1);
plot(resi_norm);grid;
%顯示選擇的字典原子
pos_num?=?pos_num(1:iter_count);
disp(pos_num);
%稀疏系數(shù)(稀疏表示)
dict?=?dictionary(:,pos_num);
y_vec?=?(dict'*dict)^(-1)*dict'*x;
disp(y_vec);
figure;plot(y_vec);

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

四、OMP的MATLAB實現(xiàn)：

%?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范數(shù)
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每列與上次殘差的內(nèi)積
????scalarproducts?=?dictionary'*residual(:,mm);
????%找到內(nèi)積中最大的列及其內(nèi)積值
????[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范數(shù)（平方和再開根號）
????resi_norm(mm+1)?=?norm(residual(:,mm+1));
?????%保存選擇的列序號
????pos_num(mm)?=?pos;
????iter_count?=?iter_count?+?1;
end
%繪出殘差的2范數(shù)曲線
resi_norm?=?resi_norm(1:iter_count+1);
plot(resi_norm);grid;
%顯示選擇的字典原子
pos_num?=?pos_num(1:iter_count);
disp(pos_num);
%稀疏系數(shù)
dict?=?dictionary(:,pos_num);
y_vec?=?(dict'*dict)^(-1)*dict'*x;
disp(y_vec);
figure;plot(y_vec);

五、MP與OMP的區(qū)別：

OMP與MP的不同根本在于殘差更新過程：OMP減去的Pe_m是e_m在所有被選擇過的原子組成的矩陣Φ_t所張成空間上的正交投影，而MP減去的Pe_m是e_m在本次被選擇的原子φ_m所張成空間上的正交投影。基于此，OMP可以保證已經(jīng)選擇過的原子不會再被選擇。

六、施密特（Schimidt）正交化與OMP 1、施密特(Schimidt)正交化的過程:

上面的的[x,y]表示向量內(nèi)積，[x,y]=x^Ty=y^Tx=[x,y]。施密特正交化公式中的b_r實際上可寫為：

分子之所以可以這么變化是由于[x,y]實際上為一個數(shù)，因此[x,y]x=x[x,y]=?xx^Ty。

2、OMP與施密特(Schimidt)正交化的關(guān)系：

結(jié)論：OMP分解過程，實際上是將所選原子依次進行Schimidt正交化，然后將待分解信號減去在正交化后的原子上各自的分量即可得殘差。其實(式3)求殘差的過程也是在進行施密特正交化。

3、驗證OMP殘差求解過程與Schmidt正交化的關(guān)系

%?驗證OMP殘差求解過程與Schmidt正交化的關(guān)系
%
clc;clear;close?all;
M?=?4;?N?=?10;
Phi?=?randn(M,N);?%?字典
for?nn?=?1:N
????Phi(:,nn)?=?Phi(:,nn)/norm(Phi(:,nn));
end
b?=?randn(M,1);?%?信號
res0?=?b;?%?初始化殘差為待稀疏信號b
%?OMP
%?選擇字典第一個原子
c1?=?Phi'*?res0;?%?求矩陣Phi各列與b的內(nèi)積
[val1,pos1]?=?max(abs(c1));?%?找到內(nèi)積中最大的列及其內(nèi)積值
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?-?Pphi*omp_res1;
omp_resb?=?b?-?Pphi*b;
%?Schimidt
y?=?Phi(:,pos2)?-?Px*Phi(:,pos2);?%?Schimidt正交化第二個向量
Py?=?y*(y'*y)^(-1)*y';
smt_res2?=?smt_res1?-?Py*b;?%?實際上是b?-?Px*b?-?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?-?Pphi*omp_res2;?
omp_resb?=?b?-?Pphi*b;
%?Schimidt
z?=?Phi(:,pos3)?-?Px*Phi(:,pos3)?-?Py*Phi(:,pos3);??%?Schimidt正交化第三個向量
Pz?=?z*(z'*z)^(-1)*z';
smt_res3?=?smt_res2?-?Pz*b;?%?實際上是b?-?Px*b?-?Py*b?-?Pz*b，上一次殘差減去b在第3列正交化所得z上的投影
%?test
norm(omp_res3-omp_resb)
norm(omp_resb-smt_res3)