徹底理解了迭代硬閾值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.