斜风细雨做小寒,淡烟疏柳媚晴滩。入淮清洛渐漫漫。
雪沫乳花浮午盏,蓼茸蒿笋试春盘。人间有味是清欢。
---- 苏轼算法
更多精彩内容请关注微信公众号 “优化与算法”微信
迭代硬阈值(Iterative Hard Thresholding)算法是求解基于 \({\ell _0}\) 范数非凸优化问题的重要方法之一,在稀疏估计和压缩感知重构等领域应用较多。IHT最初由Blumensath, Thomas等人提出,后来许多学者在IHT算法的基础上不断发展出一些改进算法,如正规化迭代硬阈值算法(Normalized Iterative Hard Thresholding, NIHT)、共轭梯度硬阈值迭代算法(Conjugate Gradient Iterative Hard Thresholding, CGIHT)、基于回溯的硬阈值迭代算法(Backtracking based Iterative Hard Thresholding)等,以及后来与贪婪算法结合产生的一些算法,如硬阈值追踪算法(Hard Thresholding Pursuit, HTP)就是结合子空间追踪算法(SP)和IHT算法的优点获得的。
笔者在此将上述这些算法原理作简要分析,并附上算法和测试代码供参考。框架
以稀疏近似问题为例,介绍IHT算法。对于如下优化问题:函数
其中 \({\bf{y}} \in {{\bf{R}}^M}\), \({\bf{A}} \in {{\bf{R}}^{M \times N}}\),且 \(M < N\), \({\left\| {\bf{x}} \right\|_0}\) 为 \(\bf{x}\) 的 \({\ell _0}\) 范数,表示 \(\bf{x}\) 中非零元素的个数。因为 \(M < N\), 显然(1)式是一个欠定问题,加上 \({\ell _0}\) 约束,式(1)同时是个组合优化问题,没法在多项式时间复杂度内求解。测试
为了求解优化问题(1),Blumensath 等人提出了著名的 IHT 算法,其迭代格式为:优化
其中 \({H_S}( \cdot )\) 表示取 \(\cdot\) 中最大的 \(S\) 个非零元素。ui
式(2)看起来并不复杂,直观上来看就是对(1)式用常数步长的梯度降低法来求解,其中在每次迭代中对更新的结果进行一次硬阈值操做。为何硬阈值操做+梯度降低能够求解(1)式的 \({\ell _0}\) 范数优化问题呢?文献[1]中做者利用Majorization-Minimization优化框架推导出了(2)式所示迭代公式,并证实了IHT算法在 \({\left\| {\bf{A}} \right\|_2} < 1\) 时,算法能收敛到局部最小值。如下对IHT算法的原理作简要分析。spa
首先简单介绍一下MM优化框架,实际上MM优化框架是处理非凸优化问题中常见的一种方法,其形式简单,思想优雅。后期关于MM思想再以专题形式进行详细的总结和讨论。3d
(1) 式的目标函数为rest
(1) 没法获得一个闭式解,只能用迭代方法来求得一个近似解。
MM优化思想:当目标函数 \(F({\bf{x}})\) 比较难优化时,一般能够选择更容易优化的替代目标函数 \(G({\bf{x}})\),当 \(G({\bf{x}})\) 知足必定的条件时,\(G({\bf{x}})\) 的最优解可以无限逼近原目标函数 \(F({\bf{x}})\) 的最优解。这样的 \(G({\bf{x}})\) 须要知足3个条件:
在知足 \({\left\| {\bf{A}} \right\|_2} < 1\) 条件时,目标函数 (4) 可用以下替代目标函数来代替:
从 (4) 式能够看出,替代目标函数 \(G({\bf{x}},{\bf{z}})\) 是知足MM优化框架中所列条件的。
对式 (4) 展开并化简:
上式的最后三项是与 \(\bf{x}\) 无关的项,对优化结果没有影响,那么使 (5) 式最小化等价于下式最小化:
对 (5) 式进行配方,能够获得:
令 \({{\bf{x}}^*} = {{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} \left( {{\bf{y}} - {\bf{Az}}} \right)\),得:
即当 \({{{\bf{x}}_i} = {{\bf{x}}^*}}\) 时, \(G'({\bf{x}},{\bf{z}})\) 取得最小值,且最小值等于 \(\sum\limits_i - {\left( {{{\bf{x}}^*}} \right)^2}\) 。
结合式 (1) 中的 \({\ell _0}\) 范数约束项,最小化 \(G'({\bf{x}},{\bf{z}})\) 变成如何使得在 \(\bf{x}\) 中非零元素个数小于 \(S\) 的状况下,最小化 \(\sum\limits_i - {\left( {{{\bf{x}}^*}} \right)^2}\)。很容易想到,只须要保留 \(\bf{x}\) 中前 \(S\) 个最大项便可,这能够经过硬阈值函数来进行操做。
硬阈值函数以下所示:
其中 \(T\) 为 \(abs\left[ {{{\bf{x}}_k} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }({\bf{y}} - {\bf{A}}{{\bf{x}}_k})} \right]\) 按从大到小排列第 \(S\) 个值的大小。
硬阈值函数以下图所示,这跟前面介绍的软阈值函数相似。
function [hat_x] = cs_iht(y,A,K,itermax,error) % Email: zhangyanfeng527@gmail.com % Reference: Blumensath T, Davies M E. Iterative hard thresholding for compressed sensing[J]. % Applied & Computational Harmonic Analysis, 2009, 27(3):265-274. % y: measurement (observe) vector % A: measurement (sensing) matrix % k: sparsity level % 注意,IHT算法须要norm(A)<1,不然结果很不稳定,甚至不收敛。 u = 1 ; % 默认步长等于1,也能够本身选择 x0 = zeros(size(A,2),1) ; % initialization with the size of original for times = 1 : itermax x_increase = A' * (y - A * x0); hat_x = x0 + u * x_increase ; [~,pos] = sort(abs(hat_x),'descend'); hat_x(pos(K + 1 : end)) = 0 ; % thresholding, keeping the larges s elements if norm(y - A * hat_x) < error || norm(x0-hat_x)/norm(hat_x) < error break ; else x0 = hat_x ; % update end end
因为IHT算法须要知足 \({\left\| {\bf{A}} \right\|_2} < 1\),所以算法对 ${\left| {\bf{A}} \right|_2} $ 的缩放很敏感,当不知足此条件时,迭代过程会严重动荡甚至不收敛。此外,IHT算法固定步长为1,没法更好的利用自适应步长加速算法的收敛,而经过调节步长来实现算法收敛加速显然不现实。
基于以上缘由,Blumensath, Thomas 等人接下来又提出了正规化硬阈值迭代算法,NIHT算法的优点是利用李普希兹连续条件来规范化步长,使算法不受矩阵 \(\bf{A}\) 缩放的影响。同时采用一维线搜索,推导出自适应步长公式,使迭代次数大幅减小。这种方式以少许可不计的计算量增长换来了算法收敛速度的大幅提高。
NIHT算法具体的原理在此不细述,感兴趣的朋友能够参考文献[2],此处供上原文中的算法伪代码:
NIHT算法程序以下:
function [x_p] = cs_niht(y,A,K,itermax,error,par_c) % Email: zhangyanfeng527@gmail.com % Reference: 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. % y: measurement (observe) vector % A: measurement (sensing) matrix % k: sparsity level % itermax: max iteration % error: error threshold % par_c: 0 < par_c < 1 % NIHT算法的好处是使用自适应的步长,再也不像IHT那样要求norm(A)<1了 x0 = zeros(size(A,2),1); % initialization with the size of original g0 = A' * y ; u = norm(g0) / norm(A * g0) ; x1 = x0 + u * g0 ; [~,pos1] = sort(abs(x1),'descend') ; x1(pos1(K + 1:end)) = 0 ; gamma0 = find(x1) ; for times = 1 : itermax g1 = A' * (y - A * x1) ; % calculate gradient u = norm(g1(gamma0)) / norm(A(:,gamma0) * g1(gamma0)) ; % calculate step size x_p = x1 + u * g1 ; [~,pos1] = sort(abs(x_p),'descend') ; x_p(pos1(K + 1 : end)) = 0 ; gamma1 = find(x_p) ; % support set W = ( norm ( x_p - x1 ) )^2 / ( norm(A * (x_p - x1)) )^2 ; if isequal(gamma0,gamma1) == 1 elseif isequal(gamma0,gamma1) == 0 % normalized step size if u <= 1 * W elseif u > 1 * W while(u > 1 * W ) u = par_c * u ; end x_p = x1 + u * g1 ; [~,pos1] = sort(abs(x_p),'descend') ; x_p(pos1(K + 1 : end)) = 0 ; gamma1 = find(x_p) ; end end if norm(y - A * x_p) < error break; else x1 = x_p ; % update gamma0 = gamma1 ; end end
尽管NIHT算法在IHT算法的基础上改进了很多地方,可是总的来讲,IHT算法和NIHT算法都是基于梯度的算法,只利用了梯度域的信息,而梯度降低法自己收敛速度受限。为了进一步提升IHT算法的收敛速度,Blanchard, Jeffrey D等人提出使用共轭梯度来代替IHT算法中的梯度,并证实了CGIHT算法的收敛速度比IHT算法快。具体能够描述为,IHT算法和NIHT算法是基于梯度的算法,其收敛速度与
有关,而CGIHT算法使用共轭梯度,其收敛速度与
有关。其中 \(K\) 是矩阵 \({{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{A}}}\) 的条件数。很显而后者小于前者,所以其收敛速度更快。
CGIHT算法的原理比较简单,这里再也不详细介绍,下面是原文中CGIHT算法的伪代码:
CGIHT算法程序以下:
function [x_p] = cs_cgiht(y,A,s,itermax,error) % Emial: zhangyanfeng527@gmail.com % Reference: Blanchard J D, Tanner J, Wei K. CGIHT: conjugate gradient iterative hard thresholding % for compressed sensing and matrix completion[J]. % Information & Inference A Journal of the Ima, 2015(4). % y: measurement (observe) vector % A: measurement (sensing) matrix % s: sparsity level % itermax: max iteration % error: error threshold x_0 = .2 * ones(size(A,2),1); % initialization with the size of original g0 = A' * ( y - A * x_0) ; [~,pos0] = sort(abs(g0),'descend') ; g0(pos0(s+1:end)) = 0 ; d0 = g0 ; u = (norm(g0) / norm(A * d0))^2 ; x_1 = x_0 + u * d0 ; [~,pos1] = sort(abs(x_1),'descend'); x_1(pos1(s + 1 : end)) = 0 ; gamma1 = find(x_1) ; %% restart CGIHT algorithm % tao = gamma1 ; % for times = 1 : itermax % g1 = A' * ( y - A * x_1) ; % if isequal(gamma0,tao) == 1 % u = (g1(gamma1)' * g1(gamma1)) / (g1(gamma1)' * A(:,gamma1)' * A(:,gamma1) * g1(gamma1)) ; % x_p = x_1 + u * g1 ; % else % belta = (norm(g1) / (norm(g0)))^2; % d1 = g1 + belta * d0 ; % u = (d1(gamma1)' * g1(gamma1)) / (d1(gamma1)' * A(:,gamma1)' * A(:,gamma1) * d1(gamma1)) ; % x_p = x_1 + u * d1 ; % end % [~,index] = sort(abs(x_p),'descend'); % x_p(index(s + 1 : end)) = 0 ; % if norm(x_p - x_1)/norm(x_p)< error || norm(y - A * x_p)<error % break; % else % gamma0 = gamma1 ; % gamma1 = index ; % x_1 = x_p ; % end % end %% CGIHT algorithm for times = 1 : itermax g1 = A' * ( y - A * x_1) ; belta = (norm(g1) / (norm(g0)))^2 ; d1 = g1 + belta * d0 ; u = (d1(gamma1)' * g1(gamma1)) / (d1(gamma1)' * A(:,gamma1)' * A(:,gamma1) * d1(gamma1)) ; x_p = x_1 + u * d1 ; [~,index] = sort(abs(x_p),'descend') ; x_p(index(s + 1 : end)) = 0 ; if norm(x_p - x_1)/norm(x_p)< error || norm(y - A * x_p)<error break; else x_1 = x_p ; g0 = g1 ; d0 = d1 ; end end end
硬阈值追踪算法(Hard Thresholding Pursuit HTP)是由SIMON FOUCART提出来的,该算法结合了子空间追踪算法(Subspace Pursuit SP)算法和IHT算法的优点,一方面具有更强的理论保证,另外一方面也表现出良好的实证效果。
与IHT算法不一样的是,HTP算法在每次迭代中采用了SP算法中用最小二乘更新稀疏系数的方法,最小二乘是一种无偏估计方法,这使得更新结果更加精确,从而能够减小迭代次数。
与SP算法不一样的是,SP算法在每次迭代中是采用测量矩阵与残差的内积(即\({{{\bf{A}}^T}({\bf{y}} - {\bf{Ax}})}\))的大小来挑选最大的 \(S\) 个支撑位置,而HTP算法则是像IHT算法同样采用 \({{\bf{x}}{\rm{ + }}{{\bf{A}}^T}({\bf{y}} - {\bf{Ax}})}\) 来挑选支撑,后者无疑比前者包含更多的信息,所以相对于SP算法,HTP算法能更快的找到支撑位置,从而减小迭代次数。
总的说来,HTP算法是结合了IHT算法和SP算法的优点,所以其表现也更好。HTP算法原理较为简单,在此不详细介绍,对SP算法的介绍在后续贪婪类算法中一并总结。
下面是HTP算法的伪代码:
HTP算法的程序以下:
function [hat_x] = cs_htp(y,A,K,itermax,error) % Emial: zhangyanfeng527@gmail.com % Reference: Foucart S. HARD THRESHOLDING PURSUIT: AN ALGORITHM FOR COMPRESSIVE SENSING[J]. % Siam Journal on Numerical Analysis, 2011, 49(6):2543-2563. % y: measurement (observe) vector % A: measurement (sensing) matrix % k: sparsity level x0 = .2 * ones(size(A,2),1) ; % initialization with the size of original u = 1 ; % stpesize=1 for times = 1 : itermax x_increase = A' * (y - A * x0); hat_x = x0 + u * x_increase ; [~,pos] = sort(abs(hat_x),'descend'); gamma = pos(1:K) ; z = A(:,gamma) \ y ; hat_x = zeros(size(A,2),1) ; hat_x(gamma) = z ; if norm(y - A * hat_x) < error || norm(hat_x-x0)/norm(hat_x) < error break; else x0 = hat_x ; end end end
如下对上述四种算法进行重构效果比较,测试程序以下:
%% CS reconstruction test % Emial: zhangyanfeng527@gmail.com clear clc close all 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稀疏的,且位置是随机的 Psi = eye(N); % x自己是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta Phi = randn(M,N); %测量矩阵为高斯矩阵 Phi = orth(Phi')'; A = Phi * Psi; %传感矩阵 % sigma = 0.005; % e = sigma*randn(M,1); % y = Phi * x + e;%获得观测向量y y = Phi * x ; %获得观测向量y %% 恢复重构信号x tic theta = cs_iht(y,A,K); x_r = Psi * theta;% x=Psi * theta toc %% 绘图 figure; plot(x_r,'k.-');%绘出x的恢复信号 hold on; plot(x,'r');%绘出原信号x hold off; legend('Recovery','Original') fprintf('\n恢复残差:'); norm(x_r-x)%恢复残差
上述测试脚本只提供了简单的实验设置,对算法重构成功率、算法迭代次数与偏差的关系(收敛速度)、重构时间、鲁棒性等测试能够自由发挥。
[1] Blumensath, Thomas, and Mike E. Davies. "Iterative hard thresholding for compressed sensing." Applied and computational harmonic analysis 27.3 (2009): 265-274.
[2] Blumensath, Thomas, and Mike E. Davies. "Normalized iterative hard thresholding: Guaranteed stability and performance." IEEE Journal of selected topics in signal processing 4.2 (2010): 298-309.
[3] Blanchard, Jeffrey D., Jared Tanner, and Ke Wei. "CGIHT: conjugate gradient iterative hard thresholding for compressed sensing and matrix completion." Information and Inference: A Journal of the IMA 4.4 (2015): 289-327.
[4] Foucart, Simon. "Hard thresholding pursuit: an algorithm for compressive sensing." SIAM Journal on Numerical Analysis 49.6 (2011): 2543-2563.
更多精彩内容请关注微信公众号 “优化与算法”