【2】施密特(Schimidt)正交化与正交匹配追踪

转载自：http://blog.csdn.net/jbb0523/article/details/45100351

        文献[1]中给出了施密特(Schimidt)正交化的过程：

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

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

        上一篇《稀疏表示与匹配追踪》中详细的解释了匹配追踪（Matching Pursuit，MP）的流程，在最后给出了正交匹配追踪（Orthogonal Matching Pursuit，OMP）的流程，并指出OMP与MP的不同根本在于残差更新过程：OMP减去的Pe_m是在所有被选择过的原子组成的矩阵φ_t所张成空间上的正交投影，而MP减去的Pe_m是在本次被选择的原子φ_p所张成空间上的正交投影。正交投影及正交投影变换矩阵的概念可以参见《压缩感知中的数学知识：投影矩阵（projectionmatrix）》。

        首先给出一个结论：

        设OMP共从冗余字典中选择了r个原子，分别是a₁，a₂，……，a_r，根据正交匹配追踪的流程可以知道待分解信号x最后剩余的残差e_romp为

         (式1)

        该残差也可以表示为

         (式2)

其中矩阵A为选择的r个原子组成的矩阵，e_(r-1)omp为选择(r-1)个原子时的残差。

        将选择的r个原子a₁，a₂，……，a_r进行施密特正交化：

则残差e_romp还可以写为

         (式3)

(式1)一般出现在稀疏分解算法中，(式2)一般出现在重构算法中，(式3)是自己琢磨出来的（受到沙威的文档中提到的施密特正交化的启发，但沙威只限于向量情况下，详情可参见[3]，此处相当于一个推广）。

        而实际上(式1)、(式2)和(式3)三种正交匹配追踪算法的残差计算方法是等价的，纯数学证明比较复杂，这里给一个MATLAB程序，可以验证三者的等价性：

[plain]  view plain  copy

1. %% 验证OMP残差求解过程与Schmidt正交化的关系  
2. M = 4;N = 10;  
3. Phi = randn(M,N);  
4. for nn = 1:N  
5.     Phi(:,nn) = Phi(:,nn)/norm(Phi(:,nn));  
6. end  
7. b = randn(M,1);  
8. e0 = b;%初始化残差为待稀疏信号b  
9. %选第1列  
10. c1 = Phi'*e0;%求出矩阵Phi每列与b的内积  
11. [val1,pos1]=max(abs(c1));%找到内积中最大的列及其内积值  
12. phit = [Phi(:,pos1)];%由所有选出的列组合的矩阵  
13. Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵  
14. e1ompe = e0 - Pphi*e0;%OMP用上一次残差减去残差在phit列空间的正交投影  
15. e1ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影  
16. x = Phi(:,pos1);%Schimidt正交化第一个向量  
17. Px = x*(x'*x)^(-1)*x';  
18. %实际上是b - Px*b  
19. e1ompsmt = e0 - Px*b;  
20. e1 = e1ompe;  
21. norm(e1ompe-e1ompb)+norm(e1ompsmt-e1ompb)  
22. %选第2列  
23. c2 = Phi'*e1;%求出矩阵Phi每列与e1的内积  
24. [val2,pos2]=max(abs(c2));%找到内积中最大的列及其内积值  
25. phit = [Phi(:,pos1) Phi(:,pos2)];%由所有选出的列组合的矩阵  
26. Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵  
27. e2ompe = e1 - Pphi*e1;%OMP用上一次残差减去残差在phit列空间的正交投影  
28. e2ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影  
29. y = Phi(:,pos2) - Px*Phi(:,pos2);%Schimidt正交化第二个向量  
30. Py = y*(y'*y)^(-1)*y';  
31. %实际上是b - Px*b - Py*b  
32. e2ompsmt = e1 - Py*b;%上一次残差减去b在第2列正交化所得z上的投影  
33. e2 = e2ompe;  
34. norm(e2ompe-e2ompb)+norm(e2ompsmt-e2ompb)  
35. %选第3列  
36. c3 = Phi'*e2;%求出矩阵Phi每列与e2的内积  
37. [val3,pos3]=max(abs(c3));%找到内积中最大的列及其内积值  
38. phit = [Phi(:,pos1) Phi(:,pos2) Phi(:,pos3)];%由所有选出的列组合的矩阵  
39. Pphi = phit*(phit'*phit)^(-1)*phit';%正交投影变换矩阵  
40. e3ompe = e2 - Pphi*e2;%OMP用上一次残差减去残差在phit列空间的正交投影  
41. e3ompb = b - Pphi*b;%OMP用待稀疏信号b减去b在phit列空间的正交投影  
42. z = Phi(:,pos3) - Px*Phi(:,pos3) - Py*Phi(:,pos3);%Schimidt正交化第三个向量  
43. Pz = z*(z'*z)^(-1)*z';  
44. %实际上是b - Px*b - Py*b - Pz*b  
45. e3ompsmt = e2 - Pz*b;%上一次残差减去b在第3列正交化所得z上的投影  
46. e3 = e3ompe;  
47. norm(e3ompe-e3ompb)+norm(e3ompsmt-e3ompb)  

程序最后一行的输出结果为零，说明三个残差是相等的。

        这里可以看出OMP与Schimidt正交化的关系：

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

        有个关键问题还是要说的，分解后在所选择各原子上的系数是多少呢？答案其实也很简单，各个系数是(A^TA)^-1A^Tx，即最小二乘解，这个解是一个列向量，每一个元素分别是组成矩阵A的各原子的线性组合系数，这个在下一篇《正交匹配追踪(OMP)在稀疏分解与压缩感知重构中的异同》也会明确再次说明。

        同理，若设MP共从冗余字典中选择了r个原子，分别是a₁，a₂，……，a_r，根据匹配追踪的流程可以知道待分解信号x每次迭代后剩余的残差e_rmp为

比较式(3)的第2个等号表示的e_romp与此处的e_rmp也可以体会出OMP与MP的区别吧。

【参考文献】

【1】同济大学数学系. 线性代数（第五版）[M].高等教育出版社，2007：114.

【2】Joel A. Tropp and AnnaC. Gilbert . Signal Recovery From Random Measurements Via Orthogonal MatchingPursuit[J]. IEEE Transactions on Information Theory, VOL. 53, NO. 12, DECEMBER2007.

【3】沙威. “压缩传感”引论.http://www.eee.hku.hk/~wsha/Freecode/Files/Compressive_Sensing.pdf