迭代软阈值(IST)算法解决CS优化目标函数

本文主要是介绍迭代软阈值(IST)算法解决CS优化目标函数，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

　　彻底理解了迭代硬阈值IHT以后，很自然的会想到：如果将软阈值(SoftThresholding)函数与Majorization-Minimization优化框架相结合形成迭代软阈值(Iterative Soft Thresholding，IST)算法（另一种常见简称为ISTA，即IterativeSoftThresholding Algorithm，另外Iterative有时也作Iterated），应该可以解决如下优化问题[1]：

此即基追踪降噪(Basis Pursuit De-Noising, BPDN)问题。

1、迭代软阈值算法流程

为了解决优化问题（将原a换为习惯的x）

执行迭代算法

其中Ψλ是对每一个数值计算软阈值函数值

这个算法称为迭代软阈值(Iterative Soft Thresholding，IST)算法。

2、迭代软阈值算法推导

基追踪降噪(Basis Pursuit De-Noising, BPDN)问题的目标函数

由于目标函数f(x)并不容易优化，根据Majorization-Minimization优化框架，可以优化替代目标函数

这里要求||Φ||2<1，则有

下面，我们对替代目标函数进行变形化简

其中，是与x无关的项；所以优化替代目标函数u(x,z)时，可以等价于优化

其中。看到这个问题熟悉么？

对！这正是软阈值(SoftThresholding)函数要解决的以下优化问题

对于标准的软阈值(Soft Thresholding)函数来说，这个优化问题的解是

注意：这里的B是一个向量，应该是逐个元素分别执行软阈值函数；其中标准的软阈值(SoftThresholding)函数是：

将符号换为此处的优化问题

则解为

注意：这里的x*是一个向量，应该是逐个元素分别执行软阈值函数；其中

然后我们根据Majorization-Minimization优化框架的流程进行迭代即可，注意z代表xn，即当前迭代值，而优化解soft(x*,λ)代表xn+1，用于下次迭代。

由于这个算法的整个过程相当于迭代执行软阈值(SoftThresholding)函数，所以把它称为迭代软阈值(Iterative Soft Thresholding)算法。

3、迭代软阈值算法MATLAB代码

在IST函数中，一共规定了三种IST跳出迭代的条件：第一个是重构结果经Phi变换后与原先y的差异，第二是最大迭代次数，第三个是重构结果x两次相邻迭代的差异。

以下为2016-08-12更新版本V1.1，相比于原先的V1.0改动之处为第1个跳出迭代循环的条件参考TwIST作了修改：

 ? 
 function [ x ] = IST_Basic( y,Phi,lambda,epsilon,loopmax )   
 %   Detailed explanation goes here     
 %Version: 1.0 written by jbb0523 @2016-08-09    
 %Version: 1.1 modified by jbb0523 @2016-08-12    
      if nargin < 5       
          loopmax = 10000;       
      end      
      if nargin < 4         
          epsilon = 1e-2;         
      end       
      if nargin < 3         
          lambda = 0.1* max ( abs (Phi'*y));        
      end       
      [y_rows,y_columns] =  size (y);         
      if y_rows<y_columns         
          y = y'; %y should be a column vector         
      end   
      soft = @(x,T)  sign (x).* max ( abs (x) - T,0);   
      n =  size (Phi,2);       
      x =  zeros (n,1); %Initialize x=0     
      f = 0.5*(y-Phi*x)'*(y-Phi*x)+lambda* sum ( abs (x)); %added in v1.1   
      loop = 0;    
      fprintf ( '\n' );   
      while 1       
          x_pre = x;   
          x = soft(x + Phi'*(y-Phi*x),lambda); %update x           
          loop = loop + 1;     
          f_pre = f; %added in v1.1   
          f = 0.5*(y-Phi*x)'*(y-Phi*x)+lambda* sum ( abs (x)); %added in v1.1   
          if abs (f-f_pre)/f_pre<epsilon %modified in v1.1   
              fprintf ( 'abs(f-f_pre)/f_pre<%f\n' ,epsilon);   
              fprintf ( 'IST loop is %d\n' ,loop);   
              break ;   
          end  
          if loop >= loopmax   
              fprintf ( 'loop > %d\n' ,loopmax);   
              break ;   
          end  
          if norm (x-x_pre)<epsilon   
              fprintf ( 'norm(x-x_pre)<%f\n' ,epsilon);   
              fprintf ( 'IST loop is %d\n' ,loop);   
              break ;   
          end  
      end    
 end 

　　原先的V1.0版本：

 ? 
 function [ x ] = IST_Basic( y,Phi,lambda,epsilon,loopmax )   
 %   Detailed explanation goes here     
 %Version: 1.0 written by jbb0523 @2016-08-09    
      if nargin < 6       
          loopmax = 10000;       
      end      
      if nargin < 5         
          epsilon = 1e-6;         
      end       
      if nargin < 4         
          lambda = 0.1* max ( abs (Phi'*y));        
      end       
      [y_rows,y_columns] =  size (y);         
      if y_rows<y_columns         
          y = y'; %y should be a column vector         
      end   
      soft = @(x,T)  sign (x).* max ( abs (x) - T,0);   
      n =  size (Phi,2);       
      x =  zeros (n,1); %Initialize x=0     
      loop = 0;    
      fprintf ( '\n' );   
      while 1       
          x_pre = x;   
          x = soft(x + Phi'*(y-Phi*x),lambda); %update x           
          loop = loop + 1;     
          if norm (y-Phi*x)<epsilon   
              fprintf ( 'norm(y-Phi*x)<%f\n' ,epsilon);   
              fprintf ( 'IST loop is %d\n' ,loop);   
              break ;   
          end  
          if loop >= loopmax   
              fprintf ( 'loop > %d\n' ,loopmax);   
              break ;   
          end  
          if norm (x-x_pre)<epsilon   
              fprintf ( 'norm(x-x_pre)<%f\n' ,epsilon);   
              fprintf ( 'IST loop is %d\n' ,loop);   
              break ;   
          end  
      end    
 end 

　　该函数非常简单，核心程序就一行，其它均为配角，为这一行服务的：

 ? 
 x = soft(x + Phi'*(y-Phi*x),lambda); %update x   

4、迭代软阈值算法测试

这里首先按传统的测试程序进行重构测试，这里的λ取值参考了文献【2】作者主页给出的代码里的方法(可以自己试一下，这个值到底是1还是2其实影响不大)：

 ? 
 clear all ; close all ; clc ;           
 M = 64; %观测值个数           
 N = 256; %信号x的长度           
 K = 10; %信号x的稀疏度           
 Index_K =  randperm (N);           
 x =  zeros (N,1);           
 x(Index_K(1:K)) = 5* randn (K,1); %x为K稀疏的，且位置是随机的    
 %x(Index_K(1:K)) = sign(5*randn(K,1));                
 Phi =  randn (M,N); %测量矩阵为高斯矩阵       
 Phi =  orth (Phi ')' ;               
 sigma = 0.005;         
 e = sigma* randn (M,1);       
 y = Phi * x + e; %得到观测向量y           
 % y = Phi * x;%得到观测向量y         
 %% 恢复重构信号x           
 tic  
 % lamda = sigma*sqrt(2*log(N));   
 lamda = 0.1* max ( abs (Phi'*y));   
 fprintf ( '\nlamda = %f\n' ,lamda)   
 %x_r =  BPDN_quadprog(y,Phi,lamda);    
 x_r = IST_Basic(y,Phi,lamda);             
 toc    
 %% 绘图           
 figure ;           
 plot (x_r, 'k.-' ); %绘出x的恢复信号           
 hold on;           
 plot (x, 'r' ); %绘出原信号x           
 hold off;           
 legend ( 'Recovery' , 'Original' )           
 fprintf ( '\n恢复残差：' );           
 fprintf ( '%f\n' , norm (x_r-x));%恢复残差   

运行结果如下：（信号为随机生成，所以每次结果均不一样）

1)图

2）CommandWindows:

lamda= 0.177639

norm(x-x_pre)<0.000001

ISTloop is 106

Elapsedtime is 0.008842 seconds.

恢复残差：1.946232

从图中可以看出，重构结果的位置基本都是正确的，但确比原始信号都要小一些，为什么呢？从Command Windows的输出结果来看，IST迭代的106次，跳出迭代的条件是x在相邻两次迭代中基本不变了（norm(x-x_pre)<0.000001），也就是说最优点x已经基本不变了，若同样执行使用MATLAB自带的quadprog的基追踪降噪BPDN_quadprog函数（参见压缩感知重构算法之基追踪降噪），重构结果是正常的，到底为什么本文的IST重构结果得不到近似的x呢？难道理论推导有问题？

直到看了文献【2】中的Debiasing：

文中提到：“In many situations, itis worthwhile to debias the solution as a postprocessing step, to eliminate theattenuation of signal magnitude due to the presence of the regularization term.”，直译过来就是“在很多场合，作为一个后处理步骤对结果进行除偏(debias)是很值得的，用于消除由于正则项的存在对信号幅度造成的衰减”，读到这里就明白了，为什么本文的IST重构结果总是与真实的x幅度差一些呢？这是因为IST求解的优化问题含有正则项(regularizationterm)，即λ||x||1。直观地来讲，因为目标函数中存在λ||x||1项，在最优化min时这个也要尽量的小，所以IST的恢复的结果比实际的x幅度要小一些（与参数λ有关）。解决这个问题的办法就是对结果进行除偏(debias，注：有道词典查不到此单词，也没看中文文献是如何翻译的，直接根据词根琢磨了一下自己瞎翻译的)。

如何除偏(debias)呢？接着来看：“Inthe debiasing step, we fix at zero those individual components or groups thatare zero at the end of the SpaRSA process, and minimize the objective over theremaining elements.”，翻译一下就是“将SpaRSA重构结果中的为零的项固定为零，然后在剩除项上对目标函数进行最小化”，简单一点说就是优化文中的式(27)（正好与本文IST解决的问题相同，后来会专门分析SpaRSA算法）：

其中AI的含义更清晰的解释如下：

若SpaRSA重构的x为：

原来的矩阵A为：

则AI等于原矩阵A只保留第3列（对应x中的元素a）和第5列（对应x中的元素b）的子矩阵：

文献【2】后面继续谈到，若要求解式(27)可以使用共轭梯度法(Conjugate gradient)。实际上，式(27)的最优解为最小二乘解（熟悉匹配追踪系列算法的应该对这个很清楚，可参见压缩感知中的数学知识：投影矩阵（projectionmatrix）），即

在SpaRSA算法中，作者给出的官方代码里包括了debias过程，由一个参数来控制是否对重构结果执行debias过程，这里为了简单，就不修改IST函数了，直接在测试代码中加入debias代码：

 ? 
 clear all ; close all ; clc ;           
 M = 64; %观测值个数           
 N = 256; %信号x的长度           
 K = 10; %信号x的稀疏度           
 Index_K =  randperm (N);           
 x =  zeros (N,1);           
 x(Index_K(1:K)) = 5* randn (K,1); %x为K稀疏的，且位置是随机的    
 %x(Index_K(1:K)) = sign(5*randn(K,1));                
 Phi =  randn (M,N); %测量矩阵为高斯矩阵       
 Phi =  orth (Phi ')' ;               
 sigma = 0.005;         
 e = sigma* randn (M,1);       
 y = Phi * x + e; %得到观测向量y           
 % y = Phi * x;%得到观测向量y         
 %% 恢复重构信号x           
 tic  
 % lamda = sigma*sqrt(2*log(N));   
 lamda = 0.1* max ( abs (Phi'*y));   
 fprintf ( '\nlamda = %f\n' ,lamda)   
 %x_r =  BPDN_quadprog(y,Phi,lamda);    
 x_r = IST_Basic(y,Phi,lamda);             
 toc    
 %Debias   
 [xsorted inds] =  sort ( abs (x_r),  'descend' );    
 AI = Phi(:,inds(xsorted(:)>1e-3));   
 xI =  pinv (AI '*AI)*AI' *y;   
 x_bias =  zeros ( length (x),1);   
 x_bias(inds(xsorted(:)>1e-3)) = xI;   
 %% 绘图           
 figure ;           
 plot (x_r, 'k.-' ); %绘出x的恢复信号           
 hold on;           
 plot (x, 'r' ); %绘出原信号x           
 hold off;           
 legend ( 'Recovery' , 'Original' )           
 fprintf ( '\n恢复残差(original)：' );           
 fprintf ( '%f\n' , norm (x_r-x));%恢复残差     
 %Debias   
 figure ;           
 plot (x_bias, 'k.-' ); %绘出x的恢复信号           
 hold on;           
 plot (x, 'r' ); %绘出原信号x           
 hold off;           
 legend ( 'Recovery-debise' , 'Original' )           
 fprintf ( '恢复残差(debias)：' );           
 fprintf ( '%f\n' , norm (x_bias-x));%恢复残差     

运行结果如下：（信号为随机生成，所以每次结果均不一样）

1)图：

第一幅图（IST重构结果与原信号对比）：

第二幅图（对IST重构结果debias后的结果与原信号对比）：

2)CommandWindows：

lamda = 0.405767

norm(x-x_pre)<0.000001

IST loop is 74

Elapsed time is 0.007175 seconds.

恢复残差(original)：3.736317

恢复残差(debias)：0.847947

从第二幅图的结果来看，debias后的结果已与BPDN_quadprog函数结果（参见压缩感知重构算法之基追踪降噪）在同一个水平了。

值得注意的是，在求最小二乘解的公式里包含矩阵求逆的过程，矩阵求逆在大规模(large scale)问题中是一个比较麻烦的问题，IST等迭代算法（包括IHTs等）本身相比于MP等算法的优势即为不含矩阵求逆，计算效率高，因此实际中求解文献【2】的式(27)不应该直接采用如上包含矩阵求逆的方式，尤其是大规模问题，否则IST算法计算量小的优势就没有了，此处仅是为了示范debias的效果。

5、结束语

值得注意的是，一般文献中出现的IST或ISTA简称中的“S”并非指的是“soft”，而是“shrinkage”，即“IteratedShrinkage/ThresholdingAlgorithm”，至于“shrinkage”是什么意思呢？我们还是先来看一下有道词典吧：

那么Iterative Soft Thresholding和Iterated Shrinkage/Thresholding有什么区别呢？你要认为它们是一样的也没问题，因为从文献中来看，很多作者的确是这么认为的；另外，还可以认为Iterative Soft Thresholding是Iterated Shrinkage/Thresholding的一种特殊情况，即每次迭代时正好是求软阈值函数时的特殊情况，而Iterated Shrinkage/Thresholding更是一种广义的称呼。

值得注意的是，本文参考文献很少，且均是与IST“不相关”的，这是因为没有哪一种文献明确提出IST，所以也不知道引用哪个了。

在下一篇中，我们从多篇文献中来看一看究竟什么是“Shrinkage”……

IST：Iterative Shrinkage/Thresholding和Iterative Soft Thresholding

6、参考文献

【1】Chen S S, Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.

【2】WrightS J, Nowak R D, Figueiredo M A T. Sparse reconstruction by separable approximation.[J]. IEEE Transactions on Signal Processing, 2009,57(7):3373-3376.

这篇关于迭代软阈值(IST)算法解决CS优化目标函数的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！