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

压缩感知重构算法之迭代硬阈值(Iterative Hard Thresholding,IHT)转载自:https://blog.csdn.net/wyw921027/article/details/52102211题目:压缩感知重构算法之迭代硬阈值(IterativeHardThresholding,IHT)本篇来介绍IHT重构算法。一般在压缩感知参考文献中,提到IHT时一般引用的都是文献【1】,但IHT实际上是在文献【2】中提出的。IHT并不是一种凸优化算法,…

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

转载自: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应用于压缩感知重构问题时进行了一个理论分析:

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

 

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

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

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

 

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

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

2、Step1:替代目标函数

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

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

这里中压缩感知重构算法之迭代硬阈值(Iterative Hard Thresholding,IHT)的M应该指的是M-sparse,S应该指的是Surrogate。这里要求:

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

 

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

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

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

 

将目标函数定义为:

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

 

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

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

 

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

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

 

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

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

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

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

 

的道理是一模一样的。

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

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

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

        接下来,式(3.5)进行了变形:

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

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

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

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

4、Step3:极值点的获得

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

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

 

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

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

 

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

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

 

        当压缩感知重构算法之迭代硬阈值(Iterative Hard Thresholding,IHT),取得最小值

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

5、Step4:迭代公式的获得

        极值点得到了,替代目标函数的极小值也得到了:

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

 

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

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

 

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

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

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

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

6、IHT算法的MATLAB代码

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

第一个版本:

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

 

[plain] view plain copy

  在CODE上查看代码片派生到我的代码片

  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

  在CODE上查看代码片派生到我的代码片

  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

  在CODE上查看代码片派生到我的代码片

  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

  在CODE上查看代码片派生到我的代码片

  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稀疏。可见算法的名称是也是一不断完善的过程啊……

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

 

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】中的一段话:

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

 

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

8.4Landweber迭代

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

 

[plain] view plain copy

  在CODE上查看代码片派生到我的代码片

  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]:

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

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

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

 

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

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

        RIHT的算法流程如下:

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

 

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

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

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

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

 

8.6稀疏度对IHT的影响

 

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

8.7作者去向?

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

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

 

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.

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/144144.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • Dataset之CIFAR-10:CIFAR-10数据集的简介、下载、使用方法之详细攻略

    Dataset之CIFAR-10:CIFAR-10数据集的简介、下载、使用方法之详细攻略Dataset之CIFAR-10:CIFAR-10数据集的简介、下载、使用方法之详细攻略目录CIFAR-10的简介1、与MNIST数据集中目比,CIFAR-10真高以下不同点2、TensorFlow官方示例的CIFAR-10代码文件3、CIFAR-10数据集的数据文件名及用途4、基于CIFAR-10数据集最新算法预测准确率对比CIFAR-10的下载1、下载CIFAR-10数据集的全部数据CIFAR-10使用方法1、使用TF读取CIFAR-10数据官网链接:TheCIFAR-10datas

    2022年10月17日
  • Java + Ajax跨域解决方案整理

    Java + Ajax跨域解决方案整理为什么会跨域呢?简单来说就是前端页面与后台服务没有部署在同一个服务器上。产生跨域的情况有:1.域名不同,端口也不同;2.域名相同但是端口不同;3.域名不同,端口相同。解决方案:一、JSONP方式1.只支持get方法,不支持postfang方法;使用时需修改前端和后端代码,用起来也不太方便,本文不推荐使用。二、使用springMVC架构的,使用版本4.2以上…

  • C++中的空类默认包含哪些类成员函数

    C++中的空类默认包含哪些类成员函数

  • 关于dbunit报Duplicate entry ‘????’ for key ‘xxx’错误的问题

    关于dbunit报Duplicate entry ‘????’ for key ‘xxx’错误的问题

  • burp suite抓包教程

    burp suite抓包教程一、打开burpsuite,点击上方Proxy,再点击openbrowser并在地址栏内输入www.baidu.com二、抓包完成后,在action菜单栏中点击SendtoRepeater三、点击Repeater,再点击Send四、点击Render,右侧出现百度页面,完成抓包流程…

  • 小程序bindtap传参_微信小程序bindtap

    小程序bindtap传参_微信小程序bindtap一边开发一边做点笔记,东西可能零散了点,一边开发一边补充。1、事件 1.bindtap绑定点击事件 2.bindinput监听输入,没输入一个字符得到一次返回值(就算是输入中文时,没敲一次键依然返回一次)2、解决小程序tabBar跳转不能带参数问题小程序这里遇到了一个难题就是如果实现tabBar栏之间的跳转的话是不能传入参数的那么我们要如何解决这个问题呢! 我的办法就是让你的传…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号