迭代硬閾值類算法總結||IHT/NIHT/CGIHT/HTP

迭代硬閾值類（IHT）算法總結

---

斜風細雨作小寒，淡煙疏柳媚晴灘。入淮清洛漸漫漫。
雪沫乳花浮午盞，蓼茸蒿筍試春盤。人間有味是清歡。
---- 蘇軾

更多精彩內容請關注微信公眾號 “優化與算法”

迭代硬閾值（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算法的優勢得到的。
筆者在此將上述這些算法原理做簡要分析，並附上算法和測試代碼供參考。

1. 迭代硬閾值算法

1.1 IHT算法描述

以稀疏近似問題為例，介紹IHT算法。對於以下優化問題：

\[\mathop {\min }\limits_{\bf{x}} \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2{\rm{ }} ~~~~~~s.t.~~~~~{\rm{ }}{\left\| {\bf{x}} \right\|_0} \le S~~~~~~~~~~~~(1) \]

其中 \({\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 算法，其迭代格式為：

\[{{\bf{x}}_{k + 1}} = {H_k}({{\bf{x}}_k} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }({\bf{y}} - {\bf{A}}{{\bf{x}}_k}))~~~~~~~~~~~~~~~~(2) \]

其中 \({H_S}( \cdot )\) 表示取 \(\cdot\) 中最大的 \(S\) 個非零元素。

式(2)看起來並不復雜，直觀上來看就是對（1）式用常數步長的梯度下降法來求解，其中在每次迭代中對更新的結果進行一次硬閾值操作。為什么硬閾值操作+梯度下降可以求解(1)式的 \({\ell _0}\) 范數優化問題呢？文獻[1]中作者利用Majorization-Minimization優化框架推導出了（2）式所示迭代公式，並證明了IHT算法在 \({\left\| {\bf{A}} \right\|_2} < 1\) 時，算法能收斂到局部最小值。以下對IHT算法的原理做簡要分析。

1.2 Majorization-Minimization優化框架

首先簡單介紹一下MM優化框架，實際上MM優化框架是處理非凸優化問題中常見的一種方法，其形式簡單，思想優雅。后期關於MM思想再以專題形式進行詳細的總結和討論。

(1) 式的目標函數為

\[{F}({\bf{x}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 ~~~~~~~~~~~~~~~~~~~~~~~~~~~(3) \]

(1) 無法得到一個閉式解，只能用迭代方法來求得一個近似解。

MM優化思想：當目標函數 \(F({\bf{x}})\) 比較難優化時，通常可以選擇更容易優化的替代目標函數 \(G({\bf{x}})\)，當 \(G({\bf{x}})\) 滿足一定的條件時，\(G({\bf{x}})\) 的最優解能夠無限逼近原目標函數 \(F({\bf{x}})\) 的最優解。這樣的 \(G({\bf{x}})\) 需要滿足3個條件：

1. 容易優化
2. 對所有 \(\bf{x}\)，\(G({\bf{x}}) \ge F({\bf{x_{k}}})\)
3. \({G_k}({{\bf{x}}_k}) = F({{\bf{x}}_k})\)
   

在滿足 \({\left\| {\bf{A}} \right\|_2} < 1\) 條件時，目標函數 (4) 可用如下替代目標函數來代替：

\[G({\bf{x}},{\bf{z}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 - \left\| {{\bf{Ax}} - {\bf{Az}}} \right\|_2^2 + \left\| {{\bf{x}} - {\bf{z}}} \right\|_2^2~~~~~~~(4) \]

從 (4) 式可以看出，替代目標函數 \(G({\bf{x}},{\bf{z}})\) 是滿足MM優化框架中所列條件的。

對式 (4) 展開並化簡：

\[\eqalign{ & G({\bf{x}},{\bf{z}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 - \left\| {{\bf{Ax}} - {\bf{Az}}} \right\|_2^2 + \left\| {{\bf{x}} - {\bf{z}}} \right\|_2^2 \cr & = \left( {\left\| {\bf{y}} \right\|_2^2 - 2{{\bf{x}}^{\mathop{\rm T}\nolimits} }{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} + \left\| {{\bf{Ax}}} \right\|_2^2} \right) - \left( {\left\| {{\bf{Ax}}} \right\|_2^2 - 2{{\bf{x}}^{\mathop{\rm T}\nolimits} }{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}} + \left\| {{\bf{Az}}} \right\|_2^2} \right) + \left( {\left\| {\bf{x}} \right\|_2^2 - 2{{\bf{x}}^T}{\bf{z}} + \left\| {\bf{z}} \right\|_2^2} \right) \cr & = \left[ {\left\| {\bf{x}} \right\|_2^2 - \left( {2{{\bf{x}}^T}{\bf{z}} + 2{{\bf{x}}^T}{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} - 2{{\bf{x}}^T}{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}}} \right)} \right] + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr & = \left[ {\left\| {\bf{x}} \right\|_2^2 - 2{{\bf{x}}^T}\left( {{\bf{z}} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} - {{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}}} \right)} \right] + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr & = \sum\limits_i {\left[ {{\bf{x}}_i^2 - 2{{\bf{x}}_i}\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr} \]

上式的最后三項是與 \(\bf{x}\) 無關的項，對優化結果沒有影響，那么使 (5) 式最小化等價於下式最小化：

\[G'({\bf{x}},{\bf{z}}) = \sum\limits_i {\left[ {{\bf{x}}_i^2 - 2{{\bf{x}}_i}\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} ~~~~~~~~~~~~~~~~~~~~(5) \]

對 (5) 式進行配方，可以得到：

\[\eqalign{ & G'({\bf{x}},{\bf{z}}) = \sum\limits_i {\left[ {{\bf{x}}_i^2 - 2{{\bf{x}}_i}\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} \cr & = \sum\limits_i {\left[ {{\bf{x}}_i^2 - 2{{\bf{x}}_i}\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right) + \left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right) - \left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} \cr & = \sum\limits_i {\left\{ {{{\left[ {{\bf{x}}_i^2 - \left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]}^2} - {{\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)}^2}} \right\}} \cr} \]

令 \({{\bf{x}}^*} = {{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} \left( {{\bf{y}} - {\bf{Az}}} \right)\)，得：

\[\sum\limits_i {\left[ {{\bf{x}}_i^2 - \left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} = \sum\limits_i {\left\{ {{{\left[ {{{\bf{x}}_i} - {{\bf{x}}^*}} \right]}^2} - {{\left( {{{\bf{x}}^*}} \right)}^2}} \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\) 個最大項即可，這可以通過硬閾值函數來進行操作。

1.3 硬閾值函數

硬閾值函數如下所示：

\[{H_S}({{\bf{x}}_i}) = \left\{ {\matrix{ {0{\rm{~~~~~if~~~}}\left| {{{\bf{x}}_i}} \right| < T} \cr {{{\bf{x}}_i}{\rm{~~~if~~~}}\left| {{{\bf{x}}_i}} \right| \ge T} \cr } } \right.\]

其中 \(T\) 為 \(abs\left[ {{{\bf{x}}_k} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }({\bf{y}} - {\bf{A}}{{\bf{x}}_k})} \right]\) 按從大到小排列第 \(S\) 個值的大小。

硬閾值函數如下圖所示，這跟前面介紹的軟閾值函數類似。

1.4 IHT算法程序

function [hat_x] = cs_iht(y,A,K,itermax,error)  
% Email: [e-mail]
% 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  

2. NIHT算法

由於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: [e-mail]
% 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  

3. CGIHT算法

盡管NIHT算法在IHT算法的基礎上改進了不少地方，但是總的來說，IHT算法和NIHT算法都是基於梯度的算法，只利用了梯度域的信息，而梯度下降法本身收斂速度受限。為了進一步提高IHT算法的收斂速度，Blanchard, Jeffrey D等人提出使用共軛梯度來代替IHT算法中的梯度，並證明了CGIHT算法的收斂速度比IHT算法快。具體可以描述為，IHT算法和NIHT算法是基於梯度的算法，其收斂速度與

\[{{K - 1} \over {K + 1}} \]

有關，而CGIHT算法使用共軛梯度，其收斂速度與

\[{{\sqrt K - 1} \over {\sqrt K + 1}} \]

有關。其中 \(K\) 是矩陣 \({{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{A}}}\) 的條件數。很顯然后者小於前者，因此其收斂速度更快。

CGIHT算法的原理比較簡單，這里不再詳細介紹，下面是原文中CGIHT算法的偽代碼：

CGIHT算法程序如下：

function [x_p] = cs_cgiht(y,A,s,itermax,error)    
% Emial: [e-mail]
% 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

4. HTP 算法

硬閾值追蹤算法（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: [e-mail]
% 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

5. 仿真

以下對上述四種算法進行重構效果比較，測試程序如下：

%% CS reconstruction test
% Emial: [e-mail]

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.

更多精彩內容請關注微信公眾號 “優化與算法”