Image Deconvolution with the Half-quadratic Splitting Method

Image Deconvolution with the Half-quadratic Splitting Method

在处理图像重建或者逆问题的时候，我们经常会看到一种称为 Half-quadratic Splitting（HQS）的方法，这是在优化领域里非常经典的一种方法，之前也断断续续地找了很多相关的资料，发现斯坦福大学的计算成像课里某一节 lecture note 把这个方法在图像反卷积中的使用介绍地非常详细。这篇文章也是基于这个 lecture note 的一个学习笔记。

Image Formation

一般来说，图像的成像过程可以用下面的式子表示：

$$b=c\ast x+\eta \text{\hspace{1em}(1)}$$

其中，$b$ 表示观察到的模糊图像，$c$ 表示一个卷积核，可以认为是光学系统的 PSF，$\eta$ 表示与信号无关的噪声，$x$ 是希望重建出来的清晰图像。

根据信号处理的理论，在空域的卷积相当于在频域的乘积，所以上面的式子可以写成：

b=F−1{F(c)⋅F(x)}+η(2)b = \\mathcal{F}^{-1} \\{ \\mathcal{F}(c) \\cdot \\mathcal{F}(x) \\} + \\eta \\tag{2} b=F−1{F(c)⋅F(x)}+η(2)

不过要注意的是，这个等价只有在卷积满足 circular boundary conditions 时才成立。

为了后续的推导方便，这里将卷积运算转换成矩阵运算，可以得到如下的式子：

$$b=c\ast x\iff \mathbf{b}=\mathbf{Cx}\text{\hspace{1em}(3)}$$

Inverse Filtering and Wiener Deconvolution

在忽略噪声的情况下，式 (2) 可以表示成：

$$\mathcal{F}(b)=\mathcal{F}(c)\cdot \mathcal{F}(x)$$

进而可以直接算出 $x$:

x~if=F−1{F(b)F(c)}(4)\\tilde{x}_{if} = \\mathcal{F}^{-1} \\{ \\frac{\\mathcal{F}(b)}{\\mathcal{F}(c)} \\} \\tag{4} x~if​=F−1{F(c)F(b)​}(4)

上面的 if 就是 inverse filtering 的缩写，如果已知卷积核的形式，就可以计算出该卷积核对应的傅里叶变换，然后除以观察图像的傅里叶变换，从而得到清晰图像在频域的值，最后在做一个反傅里叶变换，得到最终的清晰图像。

这种方法看起来直观高效，不过有一个问题，就是 $\mathcal{F}(c)$ 的值很小趋近于 0 的时候，将会放大观察图像中的噪声，维纳滤波通过加入一个阻尼系数，可以避免这个问题：

x~wf=F−1{∣F(c)∣2∣F(c)∣2+1SNRF(b)F(c)}(5)\\tilde{x}_{wf} = \\mathcal{F}^{-1} \\left\\{ \\frac{|\\mathcal{F}(c)|^{2}}{|\\mathcal{F}(c)|^{2} + \\frac{1}{SNR}} \\frac{\\mathcal{F}(b)}{\\mathcal{F}(c)} \\right\\} \\tag{5} x~wf​=F−1{∣F(c)∣2+SNR1​∣F(c)∣2​F(c)F(b)​}(5)

SNR 表示信噪比，如果观测图像中没有噪声，那么 SNR 会很高，这种情况下，维纳滤波就变成了 inverse filter。

维纳滤波相比于直接的逆滤波，可以得到更加合理的反卷积结果。不过，这类方法面临的主要问题就是无法利用自然图像的先验信息，接下来，我们将介绍如何将先验信息与优化方法结合。

Nature Image Priors

首先，我们来看什么是自然图像先验，自然图像先验一般是一种数学模型，告诉我们在自然图像中，像素的分布更倾向于什么形态，在求解病态逆问题的时候，可能存在无数种解都符合观测值，而自然图像先验，会帮助我们挑选那些看起来更合理的解。

自然图像先验可以对图像的分布进行建模，这里，我们用正则化来表示对图像的某些性质，正则化一般对应图像先验的负对数。通常的正则化，包括平滑性，稀疏性，稀疏梯度，自相似性等。比如，稀疏性，可以用一阶范数表示为：∣x∣=∣∑xi∣|\\mathbf{x}| = |\\sum x_i|∣x∣=∣∑xi​∣，而平滑性，可以用梯度的二阶范数表示：∣Δx∣2|\\Delta \\mathbf{x}|_2∣Δx∣2​，不过，在求解图像复原的逆问题中，应用最广泛的还是称为 total variation (TV) 的正则，TV 的建立，是基于下面的观察，大部分自然图像，都可以看到很多区域都是平滑的，只有在不同物体的交界处才会出现边界，在同一个物体的内部，可以认为像素值区域平滑或者相似，而在不同物体的边界处，会有一个陡然的变化。

为了对 TV 建模，需要计算图像的一阶导数，这里，将图像水平方向的一阶导数表示为 Dxx\\mathbf{D}_x\\mathbf{x}Dx​x，而垂直方向的一阶导数表示为 Dyy\\mathbf{D}_y\\mathbf{y}Dy​y，具体的数学形式如下所示：

Dxx⇔x∗dx,dx=[0000−11000]\\mathbf{D}_x\\mathbf{x} \\Leftrightarrow \\mathbf{x} \\ast d_x, \\quad d_x = \\begin{bmatrix} 0 & 0 & 0\\\\ 0 & -1 & 1 \\\\ 0 & 0 & 0 \\end{bmatrix} Dx​x⇔x∗dx​,dx​=​000​0−10​010​​

Dyx⇔x∗dy,dy=[0000−10010]\\mathbf{D}_y\\mathbf{x} \\Leftrightarrow \\mathbf{x} \\ast d_y, \\quad d_y = \\begin{bmatrix} 0 & 0 & 0\\\\ 0 & -1 & 0 \\\\ 0 & 1 & 0 \\end{bmatrix} Dy​x⇔x∗dy​,dy​=​000​0−11​000​​

这个图像的一阶导数，既可以用矩阵直接计算，也可以用卷积的形式计算。

图像的 TV 正则可以表示为图像梯度的一阶范数，TV 正则项有两种表达形式，一种称为各向异性 的 TV，另外一种称为各向同性的正则，分别如下所示：

• 各向异性的 TV 正则

TVanisotropic(x)=∥Dxx∥1+∥Dyx∥1=∑i=1N∣(Dxx)i∣+∣(Dyx)i∣=∑i=1N(Dxx)i2+(Dyx)i2TV_{anisotropic}(\\mathbf{x}) = \\left \\| \\mathbf{D}_x\\mathbf{x} \\right \\|_{1} + \\left \\| \\mathbf{D}_y\\mathbf{x} \\right \\|_{1} = \\sum_{i=1}^{N} \\left| (\\mathbf{D}_x\\mathbf{x})_{i} \\right| + \\left| (\\mathbf{D}_y\\mathbf{x})_{i} \\right| = \\sum_{i=1}^{N} \\sqrt{(\\mathbf{D}_x\\mathbf{x})^{2}_{i}}+\\sqrt{(\\mathbf{D}_y\\mathbf{x})^{2}_{i}}TVanisotropic​(x)=∥Dx​x∥1​+∥Dy​x∥1​=i=1∑N​∣(Dx​x)i​∣+∣(Dy​x)i​∣=i=1∑N​(Dx​x)i2​​+(Dy​x)i2​​

• 各向同性的 TV 正则

TVisotropic(x)=∥Dx∥2,1=∑i=1N(Dxx)i2+(Dyx)i2TV_{isotropic}(\\mathbf{x}) = \\left \\| \\mathbf{D}\\mathbf{x} \\right \\|_{2,1} = \\sum_{i=1}^{N} \\sqrt{(\\mathbf{D}_x\\mathbf{x})^{2}_{i} + (\\mathbf{D}_y\\mathbf{x})^{2}_{i}} TVisotropic​(x)=∥Dx∥2,1​=i=1∑N​(Dx​x)i2​+(Dy​x)i2​​

这两种正则的差异，从表达式中可以看出，对于各向同性来说，每个像素 $i$ 的梯度约束是两个方向同时约束的，而对于各向异性来说，两个方向的梯度约束是分开约束的。

Regularized Deconvolution with Half-quadratic Splitting

The Half-quadratic Splitting Method

接下来，我们介绍 HQS 这种优化方法，在介绍用这种方法求解图像复原的逆问题之前，我们先探讨这种方法的一般形式，考虑如下的一个优化问题：

$$\mathbf{b}=\mathbf{Ax}+\eta$$

其中，x∈RN\\mathbf{x} \\in R^{N}x∈RN 是一个待求解的向量，b∈RM\\mathbf{b} \\in R^{M}b∈RM 表示一个观测向量，$\eta$ 表示噪声，A∈RM×N\\mathbf{A} \\in R^{M \\times N}A∈RM×N 表示转换矩阵，这个问题的一般求解形式如上所示：

min⁡x12∥Ax−b∥22+λΨ(x)\\min_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{A}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\lambda \\Psi(\\mathbf{x}) xmin​21​∥Ax−b∥22​+λΨ(x)

其中，前面的第一项一般称为数据项，第二项称为正则项，直接求解上面的式子，有的时候并不能很好地收敛，一种更为鲁棒的求解方式，应该是将上面的优化函数，改写成下面这种形式：

min⁡x12∥Ax−b∥22+λΨ(z)subject toDx−z=0\\min_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{A}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\lambda \\Psi(\\mathbf{z}) \\\\ \\text{subject to} \\quad \\mathbf{D}\\mathbf{x} - \\mathbf{z} = 0 xmin​21​∥Ax−b∥22​+λΨ(z)subject toDx−z=0

这个优化函数中，引入了一个中间变量 $\mathbf{z}$，这个额外的中间变量，可以将上面的优化函数拆成两部分，一部分是数据项，另外一部分是正则项，这两项依赖的变量是相互独立的，$\mathbf{x},\mathbf{z}$ 之间靠一个约束表达式联系，如果将这个约束项合入优化函数，则整个优化函数可以写成：

Lp(x,z)=f(x)+g(z)+p2∥Dx−z∥22L_p(\\mathbf{x}, \\mathbf{z}) = f(\\mathbf{x}) + g(\\mathbf{z}) + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} Lp​(x,z)=f(x)+g(z)+2p​∥Dx−z∥22​

优化函数写成这种形式，可以通过相互迭代的方式求解，求解 $\mathbf{x}$ 与 求解 $\mathbf{z}$ 可以分开进行：

x←proxf,p(z)=arg min⁡xLp(x,z)=arg min⁡xf(x)+p2∥Dx−z∥22(13)\\mathbf{x} \\gets \\mathbf{prox}_{f, p} (\\mathbf{z}) = \\argmin_{\\mathbf{x}} L_p(\\mathbf{x}, \\mathbf{z}) = \\argmin_{\\mathbf{x}} f(\\mathbf{x}) + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\tag{13} x←proxf,p​(z)=xargmin​Lp​(x,z)=xargmin​f(x)+2p​∥Dx−z∥22​(13)

z←proxf,p(Dx)=arg min⁡zLp(x,z)=arg min⁡zg(z)+p2∥Dx−z∥22(14)\\mathbf{z} \\gets \\mathbf{prox}_{f, p} (\\mathbf{D}\\mathbf{x}) = \\argmin_{\\mathbf{z}} L_p(\\mathbf{x}, \\mathbf{z}) = \\argmin_{\\mathbf{z}} g(\\mathbf{z}) + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\tag{14} z←proxf,p​(Dx)=zargmin​Lp​(x,z)=zargmin​g(z)+2p​∥Dx−z∥22​(14)

从上面的求解过程可以看出，当我们更新 $\mathbf{x}$ 的时候，只需要考虑 $f(\mathbf{x})$，而当我们更新 $\mathbf{z}$ 的时候，只需要考虑 $g(\mathbf{z})$，而不需要同时考虑这两个函数。这个性质，可以构建一个非常灵活的框架，让我们可以灵活地与各种正则函数相结合。这种方式也被称为 plug-and-play （即插即用）。

虽然 HQS 可以用于解决各种逆问题，不过我们这里还是讨论比较特殊的一种图像解卷积问题，我们讨论一种已知固定卷积核的情况，这样对应的矩阵是一个循环 Toeplitz 矩阵。先定义如下的关系：

c∗x=F−1{F(c)⋅F(x)}=CxF−1{F(c)∗⋅F(x)}=CTxF−1{F(b)F(c)}=C−1xc * x = \\mathcal{F}^{-1} \\{ \\mathcal{F}(c) \\cdot \\mathcal{F}(x) \\} = \\mathbf{C}\\mathbf{x} \\\\ \\mathcal{F}^{-1} \\{ \\mathcal{F}(c)^{*} \\cdot \\mathcal{F}(x) \\} = \\mathbf{C}^{T}\\mathbf{x} \\\\ \\mathcal{F}^{-1} \\{ \\frac{\\mathcal{F}(b)}{\\mathcal{F}(c)} \\} = \\mathbf{C}^{-1}\\mathbf{x} c∗x=F−1{F(c)⋅F(x)}=CxF−1{F(c)∗⋅F(x)}=CTxF−1{F(c)F(b)​}=C−1x

Standard Form of HQS with TV and Denoising Regularizers

接下来，我们考虑基于 TV 正则的 HQS 的优化方法，由上面的表达式，我们可以将带 TV 正则的优化函数写成如下形式：

min⁡x12∥Cx−b∥22+λΨ(z)subject toDx−z=0\\min_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{C}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\lambda \\Psi(\\mathbf{z}) \\\\ \\text{subject to} \\quad \\mathbf{D}\\mathbf{x} - \\mathbf{z} = 0 xmin​21​∥Cx−b∥22​+λΨ(z)subject toDx−z=0

其中，D=[DxT,DyT]∈R2N×N\\mathbf{D} = [\\mathbf{D}_{x}^{T}, \\mathbf{D}_{y}^{T}] \\in R^{2N \\times N}D=[DxT​,DyT​]∈R2N×N 表示 $x,y$ 方向的一阶导数，这里的 z∈R2N\\mathbf{z} \\in R^{2N}z∈R2N 比 x∈RN\\mathbf{x} \\in R^{N}x∈RN 要大一倍，因为每个像素，有 $x,y$ 两个方向的梯度。

对于更为一般的情况，我们可以使用一个简单的正则项，将待求解的图像投影到一个灵活的自然图像空间中，整个的 HQS 的形式可以写成如下所示：

min⁡x12∥Cx−b∥22+λΨ(z)subject tox−z=0\\min_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{C}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\lambda \\Psi(\\mathbf{z}) \\\\ \\text{subject to} \\quad \\mathbf{x} - \\mathbf{z} = 0 xmin​21​∥Cx−b∥22​+λΨ(z)subject tox−z=0

Efficient Implementation of the x-Update using Inverse Filtering

前面介绍过，HQS 的方法，会交替迭代更新 $\mathbf{x},\mathbf{z}$，我们先来看 $\mathbf{x}$ 的更新，

proxf,p(z)=arg min⁡xf(x)+p2∥Dx−z∥22=arg min⁡x12∥Cx−b∥22+p2∥Dx−z∥22(20)\\mathbf{prox}_{f, p} (\\mathbf{z}) = \\argmin_{\\mathbf{x}} f(\\mathbf{x}) + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} = \\argmin_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{C}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\tag{20} proxf,p​(z)=xargmin​f(x)+2p​∥Dx−z∥22​=xargmin​21​∥Cx−b∥22​+2p​∥Dx−z∥22​(20)

将上面的表达式展开，可以得到：

12∥Cx−b∥22+p2∥Dx−z∥22=12(Cx−b)T(Cx−b)+p2(Dx−z)(Dx−z)T=12(xTCTCx−2xTCTb+bTb)+p2(xTDTDx−2xTDTz+zTz)\\frac{1}{2} \\left \\| \\mathbf{C}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\frac{p}{2} \\left \\| \\mathbf{D}\\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\\\ = \\frac{1}{2}(\\mathbf{C}\\mathbf{x} - \\mathbf{b})^{T}(\\mathbf{C}\\mathbf{x} - \\mathbf{b}) + \\frac{p}{2}(\\mathbf{D}\\mathbf{x} - \\mathbf{z})(\\mathbf{D}\\mathbf{x} - \\mathbf{z})^{T} \\\\ = \\frac{1}{2}(\\mathbf{x}^{T}\\mathbf{C}^{T}\\mathbf{C}\\mathbf{x} - 2 \\mathbf{x}^{T}\\mathbf{C}^{T} \\mathbf{b} + \\mathbf{b}^{T}\\mathbf{b} ) + \\frac{p}{2}(\\mathbf{x}^{T}\\mathbf{D}^{T}\\mathbf{D}\\mathbf{x} - 2 \\mathbf{x}^{T}\\mathbf{D}^{T} \\mathbf{z} + \\mathbf{z}^{T}\\mathbf{z} ) 21​∥Cx−b∥22​+2p​∥Dx−z∥22​=21​(Cx−b)T(Cx−b)+2p​(Dx−z)(Dx−z)T=21​(xTCTCx−2xTCTb+bTb)+2p​(xTDTDx−2xTDTz+zTz)

将上面的表达式对 $\mathbf{x}$ 求导，可以得到：

CTCx−CTb+pDTDx−pDTz\\mathbf{C}^{T}\\mathbf{C}\\mathbf{x} - \\mathbf{C}^{T} \\mathbf{b} + p \\mathbf{D}^{T}\\mathbf{D}\\mathbf{x} - p \\mathbf{D}^{T} \\mathbf{z} CTCx−CTb+pDTDx−pDTz

让导数为 0 ，进而可以求得 $\mathbf{x}$：

(CTC+pDTD)−1(CTb+pDTz)(24)( \\mathbf{C}^{T}\\mathbf{C} + p\\mathbf{D}^{T}\\mathbf{D} )^{-1}(\\mathbf{C}^{T} \\mathbf{b} + p \\mathbf{D}^{T} \\mathbf{z} ) \\tag{24} (CTC+pDTD)−1(CTb+pDTz)(24)

对于满足 circular boundary conditions 的 2D 图像解卷积问题，上面的式子，可以变换到傅里叶域，然后再进行求解，上面所有的矩阵相乘的形式，都可以找到对应的傅里叶域的形式。

Special Case of TV Regularizer

如果正则项是 TV 项，那么 $\mathbf{D}$ 就是有限差分算子，上面的公式 (24) 可以写成如下形式：

(CTC+pDTD)⇔F−1{F{c}∗⋅F{c}+p(F{dx}∗⋅F{dx}+F{dy}∗⋅F{dy})}( \\mathbf{C}^{T}\\mathbf{C} + p\\mathbf{D}^{T}\\mathbf{D} ) \\Leftrightarrow \\mathcal{F}^{-1} \\{ \\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{c\\} + p (\\mathcal{F}\\{d_x\\}^{*} \\cdot \\mathcal{F}\\{d_x\\} + \\mathcal{F}\\{d_y\\}^{*} \\cdot \\mathcal{F}\\{d_y\\}) \\} (CTC+pDTD)⇔F−1{F{c}∗⋅F{c}+p(F{dx​}∗⋅F{dx​}+F{dy​}∗⋅F{dy​})}

(CTb+pDTz)⇔F−1{F{c}∗⋅F{b}+p(F{dx}∗⋅F{z1}+F{dy}∗⋅F{z2})}(\\mathbf{C}^{T} \\mathbf{b} + p \\mathbf{D}^{T} \\mathbf{z} ) \\Leftrightarrow \\mathcal{F}^{-1} \\{ \\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{b\\} + p (\\mathcal{F}\\{d_x\\}^{*} \\cdot \\mathcal{F}\\{z_1\\} + \\mathcal{F}\\{d_y\\}^{*} \\cdot \\mathcal{F}\\{z_2\\}) \\} (CTb+pDTz)⇔F−1{F{c}∗⋅F{b}+p(F{dx​}∗⋅F{z1​}+F{dy​}∗⋅F{z2​})}

由此，可以得到公式(20) 的解为：

proxf,p(z)=F−1(F{c}∗⋅F{b}+p(F{dx}∗⋅F{z1}+F{dy}∗⋅F{z2})F{c}∗⋅F{c}+p(F{dx}∗⋅F{dx}+F{dy}∗⋅F{dy}))\\mathbf{prox}_{f, p} (\\mathbf{z}) = \\mathcal{F}^{-1} \\left( \\frac{\\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{b\\} + p (\\mathcal{F}\\{d_x\\}^{*} \\cdot \\mathcal{F}\\{z_1\\} + \\mathcal{F}\\{d_y\\}^{*} \\cdot \\mathcal{F}\\{z_2\\})}{\\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{c\\} + p (\\mathcal{F}\\{d_x\\}^{*} \\cdot \\mathcal{F}\\{d_x\\} + \\mathcal{F}\\{d_y\\}^{*} \\cdot \\mathcal{F}\\{d_y\\})} \\right) proxf,p​(z)=F−1(F{c}∗⋅F{c}+p(F{dx​}∗⋅F{dx​}+F{dy​}∗⋅F{dy​})F{c}∗⋅F{b}+p(F{dx​}∗⋅F{z1​}+F{dy​}∗⋅F{z2​})​)

Special Case of Denoising Reg

对于更为一般的正则项，$\mathbf{D}$ 可以认为是一个单位矩阵，公式 (20) 的求解将变得更为简单：

proxf,p(z)=F−1(F{c}∗⋅F{b}+pF{z}F{c}∗⋅F{c}+p)\\mathbf{prox}_{f, p} (\\mathbf{z}) = \\mathcal{F}^{-1} \\left( \\frac{\\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{b\\} + p \\mathcal{F}\\{z\\} } {\\mathcal{F}\\{c\\}^{*} \\cdot \\mathcal{F}\\{c\\} + p } \\right) proxf,p​(z)=F−1(F{c}∗⋅F{c}+pF{c}∗⋅F{b}+pF{z}​)

Updating z with the TV Regularizer

公式(14) 关于 $\mathbf{z}$ 的更新，可以表示成如下的形式：

proxg,p(v)=Sλ/p(v)=arg min⁡zλ∣z∣1+p2∥v−z∥22\\mathbf{prox}_{g, p} (\\mathbf{v}) = \\mathcal{S}_{\\lambda/p}(\\mathbf{v}) = \\argmin_{\\mathbf{z}} \\lambda \\left| \\mathbf{z} \\right|_{1} + \\frac{p}{2} \\left \\| \\mathbf{v} - \\mathbf{z} \\right \\|_{2}^{2} proxg,p​(v)=Sλ/p​(v)=zargmin​λ∣z∣1​+2p​∥v−z∥22​

其中，$\mathbf{v}=\mathbf{Dx}$，$\mathcal{S}$ 是一个分段函数：

Sk(v)={v−kv>k0∣v∣<kv+kv<−k\\mathcal{S}_{k}(v) = \\left\\{\\begin{matrix} v - k & v > k \\\\ 0 & |v| < k \\\\ v + k & v < -k \\end{matrix}\\right. Sk​(v)=⎩⎨⎧​v−k0v+k​v>k∣v∣<kv<−k​

Isotropic TV Norm

对于各向同性的 TV 正则，正则项可以表示为：

g(z)=λ∥Dx∥2,1=λ∑i=1N(Dxx)i2+(Dyx)i2g(\\mathbf{z}) = \\lambda \\left \\| \\mathbf{D}\\mathbf{x} \\right \\|_{2,1} = \\lambda \\sum_{i=1}^{N} \\sqrt{(\\mathbf{D}_x\\mathbf{x})^{2}_{i} + (\\mathbf{D}_y\\mathbf{x})^{2}_{i}} g(z)=λ∥Dx∥2,1​=λi=1∑N​(Dx​x)i2​+(Dy​x)i2​​

那么，带各向同性正则项的解卷积问题可以表示成：

min⁡x12∥Cx−b∥22+λ∑i=1N∥[zizi+N]∥2subject toDx−z=0\\min_{\\mathbf{x}} \\frac{1}{2} \\left \\| \\mathbf{C}\\mathbf{x} - \\mathbf{b} \\right \\|_2^{2} + \\lambda \\sum_{i=1}^{N} \\left \\| \\begin{bmatrix} z_i\\\\ z_{i+N} \\end{bmatrix} \\right \\|_{2} \\\\ \\text{subject to} \\quad \\mathbf{D}\\mathbf{x} - \\mathbf{z} = 0 xmin​21​∥Cx−b∥22​+λi=1∑N​​[zi​zi+N​​]​2​subject toDx−z=0

ziz_izi​ 表示 $\mathbf{z}$ 的第 $i$ 个元素，其中，$1\le i\le N$ 表示水平方向的有限差分，$N+1\le i\le 2N$ 表示垂直方向的有限差分。

$\mathbf{z}$ 的更新可以表示成：

z←proxf,p(v)=arg min⁡zλ∑i=1N∥[zizi+N]∥2+p2∥v−z∥22v=Dx\\mathbf{z} \\gets \\mathbf{prox}_{f, p} (\\mathbf{v}) = \\argmin_{\\mathbf{z}} \\lambda \\sum_{i=1}^{N} \\left \\| \\begin{bmatrix} z_i\\\\ z_{i+N} \\end{bmatrix} \\right \\|_{2} + \\frac{p}{2} \\left \\| \\mathbf{v} - \\mathbf{z} \\right \\|_{2}^{2} \\quad \\mathbf{v} = \\mathbf{D}\\mathbf{x} z←proxf,p​(v)=zargmin​λi=1∑N​​[zi​zi+N​​]​2​+2p​∥v−z∥22​v=Dx

最终 $\mathbf{z}$ 的更新可以表示为：

[zizi+N]←Sλ/p([vivi+N]),1≤i≤N\\begin{bmatrix} z_i\\\\ z_{i+N} \\end{bmatrix} \\gets \\mathcal{S}_{\\lambda /p} \\left( \\begin{bmatrix} v_i\\\\ v_{i+N} \\end{bmatrix} \\right), \\quad 1 \\leq i \\leq N [zi​zi+N​​]←Sλ/p​([vi​vi+N​​]),1≤i≤N

Updating z with DnCNN or any Gaussian Denoiser as the Regularizer

如果我们进一步审视公式 (14) ，不考虑矩阵 $\mathbf{D}$ 的情况下，

arg min⁡zg(z)+p2∥x−z∥22=arg min⁡zλΨ(z)+p2∥x−z∥22=arg min⁡zΨ(z)+p2λ∥x−z∥22(36)\\argmin_{\\mathbf{z}} g(\\mathbf{z}) + \\frac{p}{2} \\left \\| \\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\\\ = \\argmin_{\\mathbf{z}} \\lambda \\Psi(\\mathbf{z}) + \\frac{p}{2} \\left \\| \\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\\\ = \\argmin_{\\mathbf{z}} \\Psi(\\mathbf{z}) + \\frac{p}{2 \\lambda} \\left \\| \\mathbf{x} - \\mathbf{z} \\right \\|_{2}^{2} \\tag{36} zargmin​g(z)+2p​∥x−z∥22​=zargmin​λΨ(z)+2p​∥x−z∥22​=zargmin​Ψ(z)+2λp​∥x−z∥22​(36)

公式 (36) 可以看成是一个降噪问题，可以等价成如下的表达式：

arg min⁡xΨ(x)+12σ2∥v−x∥22\\argmin_{\\mathbf{x}} \\Psi(\\mathbf{x}) + \\frac{1}{2 \\sigma^{2}} \\left \\| \\mathbf{v} - \\mathbf{x} \\right \\|_{2}^{2} xargmin​Ψ(x)+2σ21​∥v−x∥22​

其中，$\mathbf{v}$ 可以看成是观测到的有噪图像，$\mathbf{x}$ 表示我们需要恢复的无噪图像，基于这个观测，对于降噪的正则项，那么对 $\mathbf{z}$ 的更新可以简单地变成一个降噪过程，这个降噪可以用传统的降噪方法，也可以用基于CNN 的降噪方法，

proxD,p(x)=D(x,σ2=λp)\\mathbf{prox}_{\\mathcal{D}, p} (\\mathbf{x}) = \\mathcal{D} \\left( \\mathbf{x}, \\sigma^{2} = \\frac{\\lambda}{p} \\right) proxD,p​(x)=D(x,σ2=pλ​)

最后，做一个总结，基于 TV 正则和基于降噪正则的 HQS 方法的主要流程可以归纳如下：