大家好，又见面了，我是你们的朋友全栈君。

转载自：https://blog.csdn.net/wyw921027/article/details/52102211

题目：压缩感知重构算法之迭代硬阈值(Iterative Hard Thresholding,IHT)

        本篇来介绍IHT重构算法。一般在压缩感知参考文献中，提到IHT时一般引用的都是文献【1】，但IHT实际上是在文献【2】中提出的。IHT并不是一种凸优化算法，它类似于OMP，是一种迭代算法，但它是由一个优化问题推导得到的。文献【1】和文献【2】的作者相同，署名单位为英国爱丁堡大学(University ofEdinburgh)，第一作者的个人主页见参考文献【3】，从个人主页来看，作者现在已到英国南安普敦大学(University of Southampton)，作者发表的论文均可以从其个人主页中下载。

        文献【1】的贡献是当把IHT应用于压缩感知重构问题时进行了一个理论分析：

1、迭代硬阈值(IHT)的提出

        值得一提的是，IHT在文献【2】中提出时并不叫Iterative Hard Thresholding，而是M-Sparse Algorithm，如下图所示：

        该算法是为了求解M-稀疏问题(M-sparse problem)式(3.1)而提出的，经过一番推导得到了迭代公式式(3.2)，其中HM(·)的含义参见式(3.3)。

        这里面最关键的问题是：式(3.2)这个迭代公式是如何推导得到的呢？

2、Step1:替代目标函数

        首先，将式(3.1)的目标函数用替代目标函数(surrogate objective fucntion)式(3.5)替换：

这里中的M应该指的是M-sparse，S应该指的是Surrogate。这里要求：

        为什么式目标函数式(3.1)可以用式(3.5) 替代呢？这得往回看一下了……

        实际上，文献【2】分别针对两个优化问题进行了讨论，本篇主要是文献中的第二个优化问题，由于两个问题有一定的相似性，所以文中在推导第二个问题时进行了一些简化，下面简单回顾一些必要的有关第一个问题的内容，第一个优化问题是：

将目标函数定义为：

        为了推导迭代公式（详见式(2.2)和式(2.3)）式(1.5)用如下替代目标函数代替：

        这里注意波浪下划线中提到的“[29]”(参见文献【4】)，surrogate objective function的思想来自这篇文件。然后注意对Φ的约束（第一个红框），之后以会有这个约束，个人认为是为了使式(2.5)后半部分大于等于零，即为了使

大于等于零（当y=z时这部分等于零）。由此自然就有了式(2.5)与式(1.5)两个目标函数的关系（第二个红框），这也很容易理解，将y=z代入式(2.5)自然可得这个关系。

        到此应该明白式(2.5)为什么可以替代式(1.5)了吧……

        而我们用式(3.5)替代目标函数

的道理是一模一样的。

        补充一点：有关对||Φ||2<1的约束文献【2】中有一处提到了如下描述：

3、Step2:替代目标函数变形

        接下来，式(3.5)进行了变形：

        这个式子是怎么来的呢？我们对式(3.5)进行一下推导：

        这里，后面三项2范数的平方是与y无关的项，因此可视为常量，若对参数y求最优化时这三项并不影响优化结果，可略去，因此就有了变形的结果，符号“∝”表示成正比例。

4、Step3:极值点的获得

        接下来文献【2】直接给出了极值点：

        注意文中提到了“landweder”，搜索一下可知经常出现的是“landweder迭代”，这个暂且不提。那么极值点是如何推导得到的呢？其实就是一个简单的配方，中学生就会的：

        令，则

        当，取得最小值

5、Step4:迭代公式的获得

        极值点得到了，替代目标函数的极小值也得到了：

        那么如何得到迭代公式式(3.2)呢？这时要注意，推导过程中有一个约束条件一直没管，即式(3.1)中的约束条件：

也就是向量y的稀疏度不大于M。综合起来说，替代函数的最小值是

那么怎么使这个最小值在向量y的稀疏度不大于M的约束下最小呢，显然是保留最大的M项（因为是平方，所以要取绝对值absolute value），剩余的置零（注意这里有个负号，所以要保留最大的M项）。

        至此，我们就得到了迭代公式式(3.2)。

6、IHT算法的MATLAB代码

         这里一共给出三个版本的IHT实现：

第一个版本：

        在作者的主页有官方版IHT算法MATLAB代码，但有些复杂，这里给出一个简化版的IHT代码，方便理解：

[plain] view plain copy

1. function [ y ] = IHT_Basic( x,Phi,M,mu,epsilon,loopmax )    
2. %IHT_Basic Summary of this function goes here    
3. %Version: 1.0 written by jbb0523 @2016-07-30    
4. %Reference:Blumensath T, Davies M E. Iterative Thresholding for Sparse Approximations[J].   
5. %Journal of Fourier Analysis & Applications, 2008, 14(5):629-654.   
6. %(Available at: http://link.springer.com/article/10.1007%2Fs00041-008-9035-z)  
7. %   Detailed explanation goes here    
8.     if nargin < 6    
9.         loopmax = 3000;    
10.     end    
11.     if nargin < 5      
12.         epsilon = 1e-3;      
13.     end     
14.     if nargin < 4      
15.         mu = 1;      
16.     end     
17.     [x_rows,x_columns] = size(x);      
18.     if x_rows<x_columns      
19.         x = x’;%x should be a column vector      
20.     end    
21.     n = size(Phi,2);    
22.     y = zeros(n,1);%Initialize y=0    
23.     loop = 0;    
24.     while(norm(x-Phi*y)>epsilon && loop < loopmax)    
25.         y = y + Phi’*(x-Phi*y)*mu;%update y    
26.         %the following two lines of code realize functionality of H_M(.)    
27.         %1st: permute absolute value of y in descending order    
28.         [ysorted inds] = sort(abs(y), ‘descend’);    
29.         %2nd: set all but M largest coordinates to zeros    
30.         y(inds(M+1:n)) = 0;    
31.         loop = loop + 1;    
32.     end    
33. end   

 第二个版本：（作者给出的官方版本）

        文件：hard_l0_Mterm.m(\sparsify_0_5\HardLab)

        链接：http://www.personal.soton.ac.uk/tb1m08/sparsify/sparsify_0_5.zip

[plain] view plain copy

1. function [s, err_mse, iter_time]=hard_l0_Mterm(x,A,m,M,varargin)  
2. % hard_l0_Mterm: Hard thresholding algorithm that keeps exactly M elements   
3. % in each iteration.   
4. %  
5. % This algorithm has certain performance guarantees as described in [1],  
6. % [2] and [3].  
7. %  
8. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
9. % Usage  
10. %  
11. %   [s, err_mse, iter_time]=hard_l0_Mterm(x,P,m,M,’option_name’,’option_value’)  
12. %  
13. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
14. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
15. %  
16. % Input  
17. %  
18. %   Mandatory:  
19. %               x   Observation vector to be decomposed  
20. %               P   Either:  
21. %                       1) An nxm matrix (n must be dimension of x)  
22. %                       2) A function handle (type “help function_format”   
23. %                          for more information)  
24. %                          Also requires specification of P_trans option.  
25. %                       3) An object handle (type “help object_format” for   
26. %                          more information)  
27. %               m   length of s   
28. %               M   non-zero elements to keep in each iteration  
29. %  
30. %   Possible additional options:  
31. %   (specify as many as you want using ‘option_name’,’option_value’ pairs)  
32. %   See below for explanation of options:  
33. %__________________________________________________________________________  
34. %   option_name    |     available option_values                | default  
35. %————————————————————————–  
36. %   stopTol        | number (see below)                         | 1e-16  
37. %   P_trans        | function_handle (see below)                |   
38. %   maxIter        | positive integer (see below)               | n^2  
39. %   verbose        | true, false                                | false  
40. %   start_val      | vector of length m                         | zeros  
41. %   step_size      | number                                     | 0 (auto)  
42. %  
43. %   stopping criteria used : (OldRMS-NewRMS)/RMS(x) < stopTol  
44. %  
45. %   stopTol: Value for stopping criterion.  
46. %  
47. %   P_trans: If P is a function handle, then P_trans has to be specified and   
48. %            must be a function handle.   
49. %  
50. %   maxIter: Maximum number of allowed iterations.  
51. %  
52. %   verbose: Logical value to allow algorithm progress to be displayed.  
53. %  
54. %   start_val: Allows algorithms to start from partial solution.  
55. %  
56. %  
57. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
58. %  
59. % Outputs  
60. %  
61. %    s              Solution vector   
62. %    err_mse        Vector containing mse of approximation error for each   
63. %                   iteration  
64. %    iter_time      Vector containing computation times for each iteration  
65. %  
66. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
67. %  
68. % Description  
69. %  
70. %   Implements the M-sparse algorithm described in [1], [2] and [3].  
71. %   This algorithm takes a gradient step and then thresholds to only retain  
72. %   M non-zero elements. It allows the step-size to be calculated  
73. %   automatically as described in [3] and is therefore now independent from   
74. %   a rescaling of P.  
75. %     
76. %     
77. % References  
78. %   [1]  T. Blumensath and M.E. Davies, “Iterative Thresholding for Sparse   
79. %        Approximations”, submitted, 2007  
80. %   [2]  T. Blumensath and M. Davies; “Iterative Hard Thresholding for   
81. %        Compressed Sensing” to appear Applied and Computational Harmonic   
82. %        Analysis   
83. %   [3] T. Blumensath and M. Davies; “A modified Iterative Hard   
84. %        Thresholding algorithm with guaranteed performance and stability”   
85. %        in preparation (title may change)   
86. % See Also  
87. %   hard_l0_reg  
88. %  
89. % Copyright (c) 2007 Thomas Blumensath  
90. %  
91. % The University of Edinburgh  
92. % Email: thomas.blumensath@ed.ac.uk  
93. % Comments and bug reports welcome  
94. %  
95. % This file is part of sparsity Version 0.4  
96. % Created: April 2007  
97. % Modified January 2009  
98. %  
99. % Part of this toolbox was developed with the support of EPSRC Grant  
100. % D000246/1  
101. %  
102. % Please read COPYRIGHT.m for terms and conditions.  
103.   
104.   
105. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
106. %                    Default values and initialisation  
107. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
108.   
109.   
110.   
111. [n1 n2]=size(x);  
112. if n2 == 1  
113.     n=n1;  
114. elseif n1 == 1  
115.     x=x’;  
116.     n=n2;  
117. else  
118.    error(‘x must be a vector.’);  
119. end  
120.       
121. sigsize     = x’*x/n;  
122. oldERR      = sigsize;  
123. err_mse     = [];  
124. iter_time   = [];  
125. STOPTOL     = 1e-16;  
126. MAXITER     = n^2;  
127. verbose     = false;  
128. initial_given=0;  
129. s_initial   = zeros(m,1);  
130. MU          = 0;  
131.   
132. if verbose  
133.    display(‘Initialising…’)   
134. end  
135.   
136. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
137. %                           Output variables  
138. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
139.   
140. switch nargout   
141.     case 3  
142.         comp_err=true;  
143.         comp_time=true;  
144.     case 2   
145.         comp_err=true;  
146.         comp_time=false;  
147.     case 1  
148.         comp_err=false;  
149.         comp_time=false;  
150.     case 0  
151.         error(‘Please assign output variable.’)          
152.     otherwise  
153.         error(‘Too many output arguments specified’)  
154. end  
155.   
156. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
157. %                       Look through options  
158. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
159.   
160. % Put option into nice format  
161. Options={};  
162. OS=nargin-4;  
163. c=1;  
164. for i=1:OS  
165.     if isa(varargin{i},’cell’)  
166.         CellSize=length(varargin{i});  
167.         ThisCell=varargin{i};  
168.         for j=1:CellSize  
169.             Options{c}=ThisCell{j};  
170.             c=c+1;  
171.         end  
172.     else  
173.         Options{c}=varargin{i};  
174.         c=c+1;  
175.     end  
176. end  
177. OS=length(Options);  
178. if rem(OS,2)  
179.    error(‘Something is wrong with argument name and argument value pairs.’)   
180. end  
181. for i=1:2:OS  
182.    switch Options{i}  
183.         case {‘stopTol’}  
184.             if isa(Options{i+1},’numeric’) ; STOPTOL     = Options{i+1};     
185.             else error(‘stopTol must be number. Exiting.’); end  
186.         case {‘P_trans’}   
187.             if isa(Options{i+1},’function_handle’); Pt = Options{i+1};     
188.             else error(‘P_trans must be function _handle. Exiting.’); end  
189.         case {‘maxIter’}  
190.             if isa(Options{i+1},’numeric’); MAXITER     = Options{i+1};               
191.             else error(‘maxIter must be a number. Exiting.’); end  
192.         case {‘verbose’}  
193.             if isa(Options{i+1},’logical’); verbose     = Options{i+1};     
194.             else error(‘verbose must be a logical. Exiting.’); end   
195.         case {‘start_val’}  
196.             if isa(Options{i+1},’numeric’) && length(Options{i+1}) == m ;  
197.                 s_initial     = Options{i+1};    
198.                 initial_given=1;  
199.             else error(‘start_val must be a vector of length m. Exiting.’); end  
200.         case {‘step_size’}  
201.             if isa(Options{i+1},’numeric’) && (Options{i+1}) > 0 ;  
202.                 MU     = Options{i+1};     
203.             else error(‘Stepsize must be between a positive number. Exiting.’); end  
204.         otherwise  
205.             error(‘Unrecognised option. Exiting.’)   
206.    end  
207. end  
208.   
209. if nargout >=2  
210.     err_mse = zeros(MAXITER,1);  
211. end  
212. if nargout ==3  
213.     iter_time = zeros(MAXITER,1);  
214. end  
215.   
216.   
217. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
218. %                        Make P and Pt functions  
219. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
220.   
221. if          isa(A,’float’)      P =@(z) A*z;  Pt =@(z) A’*z;  
222. elseif      isobject(A)         P =@(z) A*z;  Pt =@(z) A’*z;  
223. elseif      isa(A,’function_handle’)   
224.     try  
225.         if          isa(Pt,’function_handle’); P=A;  
226.         else        error(‘If P is a function handle, Pt also needs to be a function handle. Exiting.’); end  
227.     catch error(‘If P is a function handle, Pt needs to be specified. Exiting.’); end  
228. else        error(‘P is of unsupported type. Use matrix, function_handle or object. Exiting.’); end  
229.   
230. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
231. %                        Do we start from zero or not?  
232. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
233.   
234.   
235.   
236. if initial_given ==1;  
237.       
238.     if length(find(s_initial)) > M  
239.         display(‘Initial vector has more than M non-zero elements. Keeping only M largest.’)  
240.       
241.     end  
242.     s                   =   s_initial;  
243.     [ssort sortind]     =   sort(abs(s),’descend’);  
244.     s(sortind(M+1:end)) =   0;  
245.     Ps                  =   P(s);  
246.     Residual            =   x-Ps;  
247.     oldERR      = Residual’*Residual/n;  
248. else  
249.     s_initial   = zeros(m,1);  
250.     Residual    = x;  
251.     s           = s_initial;  
252.     Ps          = zeros(n,1);  
253.     oldERR      = sigsize;  
254. end  
255.   
256.   
257. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
258. %                 Random Check to see if dictionary norm is below 1   
259. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
260.   
261.           
262.         x_test=randn(m,1);  
263.         x_test=x_test/norm(x_test);  
264.         nP=norm(P(x_test));  
265.         if abs(MU*nP)>1;  
266.             display(‘WARNING! Algorithm likely to become unstable.’)  
267.             display(‘Use smaller step-size or || P ||_2 < 1.’)  
268.         end  
269.   
270.   
271. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
272. %                        Main algorithm  
273. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
274. if verbose  
275.    display(‘Main iterations…’)   
276. end  
277. tic  
278. t=0;  
279. done = 0;  
280. iter=1;  
281.   
282. while ~done  
283.       
284.     if MU == 0  
285.   
286.         %Calculate optimal step size and do line search  
287.         olds                =   s;  
288.         oldPs               =   Ps;  
289.         IND                 =   s~=0;  
290.         d                   =   Pt(Residual);  
291.         % If the current vector is zero, we take the largest elements in d  
292.         if sum(IND)==0  
293.             [dsort sortdind]    =   sort(abs(d),’descend’);  
294.             IND(sortdind(1:M))  =   1;      
295.          end    
296.   
297.         id                  =   (IND.*d);  
298.         Pd                  =   P(id);  
299.         mu                  =   id’*id/(Pd’*Pd);  
300.         s                   =   olds + mu * d;  
301.         [ssort sortind]     =   sort(abs(s),’descend’);  
302.         s(sortind(M+1:end)) =   0;  
303.         Ps                  =   P(s);  
304.           
305.         % Calculate step-size requirement   
306.         omega               =   (norm(s-olds)/norm(Ps-oldPs))^2;  
307.   
308.         % As long as the support changes and mu > omega, we decrease mu  
309.         while mu > (0.99)*omega && sum(xor(IND,s~=0))~=0 && sum(IND)~=0  
310. %             display([‘decreasing mu’])  
311.                       
312.                     % We use a simple line search, halving mu in each step  
313.                     mu                  =   mu/2;  
314.                     s                   =   olds + mu * d;  
315.                     [ssort sortind]     =   sort(abs(s),’descend’);  
316.                     s(sortind(M+1:end)) =   0;  
317.                     Ps                  =   P(s);  
318.                     % Calculate step-size requirement   
319.                     omega               =   (norm(s-olds)/norm(Ps-oldPs))^2;  
320.         end  
321.           
322.     else  
323.         % Use fixed step size  
324.         s                   =   s + MU * Pt(Residual);  
325.         [ssort sortind]     =   sort(abs(s),’descend’);  
326.         s(sortind(M+1:end)) =   0;  
327.         Ps                  =   P(s);  
328.           
329.     end  
330.         Residual            =   x-Ps;  
331.   
332.           
333.      ERR=Residual’*Residual/n;  
334.      if comp_err  
335.          err_mse(iter)=ERR;  
336.      end  
337.        
338.      if comp_time  
339.          iter_time(iter)=toc;  
340.      end  
341.   
342. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
343. %                        Are we done yet?  
344. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
345.        
346.          if comp_err && iter >=2  
347.              if ((err_mse(iter-1)-err_mse(iter))/sigsize<STOPTOL);  
348.                  if verbose  
349.                     display([‘Stopping. Approximation error changed less than ‘ num2str(STOPTOL)])  
350.                  end  
351.                 done = 1;   
352.              elseif verbose && toc-t>10  
353.                 display(sprintf(‘Iteration %i. — %i mse change’,iter ,(err_mse(iter-1)-err_mse(iter))/sigsize))   
354.                 t=toc;  
355.              end  
356.          else  
357.              if ((oldERR – ERR)/sigsize < STOPTOL) && iter >=2;  
358.                  if verbose  
359.                     display([‘Stopping. Approximation error changed less than ‘ num2str(STOPTOL)])  
360.                  end  
361.                 done = 1;   
362.              elseif verbose && toc-t>10  
363.                 display(sprintf(‘Iteration %i. — %i mse change’,iter ,(oldERR – ERR)/sigsize))   
364.                 t=toc;  
365.              end  
366.          end  
367.            
368.     % Also stop if residual gets too small or maxIter reached  
369.      if comp_err  
370.          if err_mse(iter)<1e-16  
371.              display(‘Stopping. Exact signal representation found!’)  
372.              done=1;  
373.          end  
374.      elseif iter>1   
375.          if ERR<1e-16  
376.              display(‘Stopping. Exact signal representation found!’)  
377.              done=1;  
378.          end  
379.      end  
380.   
381.      if iter >= MAXITER  
382.          display(‘Stopping. Maximum number of iterations reached!’)  
383.          done = 1;   
384.      end  
385.    
386. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
387. %                    If not done, take another round  
388. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
389.      
390.      if ~done  
391.         iter=iter+1;   
392.         oldERR=ERR;          
393.      end  
394. end  
395.   
396.   
397. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
398. %                  Only return as many elements as iterations  
399. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
400.   
401. if nargout >=2  
402.     err_mse = err_mse(1:iter);  
403. end  
404. if nargout ==3  
405.     iter_time = iter_time(1:iter);  
406. end  
407. if verbose  
408.    display(‘Done’)   
409. end  

第三个版本：

        文件：Demo_CS_IHT.m（部分）

        链接：http://www.pudn.com/downloads518/sourcecode/math/detail2151378.html

[plain] view plain copy

1. function hat_x=cs_iht(y,T_Mat,m)  
2. % y=T_Mat*x, T_Mat is n-by-m  
3. % y – measurements  
4. % T_Mat – combination of random matrix and sparse representation basis  
5. % m – size of the original signal  
6. % the sparsity is length(y)/4  
7.   
8. hat_x_tp=zeros(m,1);         % initialization with the size of original   
9. s=floor(length(y)/4);        % sparsity  
10. u=0.5;                       % impact factor  
11.   
12. % T_Mat=T_Mat/sqrt(sum(sum(T_Mat.^2))); % normalizae the whole matrix  
13.   
14. for times=1:s  
15.       
16.     x_increase=T_Mat’*(y-T_Mat*hat_x_tp);  
17.       
18.     hat_x=hat_x_tp+u*x_increase;  
19.       
20.     [val,pos]=sort((hat_x),’descend’);  % why? worse performance with abs()  
21.       
22.     hat_x(pos(s+1:end))=0;   % thresholding, keeping the larges s elements  
23.   
24.     hat_x_tp=hat_x;          % update  
25.   
26. end  

7、单次重构代码

 %压缩感知重构算法测试      

[plain] view plain copy

1. clear all;close all;clc;        
2. M = 64;%观测值个数        
3. N = 256;%信号x的长度        
4. K = 10;%信号x的稀疏度        
5. Index_K = randperm(N);        
6. x = zeros(N,1);        
7. x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的，且位置是随机的        
8. Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta        
9. Phi = randn(M,N);%测量矩阵为高斯矩阵    
10. Phi = orth(Phi’)’;      
11. A = Phi * Psi;%传感矩阵      
12. % sigma = 0.005;      
13. % e = sigma*randn(M,1);    
14. % y = Phi * x + e;%得到观测向量y        
15. y = Phi * x;%得到观测向量y      
16. %% 恢复重构信号x        
17. tic        
18. theta = IHT_Basic(y,A,K);   
19. % theta = cs_iht(y,A,size(A,2));  
20. % theta = hard_l0_Mterm(y,A,size(A,2),round(1.5*K),’verbose’,true);  
21. x_r = Psi * theta;% x=Psi * theta        
22. toc        
23. %% 绘图        
24. figure;        
25. plot(x_r,’k.-‘);%绘出x的恢复信号        
26. hold on;        
27. plot(x,’r’);%绘出原信号x        
28. hold off;        
29. legend(‘Recovery’,’Original’)        
30. fprintf(‘\n恢复残差：’);        
31. norm(x_r-x)%恢复残差       

        这里就不给出重构结果了，给出仿真结论：本人编的IHT基本版能够正常工作，但偶尔会重构失败；第二个版本hard_l0_Mterm.m重构效果很好；第三个版本Demo_CS_IHT.m重构效果很差，估计是作者疑问（why? worse performance with abs()），没有加abs取绝对值的原因吧……

8、结束语

8.1有关算法的名字

        值得注意的是，在文献【2】中将式(2.2)称为iterative hard-thresholding algorithm，而将式(3.2)称为M-sparse algorithm，在文献【1】中又将式(3.2)称为Iterative Hard Thresholding algorithm (IHTs)，一般简称IHT的较多，多余的s指的是s稀疏。可见算法的名称是也是一不断完善的过程啊……

8.2与GraDeS算法的关系

        如果你学习过GraDeS算法（参见http://blog.csdn.net/jbb0523/article/details/52059296），然后再学习本算法，是不是有一种似曾相似的感觉？

        没错，这两个算法的迭代公式几乎是一样的，尤其是文献【1】中的式(12)（如上图第二个红框）进一步拓展了该算法的定义。这个就跟CoSaMP与SP两个算法一样，在GraDeS的提出文献【5】中开始部分还提到了IHT，但后面就没提了，不知道作者是怎么看待这个问题的。如果非说二者有区别，那就是GraDeS的参数γ=1+δ2s，且δ2s<1/3。

        所以，有想法得赶紧写成论文发表出来，否则被抢了先机那就……

8.3重构效果问题

        另外，在GraDeS算法中提到该算法的重构效果不好，这里注意文献【2】中的一段话：

        也就是说，IHT作者也意识到了该种算法的问题，并提出了两种应用策略(two strategies for asuccessful application of the methods)。

8.4Landweber迭代

        在网上搜索“Landweber迭代”时找到了一段程序[6]：

[plain] view plain copy

1. function [x,k]=Landweber(A,b,x0)  
2. alfa=1/sum(diag(A*A’));  
3. k=1;  
4. L=200;  
5. x=x0;  
6. while k<L  
7.     x1=x;  
8.     x=x+alfa*A’*(b-A*x);  
9.     if norm(b-A*x)/norm(b)<0.005  
10.         break;  
11.     elseif norm(x1-x)/norm(x)<0.001  
12.         break;  
13.     end  
14.     k=k+1;  
15. end  

注意该程序的迭代部分“x=x+alfa*A’*(b-A*x);”，除了多了一些alfa系数外，这跟IHT不是基本一样么？或者说与GraDeS有什么区别？

        有关LandWeber迭代可参见文献：“Landweber L. An iteration formula for Fredholm integral equations of the first kind[J]. American journal of mathematics, 1951, 73(3): 615-624.”，此处不再多述。

8.5改进算法

        作者后来又提出了两个关于IHT的改进算法，分别是RIHT(Normalized IHT)[7]和AIHT(Accelerated IHT)[8]。

        提出RIHT主要是由于IHT有一些缺点[7]：

        新算法RIHT将会有如下优点：

        之所以作者提供的软件包（第二个版本IHT）重构效果更好是由于最新版的hard_l0_Mterm.m (\sparsify_0_5\HardLab)程序中已经更新成了RIHT。

        RIHT的算法流程如下：

        将IHT改进为AIHT后会有如下优点[8]：

        值得注意的是，AIHT应该是一类算法的总称（虽然作者只阐述了两种实现策略），这个类似于FFT是所有DFT快速算法的总称：

8.6稀疏度对IHT的影响

        自己可以试一下，IHT输入参数中的稀疏度并不是很关键，若实际稀疏度为K，则稀疏度这个输入参数只要不小于K就可以了，重构效果都挺不错的，比如第三个版本的IHT程序，作者直接将稀疏度定义为信号y长度的四分之一。

8.7作者去向？

        细心的人会发现，文献【8】的暑名单位为剑桥大学(University of Oxford)，并不是作者主页所在的南安普敦大学(University of Southampton)，在文献【8】的最后南提到：

Previous position?难道作者跳到Oxford了？

9、参考文献

【1】Blumensath T, Davies M E.Iterative hard thresholding for compressed sensing[J]. Applied & Computational HarmonicAnalysis, 2008, 27(3):265-274. (Available at:http://www.sciencedirect.com/science/article/pii/S1063520309000384)

【2】Blumensath T, Davies M E.Iterative Thresholding for Sparse Approximations[J]. Journal of Fourier Analysis & Applications,2008, 14(5):629-654. (Available at:http://link.springer.com/article/10.1007%2Fs00041-008-9035-z)

【3】Homepageof Blumensath T :http://www.personal.soton.ac.uk/tb1m08/index.html

【4】Lange, K., Hunter, D.R., Yang, I.. OptimizationTransfer Using Surrogate Objective Functions[J]. Journal of Computational &Graphical Statistics, 2000, 9(1):1-20. (Available at: http://sites.stat.psu.edu/~dhunter/papers/ot.pdf)

【5】GargR, Khandekar R. Gradient descent with sparsification: an iterative algorithmfor sparse recovery with restricted isometry property[C]//Proceedings of the26th Annual InternationalConference on Machine Learning. ACM, 2009: 337-344

【6】shasying2. landweber迭代方法.http://download.csdn.net/detail/shasying2/5092828

【7】Blumensath T, Davies M E.Normalized Iterative Hard Thresholding: Guaranteed Stability and Performance[J]. IEEE Journal of Selected Topics in Signal Processing, 2010,4(2):298-309.

【8】Blumensath T. Accelerated iterative hard thresholding[J]. Signal Processing, 2012, 92(3):752-756.