TV正则化：半二次方分裂法

$${\bf{O}} = \mathop {\arg \min }\limits_{\bf{O}} \left|\left| {\bf{O} – \bf{S}} \right|\right|_2^2 + {\left|\left| {\nabla \bf{O}} \right|\right|_1}$$

2008年Wang Yilun等在SIAM Journal on Imaging Sciences发表了一个题为A new alternating minimization algorithm for total variation image reconstruction的工作。该工作中提出了一种使用半二次方分裂（Half-quadratic splitting，HQS）的方法。这篇文章将介绍如何使用HQS求解TV正则化。

TV正则化需要求解

${\bf{O}} = \mathop {\arg \min }\limits_{\bf{O}} \left|\left| {\bf{O} – \bf{S}} \right|\right|_2^2 + {\left|\left| {\nabla \bf{O}} \right|\right|_1} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ (1)

因为TV项：${\left|\left| {\nabla \bf{O}} \right|\right|_1}$的存在使得直接求解原问题变得困难。使用半二次方分裂法，引入一个额外变量${\bf{G}} = \nabla {\bf{O}}$替换原式的梯度项，并增加一个额外的数据保真项：${\lambda _0}\| {\nabla {\bf{O}} -{\bf{G}}}\|_2^2$，${\lambda _0} \to + \infty$ 将原问题简化为以下两个子问题分别求解：

${{\bf{O}} = \mathop {\arg \min }\limits_{\bf{O}} \| {{\bf{O}} – {\bf{S}}} \|_2^2 + {\lambda _0}\| {\nabla {\bf{O}} – {\bf{G}}} \|_2^2}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ (2)

${{\bf{G}} = \mathop {\arg \min }\limits_{\bf{G}} {\lambda _0}\| {\nabla {\bf{O}} – {\bf{G}}} \|_2^2 + \lambda {{\| {\bf{G}} \|}_1}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ (3)

其中${\lambda _0}$是一个非常大的正实数以实现新的数据保真项的约束。

关于$\bf{G}$的子问题（3）可以用软阈值求解，其解析解为：

${\bf{G}} = {\rm{sign}}\left( {\nabla {\bf{O}}} \right) \circ \max \left( {\left| {\nabla {\bf{O}}} \right| – \lambda /{\lambda _0},0} \right)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ (4)

关于$\bf{O}$的子问题（2）是纯二次函数型，可以令（2）对$\bf{O}$的导数为0求解得到$\bf{O}$为：

${\bf{O}} = \frac{{{\bf{S}} + {\lambda _0}{\nabla ^T}\nabla {\bf{G}}}}{{1 + {\lambda _0}{\nabla ^T}\nabla }}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $(5)

公式（5）可以使用傅里叶变换快速求解。

如上文提到${\lambda _0}$必须是一个非常大的正实数以保证数据保真项。但是如果一开始就给${\lambda _0}$赋值一个很大的数，收敛速度不是特别理想。因此Wang Yilun等提出了一种让${\lambda _0}$逐步增大的迭代优化思路。HQS求解的完整过程可以列为如下形式：

Half-quadratic splitting methods for TV-regularization

Initializing: ${\lambda}$, ${\lambda _0}={\lambda }$, $a$, ${\lambda _{\max }}$. ${\bf{O}} =\bf{S}$.
When ${\lambda _0}<{\lambda _{\max }}$ do
$\ \ \ \ \ \ \ \ $Update $\bf{G}$ using Eq. (4)
$\ \ \ \ \ \ \ \ $Update $\bf{O}$ using Eq. (5)
$\ \ \ \ \ \ \ \ $${\lambda _0} = a{\lambda _0}$
End When

MATLAB代码如下

% input lambda = 0.01; image: s
% output: denoised/smoothed image o

dx = psf2otf([-1,1],size(s)); % x-direction differential operator
dy = psf2otf([-1;1],size(s)); % y-direction differential operator

lambda0 = lambda;
lambda_max = 1e5;

while lambda0 < lambda_max
   %% g sub-problem
    gx = real(ifft2(fft2(o).*dx));
    gy = real(ifft2(fft2(o).*dy));
        
    % Anisotropy TV
    gx = sign(gx) .* max(abs(gx) - lambda/lambda0,0);
    gy = sign(gy) .* max(abs(gy) - lambda/lambda0,0);

    % isotropy TV
    ss = sqrt(gx.^2 + gy.^2);
    gx = gx./(ss + eps) .* max(ss - lambda/lambda0,0);
    gy = gy./(ss + eps) .* max(ss - lambda/lambda0,0);
        
    %% o sub=problem
    Numerator =   fft2(s) + lambda0 * (fft2(gx).*conj(dx) + fft2(gy).*conj(dy));
    Denominator = 1 + lambda0 * DTD;
        
    o = real(ifft2(Numerator./Denominator));
    lambda0 = lambda0 * 2;
end