半二次图像复原中结构化方程组的预处理方法研究

A Review of Preconditioners for Structured Systems Arising from Half-Quadratic Image Restoration

摘要: 半二次正则化最小二乘是实现高质量图像复原的重要模型之一。在利用牛顿迭代方法等优化方法求解该模型的过程中，每一步都涉及结构化方程组的求解。预处理共轭梯度法(PCG)是求解此类方程组的有效方法，而其收敛速度取决于预处理后矩阵的特征值性质。构造合适的预处理矩阵对于提高图像复原的性能具有重要的意义。近年来，结合半二次图像复原中方程组的结构化特点，学者们基于矩阵的Schur补近似等策略构造出了一系列的预处理矩阵，并给出了相应的特征值分析。数值结果表明，这些预处理方法有效地降低了图像复原的计算成本。针对半二次图像复原中的结构化方程组，本文整理了近几年出现的预处理方法，并从不同侧面进行对比分析，旨在为进一步的预处理方法改进和研究提供思路参考。

Abstract: Half-quadratic regularized least square method is one of the important models for achieving high-quality image restoration. In the process of solving this model using optimization methods such as Newton method, each step involves solving a structured system. The Preconditioned Conju-gate Gradient (PCG) method is an efficient method for solving such systems while its convergence rate depends on the nature of the eigenvalues of the preconditioner. Constructing a suitable pre-conditioner is meaningful to improving the performance of image restoration. In recent years, com-bined with the structured characteristics of the system in half-quadratic image restoration, scholars have constructed a series of preconditioners based on the Schur complement approximation of the matrices and other strategies, and the corresponding eigenvalue analyses are given. The numerical results show that these preconditioners effectively reduce the computational cost of image restora-tion. For the structured system in half-quadratic image restoration, this paper organizes the pre-conditioners that appeared in recent years and makes a comparative analysis from different sides, aiming to provide ideas reference for further improvement and research of preconditioners.

1. 问题背景

图像复原是对退化图像进行处理以恢复图像原始信息的过程，在生物医学、天文成像等许多工程和科学领域发挥重要作用 [1] [2] [3] [4] 。在图像形成过程中，原始图像经常被噪声和模糊污染。因此，为了从退化图像中恢复出理想的图像，需要进行图像复原，常见的方法有：图像去噪、平滑、填补和去模糊等等 [4] [5] 。近年来，基于最小二乘正则化模型的代数方法逐渐成为实现高质量图像复原的主流方法 [6] [7] 。

1.1. 图像复原

由于引起图像退化的因素与性质各不相同，所以图像复原是一个复杂的数学过程，本质上是一个典型求解不适定的反问题 [8] 。对于一个 $n\times n$ 的图像X，其退化可以建模为

$b=Ax+\eta ,$ (1)

其中向量 $x\in R^{n^2}$ 和 $b\in R^{n^2}$ 分别是由原始图像X与观测图像B按行拉直而成，向量 $\eta \in R^{n^2}$ 为高斯白噪声， $A\in R^{n^2\times n^2}$ 为空间不变模糊矩阵。图像复原的主要任务是基于一些先验信息，从观测图像B中估计原始图像X。除了滤波方法以外，通常可以通过求解最小化问题来实现

$\hat{x}=\underset{x\in R^{n^2}}{\min }J\left(x\right)\text{.}$ (2)

图像复原的效果很大程度上取决于所建立的数学模型。即使对同一幅图像，利用不同数学模型求出的结果可能会有较大的差别。建立好模型后，图像复原就成了数学上反问题的求解。最简单的求解方法就是用最小二乘法，此时的代价函数为

$J\left(x\right)={\Vert Ax-b\Vert }_2^2\text{.}$ (3)

这里，b指的是带有噪音的退化图像，x指的是将要被复原的图像。上式目的是不让被复原的图像偏离退化图像太多，但这可能会造成图像的过度拟合。因为没有约束项，模型得到的结果也可能不唯一。因此，最小二乘法是一个不适定的方法，通常需要对其添加具有约束作用的正则化项以优化图像恢复模型。

1.2. 优化模型

正则化是缓解反问题不适定性、约束解特征的重要方式。对模型添加惩罚项 [9] [10] ，得到要求解的最小化问题

$\hat{x}=\arg \underset{x\in R_n^2}{\min }J\left(x\right)=\arg \underset{x\in R_n^2}{\min }\left\{{\Vert Ax-b\Vert }_2^2+\beta R\left(x\right)\right\},$ (4)

其中， ${\Vert Ax-b\Vert }_2^2$ 为数据保真项，用于保留原始图像的主要信息，正则项为形式 $R\left(x\right)={\displaystyle {\sum }_{i=1}^r\varphi \left(g_i^{\text{T}}x\right)}$ ，其中 $g_i^{\text{T}}:R^{n^2}\to R$ 对于 $i=1,\cdots ,r$ 是线性算子，G为 $g_i^{\text{T}}$ 的 $r\times n^2$ 矩阵，我们考虑 $\varphi \cancel{\equiv }0$ 为(连续可微函数)凸保边势函数的情况 [9] [11] [12] ，满足 ${\varphi }^{″}\left(t\right)>0$ ， $A^{\text{T}}A$ 可逆且 $\ker \left(A^{\text{T}}A\right)\cap \ker \left(G^{\text{T}}G\right)=\left\{0\right\}$ 。 $\beta$ 是平衡这两项的正则化参数。

不同的R项对应不同的正则化方法。典型的吉洪诺夫(Tikhonov)正则化 [13] [14] 能缓解反问题的不适定性，但通常会使结果变得光滑，不利于刻画图像边界。稀疏正则化 [15] [16] 在边界连续性和边缘结构上的刻画更具有优势，因为其L1范数能使结果具有稀疏表征的效果，更符合实际画质情况。全变差(TV)正则化 [17] [18] 则在实现噪声压制的同时，保留图像原本的结构特征，具有恢复原始信号边缘的良好性能，但其高非线性和不可微性给计算带来挑战。而半二次(HQ)正则化 [9] [11] [19] 不止能很好地保留恢复图像中的图像细节，也能有效降低最小化问题的计算成本，其正则项一般有加性和乘性两种形式，本文将着重研究半二次正则化的图像复原模型。

1.3. 预处理

图像复原中通常使用牛顿方法求解最小化问题，而牛顿法的每一次迭代都需要求解一个方程组。不同正则化模型对应不同特点的结构化方程组。这些方程组的系数矩阵一般是对称正定的，因此通常选择预处理共轭梯度法(PCG)来进行求解 [20] [21] 。然而，方程组收敛性质取决于特征值分布情况，当系数矩阵的谱分布不够聚集时，Krylov子空间法的收敛速度会相当慢。为了解决这个问题，需要对线性方程组进行预处理。近几十年来，方程组预处理一直是许多领域活跃的研究课题 [22] [23] [24] [25] 。

不同的预处理方法取决于方程组的构造和对系数矩阵的近似，对于半二次图像复原中的结构化方程组，学者们通常基于系数矩阵的不完全分解和Schur补近似来构造不同的预处理矩阵，并结合到PCG法中求解该线性方程组。针对半二次图像复原中出现的结构化方程组，本文介绍了近些年来提出的不同预处理方法，并进行了对比分析和总结。

本文具体结构如下：在第2部分，介绍了半二次正则化最小二乘的基本模型，并分别给出了其加性形式和乘性形式对应的结构化方程组。在第3部分，整理并对比分析了关于这些方程组的预处理方法。最后，本文总结了半二次图像复原中的预处理方法，并为进一步的方法改进和研究提出了建议和参考思路。

2. 半二次(HQ)正则化模型

2.1. 模型描述

最小化问题(4)通常具有很高的计算复杂度。此时，通过给正则化项添加辅助变量的方法，可以将代价函数 $J\left(x\right)$ 转化为其增广形式 [26] ，以此求解来降低优化问题的计算成本。

设辅助变量为 $v\in R^r$ ，则增广代价函数 $\tilde{J}:R^{n^2}\times R^r\to R$ 的形式

$\tilde{J}\left(x,\nu \right)={\Vert Ax-b\Vert }_2^2+\beta {\displaystyle {\sum }_{i=1}^r\left(Q\left(g_i^{\text{T}}x,\nu \left(i\right)\right)+\psi \left(\nu \left(i\right)\right)\right)},$ (5)

其中 $Q\left(.,s\right):R\to R$ 对于任意 $s\in R$ 是二次的， $\psi \left(s\right)$ 是函数 $\varphi \left(t\right)$ 的对偶势函数， $\psi :R\to R$ 满足

$\varphi \left(t\right)=\underset{s\in R}{\min }\left\{Q\left(t,S\right)+\psi \left(s\right)\right\},\forall t\in R,$ (6)

当势函数 $\varphi \left(t\right)$ 给定时，利用凸共轭理论立即可以确定对偶势函数 $\psi \left(s\right)$ ，条件(6)保证了

$J\left(x\right)=\underset{v\in R^r}{\min }\tilde{J}\left(x,\nu \right),\forall x\in R^{n^2},$ (7)

从而将(4)中的最小化问题转化为

$\left(\hat{x},\hat{v}\right)=\underset{x\in R^{n^2},v\in R^r}{\min }\tilde{J}\left(x,\nu \right)\text{.}$ (8)

对于不同的问题，二次函数Q可以取两种形式，分别是加性形式和乘性形式。我们将对两种形式下对应图像复原方程组的预处理矩阵进行研究。

2.2. 加性形式中的方程组

Q的加性形式为 [26] ，

$Q\left(t,s\right)=\frac{1}{2}{\left(t-s\right)}^2,t\in R,s\in R,$ (9)

此时，使用牛顿法求解(8)中的最小化问题，其迭代形式为，

$\left[\begin{array}{c}{x}_{k+1}\\ {\nu }_{k+1}\end{array}\right]=\left[\begin{array}{c}{x}_k\\ {\nu }_k\end{array}\right]+d_k,$ (10)

在每次牛顿迭代中，我们需要求解一个2 × 2块的结构化线性方程组，其中d_k是线性方程组的解，

$H\left(x_k,{\nu }_k\right)d_k=-\nabla \tilde{J}\left(x_k,{\nu }_k\right),$ (11)

这里 $\tilde{J}\left(x_k,{\nu }_k\right)$ 的Hessian矩阵如下，

$H\left(x_k,{\nu }_k\right)=\left[\begin{array}{cc}2A^{\text{T}}A+\beta G^{\text{T}}G& -\beta G^{\text{T}}\\ -\beta G& \beta \Lambda \left({\nu }_k\right)\end{array}\right]=\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right],$ (12)

其中 $\Lambda \left({\nu }_k\right)$ 是一个对角矩阵，其对角元为 ${\psi }^{″}\left({\nu }_i\right)>0,i=1,\cdots ,r$ 。

2.3. 乘性形式中的方程组

Q项的乘法形式为 [9] [27] ，

$Q\left(t,s\right)=\frac{1}{2}{t}^{2}s,t\in R,s\in R_{+},$ (13)

在每次牛顿迭代中，我们需要求解一个2×2块的结构化线性方程组

$H\left(x_k,{\nu }_k\right)d_k=\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{c}{d}_1\\ d_2\end{array}\right]=\left[\begin{array}{c}{r}_1\\ r_2\end{array}\right]\equiv r,$ (14)

其中d是线性方程组的解。对应 $\tilde{J}\left(x_k,{\nu }_k\right)$ 的Hessian矩阵为，

$H\left(x_k,{\nu }_k\right)=\left[\begin{array}{cc}2A^{\text{T}}A+\beta G^{\text{T}}\Lambda \left({\nu }_k\right)G& 2\beta G^{\text{T}}\Lambda \left(Gx\right)\\ -2\beta \Lambda \left(Gx\right)G& \beta \Lambda \left({\psi }^{″}\left({\nu }_k\right)\right)\end{array}\right]\equiv \left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right],$ (15)

其中 $\Lambda \left({\nu }_k\right)$ 是对角矩阵， $H\equiv H\left(x_k,{\nu }_k\right)$ 是对称矩阵， $H_{11},H_{22}$ 是对称正定的，并满足 $A^{\text{T}}A$ 是可逆的且 ${\varphi }^{″}\left(t\right)>0$ 对于 $\forall t\in R$ 成立。

3. 预处理方法

针对半二次图像复原模型的加性形式和乘性形式，本部分关于其中的结构化方程组，分别整理了相应的预处理方法，并从不同侧面进行了对比分析。

3.1. 加性形式

3.1.1. 预处理形式及性质

关于方程组(11)，Huang和Lu [11] 提出了两个块预处理矩阵。第一个块预处理矩阵为，

$M_1=\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}+H_{21}{H}_{11}^{-1}{H}_{12}\end{array}\right],$ (16)

它具有以下的分解形式，

$M_1=\left[\begin{array}{cc}I& 0\\ H_{21}{H}_{11}^{-1}& I\end{array}\right]\left[\begin{array}{cc}{H}_{11}& 0\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}I& H_{11}^{-1}{H}_{12}\\ 0& I\end{array}\right]\text{.}$ (17)

理论分析显示了预处理矩阵 $M_1$ 的谱性质。

定理1 预处理矩阵 $H{M}_1^{-1}$ 至少有 $n^2$ 个单位特征值和 $n^2$ 个单位特征值对应的线性无关的特征向量。 $H{M}_1^{-1}$ 其他的特征值分布在 $\left[1-\delta ,1\right]$ ，其中 $\delta =\max {\left(1+{\psi }^{″}\left(v_k\right)\right)}^{-1}$ 。

类似地，第二个块预处理矩阵构建如下，

$M_2=\left[\begin{array}{cc}{H}_{11}+H_{12}{H}_{22}^{-1}{H}_{21}& H_{12}\\ H_{21}& H_{22}\end{array}\right],$ (18)

具有以下分解形式，

$M_2=\left[\begin{array}{cc}I& H_{12}{H}_{22}^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{H}_{11}& 0\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}I& 0\\ H_{22}^{-1}{H}_{21}& I\end{array}\right]\text{.}$ (19)

预处理矩阵 $M_2$ 具有以下谱性质。

定理2 预处理矩阵 $H{M}_2^{-1}$ 至少有 $2n^2$ 个单位特征值和 $2n^2$ 个单位特征值对应的线性无关的特征向量。 $H{M}_2^{-1}$ 其他的特征值分布在 $\left[1-\delta ,1\right]$ ，其中 $\delta =\max {\left(1+{\psi }^{″}\left(v_k\right)\right)}^{-1}$ 。

Wang等 [28] 基于Hessian矩阵H的分解形式，对H₁₁的Schur补S₁的逆矩阵进行近似，构建出第一个预处理矩阵P₁，其逆矩阵形式如下，

$P_1^{-1}=\left[\begin{array}{cc}I& -H_{11}^{-1}{H}_{12}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{H}_{11}^{-1}& 0\\ 0& {\hat{S}}_1^{-1}\end{array}\right]\left[\begin{array}{cc}I& 0\\ -H_{21}{H}_{11}^{-1}& I\end{array}\right].$ (20)

其中 $S_1=H_{22}^{\frac{1}{2}}\left(1-T_1\right)H_{22}^{\frac{1}{2}}$ ， $T_1=H_{22}^{\frac{1}{2}}{H}_{21}{H}_{11}^{-1}{H}_{12}{H}_{22}^{\frac{1}{2}}$ ，因H为对称正定的，显然S₁为对称正定矩阵。而近似Schur补的逆通过取其泰勒展开式的线性部分近似 ${\hat{S}}_1^{-1}=H_{22}^{-1}+H_{22}^{-1}{H}_{21}{H}_{11}^{-1}{H}_{12}{H}_{22}^{-1}$ 。

预处理矩阵P₁具有以下谱性质。

定理3 预处理矩阵 $H{P}_1^{-1}$ 至少有 $n^2$ 个单位特征值和 $n^2$ 个单位特征值对应的线性无关的特征向量。 $H{P}_1^{-1}$ 其他的特征值分布在 $\left[1-{\delta }^2,1\right]$ ，其中 $\delta =\max {\left(1+{\psi }^{″}\left(v_k\right)\right)}^{-1}$ 。

类似地，对H₂₂的Schur补S₂的逆矩阵进行近似，得到第二个预处理矩阵P₂其逆矩阵形式如下，

$P_2^{-1}=\left[\begin{array}{cc}I& 0\\ -H_{22}^{-1}{H}_{21}& I\end{array}\right]\left[\begin{array}{cc}{\hat{S}}_2^{-1}& 0\\ 0& H_{22}^{-1}\end{array}\right]\left[\begin{array}{cc}I& -H_{12}{H}_{22}^{-1}\\ 0& I\end{array}\right].$ (21)

其中 $S_2=H_{11}^{\frac{1}{2}}\left(1-T_2\right)H_{11}^{\frac{1}{2}}$ ， $T_2=H_{11}^{\frac{1}{2}}{H}_{12}{H}_{22}^{-1}{H}_{21}{H}_{11}^{\frac{1}{2}}$ ，因H为对称正定的，显然S₂为对称正定矩阵。而近似Schur补的逆通过取其泰勒展开式的线性部分实现 ${\hat{S}}_2^{-1}=H_{11}^{-1}+H_{11}^{-1}{H}_{12}{H}_{22}^{-1}{H}_{21}{H}_{11}^{-1}$ 。

预处理矩阵P₂具有以下谱性质。

定理4 预处理矩阵 $H{P}_2^{-1}$ 至少有 $2n^2$ 个单位特征值和 $2n^2$ 个单位特征值对应的线性无关的特征向量。 $H{P}_2^{-1}$ 其他的特征值分布在 $\left[1-{\delta }^2,1\right]$ ，其中 $\delta =\max {\left(1+{\psi }^{″}\left(v_k\right)\right)}^{-1}$ 。

Zhao和Huang [29] 构造了一个限制性预处理矩阵Z如下，

$Z=\left[\begin{array}{cc}I& 0\\ H_{21}{H}_{11}^{-1}& I\end{array}\right]\left[\begin{array}{cc}{H}_{11}& 0\\ 0& {\tilde{S}}_1\end{array}\right]\left[\begin{array}{cc}I& H_{11}^{-1}{H}_{12}\\ 0& I\end{array}\right]\text{.}$ (22)

其中 ${\tilde{S}}_1=H_{22}-dI$ ，选定合适的常数d使得dI近似D，从而 ${\tilde{S}}_1$ 近似替代了Hessian矩阵关于H₁₁的Schur补形式 $S_1=H_{22}-H_{21}{H}_{11}^{-1}{H}_{12}=H_{22}-UD{U}^{\text{T}}$ ，

$D={\beta }^2{\Sigma }^2{\left(2{\Lambda }^2+\beta \tau {\Sigma }^2\right)}^{-1}$

是一个对角矩阵，第i个对角项是 $d_i$ 如下，

$d_i=\{\begin{array}{l}\frac{{\beta }^2{\sigma }^2}{2{\lambda }_i^2+\beta \tau {\sigma }_i^2},i=1,2,\cdots ,r\\ 0,i=r+1,\cdots ,mn.\end{array}$

其中 $A=F\Lambda F^{\text{T}}$ ， $G^{\text{T}}G=F{\Sigma }^2{F}^{\text{T}}$ ，F是快速傅里叶变换矩阵， $\Lambda =\text{diag}\left({\lambda }_1,\cdots ,{\lambda }_{mn}\right)$ ， $\Sigma =\text{diag}\left({\sigma }_1,\cdots ,{\sigma }_r,0,\cdots ,0\right)$ ，r是矩阵 $G^{\text{T}}G$ 的秩。

预处理矩阵Z具有以下谱性质。

定理5 对于到 $\psi \left(v\right)$ 和D，有

$v_{\min }=\underset{i}{\min }\left\{{\psi }^{″}\left(v_i\right)\right\},d_{\min }=\underset{i}{\min }\left\{d_i\right\},d_{\max }=\underset{i}{\max }\left\{d_i\right\},$

其中 $d_i\left(i=1,\cdots ,mn\right)$ 是矩阵D的对角线元素。假设 $d_{\max }<\beta /\tau$ ， ${\tilde{S}}_1$ 中的d满足 $d_{\min }<d<d_{\max }$ 。矩阵 $H^{-1}$ 至少有mn个单位特征值，且其他特征值分布于区间

$\left[1-\frac{d_{\max }-d}{\left(\beta /\tau -d\right)+\beta v_{\min }},1+\frac{d-d_{\min }}{\left(\beta /\tau -d\right)+\beta v_{\min }}\right]\text{.}$

3.1.2. 对比分析

方程组(11) Hessian矩阵的完全分解形式为，

$H=\left[\begin{array}{cc}I& 0\\ H_{21}{H}_{11}^{-1}& I\end{array}\right]\left[\begin{array}{cc}{H}_{11}& 0\\ 0& S_1\end{array}\right]\left[\begin{array}{cc}I& H_{11}^{-1}{H}_{12}\\ 0& I\end{array}\right]\text{.}$ (23)

其中 $S_1=H_{22}-H_{21}{H}_{11}^{-1}{H}_{12}$ 为H₁₁的Schur补。而另一种完全分解形式为，

$H=\left[\begin{array}{cc}I& H_{12}{H}_{22}^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{S}_2& 0\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}I& 0\\ H_{22}^{-1}{H}_{21}& I\end{array}\right]\text{.}$ (24)

其中 $S_2=H_{11}-H_{12}{H}_{22}^{-1}{H}_{21}$ 为H₂₂的Schur补形式。

通过对比分析，我们发现，上述几种预处理方法本质上都是基于Schur补近似来构造的预处理矩阵，只是近似的策略不同。从Schur补近似形式的角度看，Huang是将S₁近似为H₂₂，将S₂近似为H₁₁，Zhao是将S₁近似为 $H_{22}-dI$ 。不同的是，Wang是基于泰勒展开对S₁和S₂的逆矩阵进行近似。从近似程度的角度看，Wang提出的预处理与Zhao提出的预处理比之Huang的预处理，更能近似系数矩阵H。上述预处理对应的谱分析和数值结果也印证了这一点，后两种预处理系数矩阵的特征值更加聚集，应用到PCG方法时所需要的迭代步数更少。

Huang论文中的近似是构造最简单计算量也最小的方法，但同时意味着其近似程度有着更高的提升空间，可以基于Shcur补进一步改进其对S₁与S₂的近似程度。另一种想法是在两种矩阵分解中，不对S₁与S₂做近似，转而对H₁₁或H₂₂做其他的近似变换，由于H₁与H₂组成不同，其各自近似效果也会有差距。而Wang论文中通过泰勒展开到二次项再取线性部分近似，对此，一种直接的改进想法就是泰勒展开到更高次项，使得矩阵的近似程度更高，但这将会引起每步更多计算量的问题。Zhao的论文中只是把对角阵近似为数量矩阵，这个近似比较粗糙，还涉及到参数的选取问题，将来可以研究最优参数的选择，或者找到更优的近似方式，来改进相应的预处理方法。当然，不论是怎样的改进构造想法，如何控制好近似程度与计算量之间的平衡是较为关键的。

3.2. 乘性形式

3.2.1. 预处理形式及性质

对于乘性形式，Huang和Lu [19] 提出了预处理矩阵 $M_0$ 如下，

$M_0=\left[\begin{array}{cc}{\tilde{H}}_{11}+H_{12}{H}_{22}^{-1}{H}_{21}& H_{12}\\ H_{21}& H_{22}\end{array}\right],$ (25)

其中， ${\tilde{H}}_{11}$ 是 $H_{11}$ 的近似，通过用单位矩阵I来近似代替 $\Lambda \left(v\right)$ 得到 ${\tilde{H}}_{11}=2A^{\text{T}}A+2\beta G^{\text{T}}G$ ，这种近似是考虑到 $2A^{\text{T}}A+2\beta G^{\text{T}}G$ 在周期边界条件下是块状循环矩阵，相关的运算可以通过快速傅里叶变换(FFT)实现。 $M_0$ 具有以下分解形式，

$M_0=\left[\begin{array}{cc}I& H_{12}{H}_{22}^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{\tilde{H}}_{11}& 0\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}I& 0\\ H_{22}^{-1}{H}_{21}& I\end{array}\right]\text{.}$ (26)

预处理矩阵 $M_0$ 具有以下谱性质。

定理6 对于对角矩阵Λ对角线元素的最大值与最小值 ${\sigma }_{\max }$ 与 ${\sigma }_{\min }$ ，矩阵 $H{M}_0^{-1}$ 至少有 $2n^2$ 个单位特征值，且其他特征值都分布在区间 $\left[\min \left({\sigma }_{\min },1\right),\max \left({\sigma }_{\max },1\right)\right]$ 内。矩阵 $H{M}_0^{-1}$ 对应于单位特征值有 $2n^2$ 个

的线性无关的特征向量，且对应于非单位特征值 $\lambda$ 的其他特征向量都可以由 $\left[\begin{array}{c}u\\ 0\end{array}\right]$ 得到，其中u是矩阵

$\left(H_{11}-H_{12}{H}_{22}^{-1}{H}_{21}\right){\tilde{H}}_{11}$ 关于 $\lambda$ 的特征向量。

然而，使用单位矩阵I对H₁₁的近似较为简略，作者进一步提出 ${\hat{H}}_{11}=2A^{\text{T}}A+2\theta \beta G^{\text{T}}G$ ( $\theta >0$ 是一个常数)，其中用 $\theta I$ 改进了对 $\Lambda \left(v\right)$ 的近似代替，此预处理矩阵 $M_{\theta }$ 如下，

$M_{\theta }=\left[\begin{array}{cc}{\hat{H}}_{11}+H_{12}{H}_{22}^{-1}{H}_{21}& H_{12}\\ H_{21}& H_{22}\end{array}\right],$ (27)

具有以下分解形式，

$M_{\theta }=\left[\begin{array}{cc}I& H_{12}{H}_{22}^{-1}\\ 0& I\end{array}\right]\left[\begin{array}{cc}{\hat{H}}_{11}& 0\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}I& 0\\ H_{22}^{-1}{H}_{21}& I\end{array}\right]\text{.}$ (28)

预处理矩阵 $M_{\theta }$ 具有以下谱性质。

定理7 对于对角矩阵Λ对角线元素的最大值与最小值 ${\sigma }_{\max }$ 与 ${\sigma }_{\min }$ ，矩阵 $H{M}_{\theta }^{-1}$ 至少有 $2n^2$ 个单位特征值，且其他特征值都分布在区间 $\left[\min \left(\frac{{\sigma }_{\min }}{\theta },\frac{1}{\theta }\right),\max \left(\frac{{\sigma }_{\max }}{\theta },\frac{1}{\theta }\right)\right]$ 内。矩阵 $H{M}_{\theta }^{-1}$ 对应于单位特征值有 $2n^2$ 个的线性无关的特征向量，且对应于非单位特征值 $\lambda$ 的其他特征向量都可以由 $\left[\begin{array}{c}u\\ 0\end{array}\right]$ 得到，其中u是矩阵 $\left(H_{11}-H_{12}{H}_{22}^{-1}{H}_{21}\right){\hat{H}}_{11}$ 关于 $\lambda$ 的特征向量。

前面所述的预处理都是基于Hessian矩阵的不完全分解对矩阵进行处理，而Zhao和Huang [30] 基于超松弛迭代(SSOR)将矩阵分解成三部分，分别是两个上下三角形矩阵和一个对角阵，

$H=L+D+L^{\text{T}}=\left[\begin{array}{cc}0& 0\\ H_{21}& 0\end{array}\right]+\left[\begin{array}{cc}{H}_{11}& 0\\ 0& H_{22}\end{array}\right]+\left[\begin{array}{cc}0& H_{12}\\ 0& 0\end{array}\right]$ (29)

对应的SSOR迭代方程组如下，

$\{\begin{array}{l}\left(D-\omega L\right)d^{k+\frac{1}{2}}=\left(\omega L^{\text{T}}+\left(1-\omega \right)D\right)d^k+\omega r\\ \left(D-\omega L^{\text{T}}\right)d^{k+1}=\left(\omega L+\left(1-\omega \right)D\right)d^{k+\frac{1}{2}}+\omega r\end{array}$ (30)

由于计算成本太大，提出了改进的SSOR迭代法，基于其中的矩阵分裂，通过定义非奇异矩阵 $Q\left(\omega \right)$ ，构造了改进的SSOR块预处理矩阵P如下，

$P\left(\omega \right)={\left(\tilde{D}+\omega L\right)}^{\text{T}}Q{\left(\omega \right)}^{-1}\left(\tilde{D}+\omega L\right),$ (31)

$Q\left(\omega \right)=\omega \left(\tilde{D}+{\tilde{D}}^{\text{T}}-\omega D\right),$ (32)

其中，非奇异矩阵 $\tilde{D}$ 是对分块对角矩阵D的近似替代， $\tilde{D}=\left[\begin{array}{cc}{\tilde{H}}_{11}& 0\\ 0& H_{22}\end{array}\right]$ ，通过简单地将 $\Lambda \left(v\right)$ 替换为标量矩阵 $\theta I$ 对系数矩阵中的H₁₁块进行修正， ${\tilde{H}}_{11}=2A^{\text{T}}A+2\theta \beta G^{\text{T}}G$ 。

当对图像施加周期性边界条件时，模糊矩阵A是带循环块的循环矩阵(BCCB)，如果G是一阶差分地离散化矩阵，A和 $G^{\text{T}}G$ 可以通过快速傅里叶变换(FFT)实现对角化，且上述定义的预处理矩阵P就变成，

$P\left(\omega ,\theta \right)=\left[\begin{array}{cc}{\tilde{H}}_{11}& \omega H_{12}\\ 0& H_{22}\end{array}\right]\left[\begin{array}{cc}\frac{1}{\omega }{\left(2{\tilde{H}}_{11}-\omega H_{11}\right)}^{-1}& 0\\ 0& \frac{1}{\omega \left(2-\omega \right)}{H}_{22}^{-1}\end{array}\right]\left[\begin{array}{cc}{\tilde{H}}_{11}& 0\\ \omega H_{21}& H_{22}\end{array}\right]$ (33)

作者对条件数进行了详细证明与分析，并给出了条件数上界与下界的具体参考 [30] 。

3.2.2. 对比分析

从近似形式的角度看，我们发现Huang论文中的预处理矩阵是基于(24)矩阵分解式，而Zhao论文中的预处理矩阵是基于(29)的矩阵分裂，前者结合Shur补进行近似，后者从定常迭代法中的SSOR迭代中沿用了矩阵分裂思想进行近似。从近似程度的角度看，Zhao在近似中引入了SSOR法中的超松弛因子 $\omega$ ，得以有效地控制近似程度，而Huang是将对角矩阵近似为数量矩阵，有效减少了计算量，但近似程度却不如Zhao提出的预处理矩阵。

根据以上对比分析，我们的一种想法是可以改进Huang中对于S₁与S₂的简单近似，或者对分解矩阵中的例如H₁₁与H₂₂做其他近似变换，亦或类似加性形式预处理中Zhao所采用的策略来近似对角元。另一种想法是将形式较为复杂的(33)作进一步简化，保留超松弛因子的灵活性，同时能够降低计算量且使得谱性质更易于分析。此外，针对半二次图像复原中方程组系数矩阵的结构化特点，可以结合更多的矩阵分裂或分解模式加以改进，或是结合更精确的边界条件进行近似构造。

4. 总结

目前，针对半二次图像复原中的结构化方程组，学者们基于Schur补近似等策略，构造了一系列的预处理方法，有效地降低了图像复原中求解方程组的计算成本。本文整理了近些年出现的预处理方法，从不同角度进行了对比分析。我们发现，对于加性形式下的结构化方程组，预处理矩阵的构造主要是基于Hessian矩阵的不完全分解和Schur近似。不同的预处理方法采用了不同的Schur近似策略。典型的有，一是直接将Schur补近似为系数矩阵的子矩阵，二是通过对角矩阵的近似和BCCB矩阵的傅里叶分解来近似Schur补，三是基于逆矩阵的分解和截断的泰勒展开来实现Schur补的近似。对于乘性形式下的结构化方程组，还可以基于SSOR等定常迭代法的矩阵分裂来完成对Hessian矩阵的近似，进而构造出含有参数的预处理矩阵。结果表明，相比原有的方程组，以上预处理后的系数矩阵的特征值更加聚集，PCG的迭代步数明显减少。

但是，目前已有的预处理方法仍然存在一定的局限性，值得进一步改进和完善。关于半二次图像复原问题中结构化方程组预处理方法的研究，本文认为有如下几个方面仍值得继续探索：

一，构造更有效精确的近似形式。在改进预处理方向上，可以尝试结合边界条件或泰勒展开到更高次项等进行更精确的近似，也可以尝试在矩阵分解基础上对其他关键部分进行近似，其中最为重要的考量是平衡好近似程度与计算成本的代价。

二，构造参数自适应的预处理方法。通过对角矩阵的近似，或者基于SSOR迭代法中的矩阵分裂等方法来构造预处理矩阵将不可避免地引入超松弛因子等参数，给预处理的谱分析和数值实验带来困难和不确定性。在预处理构造中，如何根据模型信息来自适应地调整每步优化迭代中的最佳参数，仍需要进一步研究。

基金项目

国家自然科学基金项目(项目编号：12001022)。

国家级大学生创新创业训练计划项目(项目编号：202310011021)。

NOTES

^*通讯作者。