$${\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