离散不适定问题的扩展迭代正则化方法

The Enriched Iterative Regularization Method for Discrete Ill-Posed Problems

摘要: Arnoldi-Tikhonov方法是求解大规模离散不适定问题的一种常用子空间迭代正则化方法，其由Arnoldi算法构建低维Krylov子空间，再对低维问题用Tikhonov正则化从而获得正则化解。但由于低维子空间信息缺失，正则化解的有时效果欠佳。为了改进正则化效果，本文通过增加一个含有特定先验信息的低维子空间来扩展Krylov子空间，提出了求解大规模离散不适定问题的一种扩展子空间迭代正则化方法。该方法通过扩展Arnoldi算法构建扩展子空间，并融合Tikhonov正则化，从而获得更优正则化解。针对经典算例，将所提算法与Arnoldi-Tikhonov算法进行了数值实验和性态比较，数值结果验证了所提算法的有效性。

Abstract: Arnoldi-Tikhonov method is one of the often used Krylov subspace iterative regularization methods for solving large scale discrete ill-posed problems. It generates the Krylov subspace through Arnoldi process and gets the regularized solution through applying Tikhonov regularization method to the projected small problem. However, due to the information deficiency of the dimension reduced subspace, the regularized solution by Arnoldi-Tikhonv method sometimes is not as good as expected. In order to improve the regularized solution, an enriched subspace iterative regularization method is proposed in this paper. The proposed new method enriches the Krylov subspace by adding a special subspace that holds some specific prior information. Numerical experiments are carried out. The numerical results show that the proposed method is more effective than the Arnoldi-Tikhonov method.

1. 引言

考虑求解如下大规模离散不适定问题：

$\underset{x}{\min }{\Vert Ax-b\Vert }_2^{2}A\in R^{m\times n}$ . (1)

此类问题经常出现于科学和工程等诸多领域[1]，通常是由连续不适定问题的离散化产生[2]，例如，具有紧致核的第一类Fredholm积分方程的离散化等。求解离散不适定问题的常用方法是Tikhonov正则化方法[3] [4]，通常采用以下形式：

$\underset{x}{\min }\left\{{\Vert Ax-b\Vert }_2^2+{\lambda }^2{\Vert x\Vert }_2^2\right\}$ (2)

其中$\lambda$ 是正则化参数[5] [6]，意图控制着${\Vert Ax-b\Vert }_2^2$ 和${\Vert x\Vert }_2^2$ 组合权重。

对于大型离散不适定问题，由于矩阵A的规模过大，限于存储量和计算量，应用Tikhonov正则化方法直接求解该问题正则化解是不现实的。

因此，需要将矩阵A投影在一个低维空间中，得到一个低维空间中的不适定问题，即：

$\underset{x}{\min }{\Vert Ax-b\Vert }_2^2\text{s}\text{.t}\text{.}x\in S_k$

这里的子空间S_k是一个维数为k的线性子空间。如果能够得到S_k的一组合适的基$\text{span}\left\{v_1,v_2,v_3,\cdots ,v_k\right\}$ ，那么这种方法可以取得很好的效果。于是，通常选取S_k是Krylov子空间[4]：

$K_k=\text{span}\left\{b,Ab,A^{2}b,\cdots ,A^{k-1}b\right\}$ .

例如，GMRES、RRGMRES方法[7]-[12]和LSQR方法[13]都是使用到了此类子空间。

对于不同的问题，上述子空间中的一个或多个可能非常适合用来计算出一个好的正则化解，也就是说，这个正则化解可以对精确解有着很好的估计，并且对于数据的扰动不敏感。将低维Krylov子空间投影与Tikhonov正则化相结合，即将Tikhonov正则化方法应用在原问题的投影问题，这就是所谓的经典子空间迭代正则化算法。特别地，如果使用Arnoldi算法投影，即为Arnoldi-Tikhonov正则化算法。

然而，对于一些矩阵A和向量b，构造Krylov子空间可能不包含可以表示精确解$x_{exact}$ 已知特征的向量。由于Arnoldi-Tikhonov正则化算法中生成的Krylov子空间的维数通常远小于矩阵A的维数，所以通过该Krylov子空间不能准确地表示这样的解。在这种情况下，使用Arnoldi-Tikhonov正则化算法就会得不到理想的效果，即得到的解不是精确解$x_{exact}$ 好的近似。

于是，本文对利用Arnoldi算法进行Krylov子空间进行了扩展，通过加入特定的先验信息，来使得到的正则化解更加准确。将这种新方法称为扩展Arnoldi算法。该方法的主要思想是：假设所需求得的解在子空间W_p中有一个重要的组成部分，这里的W_p是一个给定的子空间，它的维数为$p<k$ ，将该子空间与上述Krylov子空间相结合。于是。可以将所需求得的解分称两个部分，其中一个部分在给出的子空间W_p中，另一部分在W_p的正交补空间$W_p^{\perp }$ 中。具体来说，在一个新的子空间$S_{p,k}$ 中计算问题的正则化解，这里的子空间$S_{p,k}$ 是上述两个子空间W_p和K_k的直和

$S_{p,k}=W_p+K_k\equiv \left\{x+y|x\in W_p\wedge y\in K_k\right\}$ ,

这里的$S_{p,k}$ 是一个线性子空间。

需要强调的一点是所选择的子空间必须是合适的，只有一个好的扩展子空间W_p才能使得到的解比原方法得到的正则化解更加精确。如果所选的W_p不够理想，那么通过该方法得到的解可能与精确解相差甚远，因此子空间W_p的选取十分重要。

基于上述思路，对Arnoldi-Tikhonov正则化算法进行扩展，称算法称为扩展Arnoldi迭代正则化算法。与传统Arnoldi-Tikhonov正则化算法类似，可以计算出扩展Arnoldi迭代正则化算法相应的Tikhonov正则化参数，并且在子空间W_p选择合适的情况下，通过该方法得到的正则化解通常优于传统Arnoldi-Tikhonov正则化算法。

本文的组织结构如下：第2节对Arnoldi算法进行扩展，得到了扩展Arnoldi算法。在第3节中，融合扩展Arnoldi算法和Tikhonov正则化，给出了扩展Arnoldi迭代正则化算法，给出了具体的算法描述。第4节进行了一系列数值实验，来展示和验证本文所提出的扩展Arnoldi迭代正则化算法的有效性和优越性。

2. 扩展Arnoldi算法

本节讨论如何实现扩展子空间算法。类似于Arnoldi算法构造Krylov子空间的标准正交基，实施推广的Hessenberg分解

$A\left[W_p,b,Ab,A^{2}b,\cdots ,A^{k-1}b\right]=V_{p+k+1}{H}_{p+k}$ .

于是，可以计算出问题的解：

$x^{\left(k\right)}=\left[W_p,b,Ab,A^{2}b,\cdots ,A^{k-1}b\right]y^{\left(k\right)}$

$y^{\left(k\right)}=\arg {\min }_y{\Vert H_{p+1}y-V_{p+k+1}^{\text{T}}b\Vert }_2^2$

对于Krylov子空间K_k的构造，可以通过Arnoldi算法来计算出其一组标准正交基；接下来将子空间W_p扩展到这组标准正交基中，使之成为一个新的扩展的子空间。当Arnoldi算法的迭代次数为k时，可以得到如下形式的矩阵分解：

$A\left[V_k,W_p\right]=\left[V_{k+1},{\tilde{V}}_k\right]\left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\end{array}\right]$ .

在上面的公式中V_k，$V_{k+1}$ ，和B_k都是与经典的Arnoldi算法相关的，其余的矩阵与加入的扩展子空间W_p有关。

具体来说，

• $A{V}_k=V_{k+1}{B}_k$ 是由Arnoldi过程经过k次迭代后得到的；

• $V_k\in {ℝ}^{n\times k}$ 有着标准正交的列向量，并且可以张成为Krylov子空间K_k；

• $V_{k+1}\in {ℝ}^{n\times \left(k+1\right)}$ 有着标准正交的列向量，并且它是由矩阵V_k扩展了一列之后形成的；

• $B_k\in {ℝ}^{\left(k+1\right)\times k}$ 是一个上Hessenberg矩阵；

• ${\tilde{V}}_k\in {ℝ}^{n\times p}$ 满足等式$A{W}_p=V_{k+1}{G}_k+{\tilde{V}}_k{F}_k$ 并且${\tilde{V}}_k^{\text{T}}{V}_{k+1}=0$；

• $F_k\in {ℝ}^{p\times p}$ 与Arnoldi过程的迭代次数有关，迭代的次数不同，F_k也会不同；

• $G_k\in {ℝ}^{\left(k+1\right)\times p}$ 与矩阵$V_{k+1}$ 有关。

由上述矩阵分解可知$\left[V_{k+1},{\tilde{V}}_k\right]$ 的列向量形成了子空间$S_{p,k}$ 的一组基。矩阵G_k和F_k分别是由$V_{k+1}$ 和${\tilde{V}}_k$ 来计算得到的，计算公式如下所示：

$G_k=V_{k+1}^{\text{T}}A{W}_p\text{,}{F}_k\in {\tilde{V}}_k^{\text{T}}A{W}_p$.

于是，原问题可转换为求解如下最小二乘问题：

$y^{\left(k\right)}=\arg {\min }_y{\Vert \left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\end{array}\right]y-\left[\begin{array}{c}{V}_{k+1}^{\text{T}}\\ {\tilde{V}}_k^{\text{T}}\end{array}\right]b\Vert }_2^2$

通过解决此最小二乘问题求得$y^{\left(k\right)}$ ，进而可以求得$x^{\left(k\right)}$ ：

$x^{\left(k\right)}=\left[V_k,W_p\right]y^{\left(k\right)}$ .

于是，给出具体的扩展算法如下算法1：

Algorithm 1. Enriched Arnoldi algorithm

算法1. 扩展Arnoldi算法

1) 设定$V_1=b/{\Vert b\Vert }_2$ ，$V_0=\left[\right]$ ，$B_0=\left[\right]$ ，$G_0=V_1^{\text{T}}A{W}_p$ ，$k=1$ ；

2) 通过Arnoldi过程来计算得出$v_{k+1}$ 使得等式$A{V}_k=V_{k+1}{B}_k$ 成立，其中$V_{k+1}=\left[V_k,v_{k+1}\right]$ ，$B_k=\left[\begin{array}{cc}{B}_{k-1}& \begin{array}{c}0\\ \times \end{array}\\ 0& \times \end{array}\right]$ ；

3) 计算$G_k=\left[\begin{array}{c}{G}_{k-1}\\ v_{k+1}^{\text{T}}A{W}_p\end{array}\right]\in {ℝ}^{\left(k+1\right)\times p}$ ；

4) 通过$V_{k+1}$ 对$A{W}_p$ 进行标准正交化，从而计算得到${\tilde{V}}_k\in {ℝ}^{n\times p}$ ；

5) 计算$F_k={\tilde{V}}_k^{\text{T}}A{W}_p\in {ℝ}^{p\times p}$；

6) 通过解决最小二乘问题${\min }_y{\Vert \left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\end{array}\right]y-\left[\begin{array}{c}{V}_{k+1}^{\text{T}}\\ {\tilde{V}}_k^{\text{T}}\end{array}\right]b\Vert }_2^2$来得到解$y^{\left(k\right)}$ ；

7) 得到$x^{\left(k\right)}=\left[V_k,W_p\right]y^{\left(k\right)}$ ；

8) 在迭代次数达到要求时停止，否则令$k=k+1$ 并重复步骤2。

对比原算法，上述扩展Arnoldi算法需要对矩阵W_p，${\tilde{V}}_k$ ，G_k和F_k进行额外的计算，但是由于选择的扩展子空间W_p维数很小，计算它们的开销可以忽略不计。

3. 扩展Arnoldi迭代正则化算法

本节融合扩展Arnoldi算法和Tikhonov正则化，给出扩展Arnoldi正则化算法。

基于扩展Arnoldi算法，得到

$A\left[V_k,W_p\right]=\left[V_{k+1},{\tilde{V}}_k\right]\left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\end{array}\right]$ .

将此分解带入上述问题(2)，进行运算处理后得到如下低维问题：

$y_{{\lambda }_k}^{\left(k\right)}=\arg {\min }_y{\Vert \left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\\ {\lambda }_k{V}_k& {\lambda }_k{W}_p\end{array}\right]y-\left[\begin{array}{c}{V}_{k+1}^{\text{T}}\\ {\tilde{V}}_k^{\text{T}}\\ 0\end{array}\right]b\Vert }_2^2$ (3)

其中$y_{{\lambda }_k}^{\left(k\right)}$ 满足如下的等式：

$x_{{\lambda }_k}^{\left(k\right)}=\left[V_k,W_p\right]y_{{\lambda }_k}^{\left(k\right)}$ .

于是对于低维子空间最小二乘问题运用Tikhonv正则化算法，就实现了扩展子空间投影正则化过程。这种算法的优势在于：一是通过扩展Arnoldi算法获得具有特定信息的投影低维子空间，实现投影降维和信息补充的双重效果；二是可以更为灵活地选用正则化参数$\lambda$ ，包括偏差原理[14]，L-曲线准则[15]-[17]以及广义交叉验证方法[18]等。本文选用广义交叉验证方法来作为正则化参数的选取方法[19]。

于是，给出扩展Arnoldi正则化算法(算法2)描述：

Algorithm 2. Enriched Arnoldi regularization algorithm

算法2. 扩展Arnoldi正则化算法

1) 设定$V_1=b/{\Vert b\Vert }_2$ ，$V_0=\left[\right]$ ，$B_0=\left[\right]$ ，$G_0=V_1^{\text{T}}A{W}_p$ ，$k=1$ ；

2) 通过Arnoldi过程来计算得出$v_{k+1}$ 使得等式$A{V}_k=V_{k+1}{B}_k$ 成立，其中$V_{k+1}=\left[V_k,v_{k+1}\right]$ ，$B_k=\left[\begin{array}{cc}{B}_{k-1}& \begin{array}{c}0\\ \times \end{array}\\ 0& \times \end{array}\right]$ ；

3) 计算$G_k=\left[\begin{array}{c}{G}_{k-1}\\ v_{k+1}^{\text{T}}A{W}_p\end{array}\right]\in R^{\left(k+1\right)\times p}$ ；

4) 通过$V_{k+1}$ 对$A{W}_p$ 进行标准正交化，计算${\tilde{V}}_k\in {ℝ}^{n\times p}$ ；

5) 计算$F_k={\tilde{V}}_k^{\text{T}}A{W}_p\in R^{p\times p}$；

6) 使用广义交叉验证方法获得正则化参数$\lambda$ ，计算${\min }_y{\Vert \left[\begin{array}{cc}{B}_k& G_k\\ 0& F_k\\ {\lambda }_k{V}_k& {\lambda }_k{W}_p\end{array}\right]y-\left[\begin{array}{c}{V}_{k+1}^{\text{T}}\\ {\tilde{V}}_k^{\text{T}}\\ 0\end{array}\right]b\Vert }_2^2$的正则化解$y_{{\lambda }_k}^{\left(k\right)}$ ；

7) 于是得到原始正则化解$x_{{\lambda }_k}^{\left(k\right)}=\left[V_k,W_p\right]y_{{\lambda }_k}^{\left(k\right)}$ ；

8) 在迭代次数达到要求时停止，否则令$k=k+1$ 并重复步骤2。

对于大规模离散不适定问题，通过扩展Arnoldi正则化算法计算得到正则化解通常优于通过经典Arnoldi-Tikhonov正则化算法得到的解。在下一节中，通过一系列数值算例来验证其有效性和优越性。

4. 数值实验

本节将给出一些数值例子验证扩展子空间迭代正则化算法的效能。

具体实验参数说明如下：考虑含有噪音来自第一类Fredholm积分方程离散化得到的离散不适定问题

$A{x}_{exact}=b_{exact}\text{+}e$

这里A分别取deriv2、foxgood和baart矩阵问题[2]；e为随机高斯白噪声向量，其噪声水平为$\eta ={\Vert e\Vert }_2/{\Vert b_{exact}\Vert }_2$ ；相对误差的计算方法为${\Vert x_{exact}-x^{\left(k\right)}\Vert }_2/{\Vert x_{exact}\Vert }_2$ ，即用计算所得正则化解与精确解在2-范数意义下的相对误差来描述该解的准确程度。

分别比较扩展Arnoldi迭代正则化方法(Enriched Arnoldi-Tikhonov，简记为E-A-T算法)和Arnoldi-Tikhonov正则化方法(Arnoldi-Tikhonov，简记为A-T算法)所得解的准确度。在数值实验中，统一使用广义交叉验证方法(GCV)作为正则化参数的选取方法。

4.1. 例1

考虑第一个问题是deriv2 (n, 2)问题，选择扩展子空间为

$W_2=\text{span}\left\{w_1,w_2\right\},w_1={\left(1,1,\cdots ,1\right)}^{\text{T}},w_2={\left(1,2,\cdots ,n\right)}^{\text{T}}$

设置n的值为32，相对噪音水平$\eta ={10}^{-6}$ ，选择Arnoldi过程的迭代次数k为20。通过上述两种方法对该问题进行了求解，结果如图1、图2所示。

表1也展示了扩展Arnoldi迭代正则化方法可以明显得提升计算得出的解的准确性。

Table 1. Comparing the relative errors of the two algorithms of the first case in example 1

表1. 例1的第一种情况中两种算法相对误差数值对比

算法	A-T算法	E-A-T算法
相对误差	5.2780e−02	3.8089e−04

Figure 1. The regular solution obtained by A-T algorithm of the first case in example 1

图1. 例1的第一种情况中通过A-T算法计算得到的正则化解

Figure 2. The regular solution obtained by E-A-T algorithm of the first case in example 1

图2. 例1的第一种情况中通过E-A-T算法计算得到的正则化解

同样对于deriv2 (n, 2)问题，设置n的值为1000，并且相对噪音水平$\eta ={10}^{-3}$ ，选择k为50。类似的，也得到了该问题相应的解曲线和相对误差，分别如图3、图4以及表2所示。

Table 2. Comparing the relative errors of the two algorithms of the second case in example 1

表2. 例1的第二种情况中两种算法相对误差数值对比

算法	A-T算法	E-A-T算法
相对误差	1.6747e−01	1.7762e−02

Figure 3. The regular solution obtained by A-T algorithm of the second case in example 1

图3. 例子的第二种情况中通过A-T算法计算得到的正则化解

Figure 4. The regular solution obtained by E-A-T algorithm of the second case in example 1

图4. 例1的第二种情况中通过E-A-T算法计算得到的正则化解

可以看到通过扩展Arnoldi迭代正则化方法得到的解表现得更好，因为它的解曲线更加光滑，并且与精确解曲线更加接近。可见，对于此类问题，本文所提算法是有效且有优势的。

4.2. 例2

考虑foxgood (n)问题，选择扩展子空间与第一个例子相同，如下所示

$W_2=\text{span}\left\{w_1,w_2\right\},w_1={\left(1,1,\cdots ,1\right)}^{\text{T}},w_2={\left(1,2,\cdots ,n\right)}^{\text{T}}$

Figure 5. The regular solution obtained by A-T algorithm in example 2

图5. 例2中通过A-T算法计算得到的正则化解

Figure 6. The regular solution obtained by E-A-T algorithm in example 2

图6. 例2中通过E-A-T算法计算得到的正则化解

Table 3. Comparing the relative errors of the two algorithms in example 2

表3. 例2中两种算法相对误差数值对比

算法	A-T算法	E-A-T算法
相对误差	1.2207e−02	4.4613e−04

设置n的值为1000，设置相对噪音水平$\eta ={10}^{-3}$ ，选择迭代次数k的值为30。同样的，通过上述两种方法对该问题进行了求解，并且将解向量的图像与相对误差的结果展示在图5、图6以及表3中。从图5中可以看出，通过Arnoldi-Tikhonov正则化方法得到的解向量的前100项与精确解向量的前100项差距较大，而从图6来看通过本文所提方法得到的正则解向量很好的修正了这一点，使得解向量整体上都与精确解向量十分相近。正因如此，通过本文的方法求得的解相对误差比Arnoldi-Tikhonov正则化方法要小得多，如表3所示。

可见，对于该问题，扩展Arnoldi迭代正则化方法依然比Arnoldi-Tikhonov正则化方法更为有效。

4.3. 例3

考虑baart (n, 1)问题，选择的扩展子空间为

$W_2=\text{span}\left\{w_1,w_2\right\},w_1={\left(1,1,\cdots ,1\right)}^{\text{T}},w_2={\left(1,2,\cdots ,n\right)}^{\text{T}},w_3={\left(1,4,\cdots ,n^2\right)}^{\text{T}}$

设置n的值为1000，相对噪音水平设置为$\eta ={10}^{-4}$ ，选择迭代次数k为30。通过上述两种方法对该问题进行求解，将解向量的图像与相对误差的结果展示在图7、图8以及表4中。从图7中可以看出，通过Arnoldi-Tikhonov正则化方法得到的解向量与精确解向量有较大差距，得到的解图像只有与精确解图像大致类似的形状，而从图8中可知通过本文方法得到的解图像与精确解图像更加相近，精确程度远好于Arnoldi-Tikhonov正则化方法。在表4中也得出此结论。

Figure 7. The regular solution obtained by A-T algorithm in example 3

图7. 例3中通过A-T算法计算得到的正则化解

Figure 8. The regular solution obtained by E-A-T algorithm in example 3

图8. 例3中通过E-A-T算法计算得到的正则化解

Table 4. Comparing the relative errors of the two algorithms in example 3

表4. 例3中两种算法相对误差数值对比

算法	A-T算法	E-A-T算法
相对误差	1.1177e−01	2.4867e−02

5. 结论

本文融合扩展Arnoldi算法和Tikhonv正则化，给出了求解大规模离散不适定问题的一种扩展迭代正则化算法，即扩展Arnoldi-Tikhonov正则化算法。该算法实现了特定信息补充和投影降维的双重效果，获得更具优势的正则化解。并通过一系列数值算例，展示和验证了所提算例的有效性和优越性。

基金项目

本文得到中央高校基本科研业务费专项资金(NG2023004)资助。

NOTES

^*通讯作者。