圖象恢復學習筆記（二）

主要關於ISTA和ADMM。

一些常用的矩陣和向量微分

矩陣A-標量t

　　　　　　

標量函數f-矢量x

　　　　　　

矢量函數g-矢量x

　　　　　　

常用導數（b,x為向量，A為矩陣）

　　　　　　

　

ISTA(Iterative Shrinkage-Thresholding Algorithm)

ISTA算法作用是求解以下形式目標函數

 　　　　　　

其中

前一項為最小二乘數據擬合項，這一部分是可微的，可以用簡單的梯度下降求解；后一項為L1范數懲罰項，作用是得到稀疏解。

ISTA的方法即將梯度下降的迭代解轉換為

　　　　　　

的形式，然后用軟閾值方法求解。

梯度下降迭代解

考慮到如果沒有懲罰項，梯度下降迭代解可以看作是使f的局部二次模型達到最小值時x的取值：

　　　　　　

其中t_k為步長。

加上懲罰項之后

　　　　　　

令，轉換為

　　　　　　

考慮到

上式轉換為

 　　　　　　

軟閾值

考慮函數

如何求f(x)為最小值時x的取值？

f(x)中后一項|x|在0點不可微，但是凸函數，可以使用次梯度（subgradient）處理，|x|的次導數為sgn(x)

　　　　　　

對f(x)求次梯度，並令次梯度為0：

　　　　　　　　　　

考慮tk=1，λ=1的簡單情況

　　　　　　　　　　

將y看作x*的函數：

　　　　　　　　　　

旋轉坐標軸：

　　　　　　　　　

軟閾值函數即：

　　　　　

ISTA算法

綜合以上兩個部分，目標函數

　　　　　　　　

可以通過以下式子更新迭代求解：

　　　　　　　　

其中T是軟閾值  

　　　　　　

ADMM最優化方法

主要思想是采用分裂的方法求解以下方程：

 　　　　　　

Chan S H , Wang X , Elgendy O A . Plug-and-Play ADMM for Image Restoration: Fixed Point Convergence and Applications[J]. IEEE Transactions on Computational Imaging, 2016, 3(1):84-98.

這篇文章介紹了ADMM（交替方向乘子）方法，這是一種在圖像恢復中廣泛應用的約束優化方法，主要優點在於具有模塊化的結果，可以作為子算法模塊插入現成的圖像降噪算法中，但收斂條件和速度還要看具體算法的實現情況。

MAP

最大后驗概率（MAP）：

　　　　　　

（1）與下式等價：

　　　　　　

其中

這是個無約束優化問題（unconstrained optimization）,可以用ADMM求解。

ADMM

ADMM的主要思想是通過分裂變量將無約束優化問題轉化為約束問題，然后交替迭代求解。

將變量x分裂為x和v，x=v，上式轉化為：

　　　　

其增廣拉格朗日函數為：

　　　　

其中，增設一個變量u，稱為拉格朗日算子（又稱對偶變量），添加了懲罰函數（ρ 項）。

根據ADMM，根據以下步驟重復迭代，交替對x,v,u進行更新：

　　　　　　

收斂到增廣拉格朗日函數的鞍點。

在交替迭代的步驟中，（5）可以看作一個反演過程，涉及前向成像模型f(x)，（6）可以看做降噪噪過程，涉及到先驗g(v)。

令σ = √ λ/ρ,（6）變為：

　　　　　　

將v(k)項視為含噪聲圖像，（8）最小化無噪聲圖像v和v(k)之間的殘差，如果先驗g(v)為全變差范數（total variation norm），（8）為標准全變差去噪問題。

ADMM舉例-lasso

http://web.stanford.edu/~boyd/papers/admm/

目標函數

　　　　　　　　

寫成 ADMM 形式

　　　　　　　　

其中　　　　

更新步驟：

　　　　　　　

Matlab代碼：

測試： randn('seed', 0); rand('seed',0); m = 1500;       % number of examples n = 5000;       % number of features p = 100/n;      % sparsity density x0 = sprandn(n,1,p); A = randn(m,n); A = A*spdiags(1./sqrt(sum(A.^2))',0,n,n); % normalize columns
b = A*x0 + sqrt(0.001)*randn(m,1); lambda_max = norm( A'*b, 'inf' );
lambda = 0.1*lambda_max;

程序代碼：
function [z, history] = lasso_via_ADMM(A, b, lambda, rho, alpha)
 
% Solve lasso problem via ADMM
% objective function: 
%       minimize 1/2*|| Ax - b ||_2^2 + \lambda || x ||_1
% The solution is returned in the vector x.
% rho is the augmented Lagrangian parameter.
% alpha is the over-relaxation parameter (typical values for alpha are
% between 1.0 and 1.8).
 
%%%   Global constants and defaults     %%%%
MAX_ITER = 1000;
ABSTOL   = 1e-4;
RELTOL   = 1e-2;
 
%%%%   Data preprocessing    %%%%
[m, n] = size(A);
Atb = A'*b;
 
%%%%    ADMM solver    %%%%
x = zeros(n,1);
z = zeros(n,1);
u = zeros(n,1);
 
for k = 1:MAX_ITER
    
    % x-update
    [m,n] = size(A);
    temp = Atb + rho*(z - u);  
    if(m<n)
        L = chol( speye(m) + 1/rho*(A*A'), 'lower' )
        L = sparse(L);
        U = sparse(L');  
        x = temp/rho - (A'*(U \ ( L \ (A*temp) )))/rho^2;
    else
        L = chol( A'*A + rho*speye(n), 'lower' );
        L = sparse(L);
        U = sparse(L');
        x = U \ (L \ q);
    end
    
    % z-update with relaxation
    zold = z;
    x_hat = alpha*x + (1 - alpha)*zold;
    z = soft_threshold(x_hat + u, lambda/rho);
 
    % u-update
    u = u + (x_hat - z);
 
    % diagnostics, reporting, termination checks
    r_norm = norm(x - z);
    s_norm  = norm(-rho*(z - zold));
    eps_pri = sqrt(n)*ABSTOL + RELTOL*max(norm(x), norm(-z));
    eps_dual = sqrt(n)*ABSTOL + RELTOL*norm(rho*u);
 
    if (r_norm < eps_pri && s_norm < eps_dual)
         break;
    end
end
end
 
%soft_threshold
function z = soft_threshold(x,kappa)
    z = sign(x).*max(abs(x)-kappa,0);  
end