斜風細雨做小寒,淡煙疏柳媚晴灘。入淮清洛漸漫漫。
雪沫乳花浮午盞,蓼茸蒿筍試春盤。人間有味是清歡。
---- 蘇軾算法
更多精彩內容請關注微信公衆號 「優化與算法」微信
迭代硬閾值(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.
更多精彩內容請關注微信公衆號 「優化與算法」