关于求解矩阵方程AXB = C的预处理Richardson迭代

On the Preprocessed Richardson Iteration for Solving Matrix Equation AXB = C

摘要: 本文类比于求解线性方程Ax = b的Richardson迭代，通过引入可调参数ω，提出了求解矩阵方程AXB = C的Richardson迭代及其Jacobi和Gauss-Seidel预处理迭代，并详细分析了它们的收敛性。此外，对于一些特殊情况，可以得到参数ω的最优选择，使得迭代矩阵的谱半径达到最小。最后，通过数值实验，我们验证了所提算法的有效性。

Abstract: In this paper, analogous to Richardson iteration for solving linear equation Ax = b, Richardson iteration and preprocessed Jacobi and Gauss-Seidel iteration are proposed for solving the matrix equation AXB = C by introducing a tunable parameter ω, and their convergence properties are analyzed in detail. Moreover, the optimal choices of the parameter ω to minimize the spectral radius of the iteration matrix are also obtained for some special cases. Finally, numerical experiments are carried out to illustrate the effectiveness of the proposed algorithms.

1. 引言

本文考虑求解矩阵方程

$AXB=C,$ (1.1)

其中$A,B,C\in R^{n\times n}$ 。矩阵方程(1.1)在数学、科学计算以及工程方面都有广泛的应用。例如，在计算机辅助几何设计(CAGD)中的曲面拟合、信号与图像处理、摄影测量等许多实际应用中，通常都需要求解方程(1.1)。另外，还包括微分方程数值解、控制系统、稳定性理论、优化算法等领域都涉及到了方程(1.1)的求解。

由于矩阵方程(1.1)在矩阵理论和应用中的重要性，如何有效地求解该方程一直是人们研究的热点。在过去的几十年里，已经提出了许多数值算法来求解该方程。当矩阵A和B阶数较小时，基于QR分解的直接法求解十分有效[1]。然而，对于大型稀疏矩阵A和B，直接法的开销非常昂贵，此时我们通常考虑迭代法来进行求解。在[2]中，作者提出了方程(1.1)的最小二乘对称解的迭代方法，对于任意初始对称矩阵，在不考虑舍入误差的情况下，可以在有限迭代步数内得到一个解，并且通过选择一类特殊的初始对称矩阵，可以得到最小范数解；在[3]中作者讨论了矩阵方程AXB = C的自反解和反自反解，并给出了解存在的一些条件；在[4]中，作者提出了一种迭代算法来获得方程(1.1)的广义中心对称解，并在有限迭代步内推导出给定矩阵的最优逼近解；在[5]中，作者提出了求解方程(1.1)的HSS迭代方法，该方法是通过拓展求解Ax = b的HSS迭代方法得到的；[6]中通过对矩阵A或B进行分裂，构造了求解矩方程(1.1)的Jacobi和Gauss-Seidel型迭代方法，建立了预条件迭代算法，并通过数值实验验证了算法的有效性。

本文通过引入可调参数ω给出了求解方程(1.1)的Richardson迭代法及其Jacobi和Gauss-Seidel预处理迭代方法，并分析了这三种迭代方法的收敛性，最后利用数值实验证明了算法的有效性。

2. 预备知识

下面，我们回顾一些基本定义和已知结果。

定义1. [7]设X为$m\times n$ 矩阵，$vec\left(X\right)$ 定义为大小为$m\cdot n$ 的列向量，是由X的列从左到右堆叠而成。设$A,B$ 分别为$m\times n$ 和$p\times q$ 矩阵。那么克罗内克积$A\otimes B$ 就是$\left(m\cdot p\right)\times \left(n\cdot q\right)$ 矩阵

$\left[\begin{array}{ccc}{a}_{11}\cdot B& \cdots & a_{1n}\cdot B\\ ⋮& & ⋮\\ a_{m1}\cdot B& \cdots & a_{mn}\cdot B\end{array}\right]$ .

引理1. [7]若$A,B,X$ 为$m\times m$ ，$n\times n$ 和$m\times n$ 矩阵，则有：

1) $vec\left(AX\right)=\left(I_n\otimes A\right)\cdot vec\left(X\right)$ ；

2) $vec\left(XB\right)=\left(B^{\text{T}}\otimes I_m\right)\cdot vec\left(X\right)$ 。

引理2. [7]克罗内克积有以下性质：

1) 若$A\cdot C,B\cdot D$ 存在，则$\left(A\otimes B\right)\cdot \left(C\otimes D\right)=\left(A\cdot B\right)\otimes \left(C\cdot D\right)$ ；

2) 若矩阵$A,B$ 可逆，则${\left(A\otimes B\right)}^{-1}=A^{-1}\otimes B^{-1}$ ；

3) ${\left(A\otimes B\right)}^{\text{T}}=A^{\text{T}}\otimes B^{\text{T}}$ 。

引理3. [7]令$\lambda \left(A\right),\mu \left(B\right)$ 分别为矩阵$A\in R^{n\times n},B\in R^{m\times m}$ 的谱，则：

$\lambda \left(A\otimes B\right)=\left\{{\lambda }_i{\mu }_j:{\lambda }_i\in \lambda \left(A\right),{\mu }_j\in \mu \left(B\right),i=1,2,\cdots ,n;j=1,2,\cdots ,m\right\}$ .

定义2. [8] [9]设$Z^{n\times n}$ 表示所有非对角元素非正的$n\times n$ 实矩阵的集合，若矩阵$A\in Z^{n\times n}$ 为非奇异矩阵，且$A^{-1}\ge 0$ ，则称矩阵A为M矩阵。

定义3. [9]对于一个$n\times n$ 矩阵$A=\left(a_{lj}\right)$ ，我们定义一个新的矩阵$\langle A\rangle =\left(a_{lj}\right)$ ，其中

${\alpha }_{lj}=\{\begin{array}{l}\left|{\alpha }_{lj}\right|,j=l\\ -\left|{\alpha }_{lj}\right|,j\ne l\end{array}$ ，

我们称这个矩阵$\langle A\rangle$ 为矩阵A的比较矩阵。

定义4. [9]如果$\langle A\rangle$ 是一个M矩阵，则我们称A是一个H矩阵。

引理4. [8] [9]如果矩阵A为对称正定矩阵或H矩阵，$A=D_1-L_1-U_1$ ，其中$D_1,L_1,U_1$ 分别是矩阵A的对角部分，严格下三角部分和严格上三角部分，则：

$\rho \left({\left(D_1-L_1\right)}^{-1}{U}_1\right)<1$ .

3. 预处理Richardson迭代

3.1. Richardson迭代

对于求解线性方程$Ax=b$ ，Richardson迭代格式如下：

$x^{k+1}=x^k+\omega \left(b-A{x}^k\right)$ .(3.1)

由于矩阵方程$AXB=C$ 可以用Kronecker积符号表示为数学上等价的矩阵向量形式：

$Tx=c$ ，(3.2)

其中$T=B^{\text{T}}\otimes A$ ，$x,c$ 分别表示$vec\left(X\right),vec\left(C\right)$ ，故类似于Richardson迭代我们引入可调参数$\omega$ 得到以下迭代格式：

$x^{k+1}=x^k+\omega \left(c-B^{\text{T}}\otimes A{x}^k\right)$ ，(3.3)

其迭代矩阵为$G_{\omega }=I-\omega B^{\text{T}}\otimes A$ 。但是在具体计算中我们不考虑上述的迭代格式，我们一般考虑与其等价的矩阵迭代格式：

$X^{k+1}=X^k+\omega \left(C-A{X}^{k}B\right)$ .(3.4)

为了加快求解速度，下面我们考虑用Jacobi预处理子和Gauss-Seidel预处理子对矩阵方程(1.1)进行预处理，得到相应的预处理Richardson迭代。

3.2. Jacobi和Gauss-Seidel预处理迭代

在[6]中，提出了求解矩阵方程(1.1)的预处理迭代方法，其描述如下：对于给定的两个预处理矩阵$P_1,P_2$ ，则有：

$AXB{P}_1=C{P}_1$ ，

$P_{2}AXB=P_{2}C$ ，

$P_{2}AXB{P}_1=P_{2}C{P}_1$ .

我们设

$A=D_1-L_1-U_1,B=D_2-L_2-U_2$ ，(3.5)

其中$D_1,D_2$ 分别是矩阵$A,B$ 的对角部分，$L_1,L_2$ 分别是矩阵$A,B$ 的严格下三角部分，$U_1,U_2$ 分别是矩阵$A,B$ 的严格上三角部分。

我们令$P_1=D_1^{-1},P_2=D_2^{-1}$ ，则矩阵方程(1.1)进行Jacobi预处理之后可得：

$AX{B}_1=C_1$ ，

$A_{1}XB=C_2$ ，

$A_{1}X{B}_1=C_3$ .

其中$B_1=B{P}_1,C_1=C{P}_1,A_1=P_{2}A,C_2=P_{2}C,C_3=P_{2}C{P}_1$ 。此时我们由(3.3)可得相应的Jacobi预处理迭代(J-Richardson)为：

$x^{k+1}=\left(I-\omega B_1^{\text{T}}\otimes A\right)x^k+\omega c_1$ ，(3.6)

$x^{k+1}=\left(I-\omega B^{\text{T}}\otimes A_1\right)x^k+\omega c_2$ ，(3.7)

$x^{k+1}=\left(I-\omega B_1^{\text{T}}\otimes A_1\right)x^k+\omega c_3$ .(3.8)

其中$c_1,c_2,c_3$ 分别表示$vec\left(C_1\right),vec\left(C_2\right),vec\left(C_3\right)$ 与迭代(3.8)等价的矩阵迭代格式为：

$X^{k+1}=X^k+\omega \left(C_3-A_1{X}^k{B}_1\right)$ .(3.9)

我们令$P_1={\left(D_1-L_1\right)}^{-1},P_2={\left(D_2-L_2\right)}^{-1}$ ，同理由(3.3)可得相应的Gauss-Seidel预处理迭代(GS-Richardson)为：

$x^{k+1}=\left(I-\omega B_2^{\text{T}}\otimes A\right)x^k+\omega c_4$ ，(3.10)

$x^{k+1}=\left(I-\omega B^{\text{T}}\otimes A_2\right)x^k+\omega c_5$ ，(3.11)

$x^{k+1}=\left(I-\omega B_2^{\text{T}}\otimes A_2\right)x^k+\omega c_6$ .(3.12)

其中$B_2=B{P}_1,A_2=P_{2}A$ ，$c_4,c_5,c_6$ 分别表示$vec\left(C{P}_1\right),vec\left(P_{2}C\right),vec\left(P_{2}C{P}_1\right)$ .与迭代(3.12)等价的矩阵迭代格式为：

$X^{k+1}=X^k+\omega \left(P_{2}C{P}_1-A_2{X}^k{B}_2\right)$ .(3.13)

4. 收敛性分析

4.1. Jacobi预处理迭代收敛性分析

首先我们回顾定理4.1 [10]，即求解线性方程$Ax=b$ 的Richardson迭代的收敛性定理。

定理4.1 若Richardson迭代(3.1)中系数矩阵A的特征值都为正实数，$\alpha$ 满足$0<\alpha <\frac{2}{{\lambda }_{\max }}$ 时，迭代(3.1)收敛，且$\alpha =2/\left({\lambda }_{\min }+{\lambda }_{\max }\right)$ 时，迭代矩阵的谱半径$\rho \left(G_{\alpha }\right)$ 取得最小值$\frac{{\lambda }_{\max }-{\lambda }_{\min }}{{\lambda }_{\min }+{\lambda }_{\max }}$ ，其中$G_{\alpha }=I-\alpha A$ 为迭代矩阵，${\lambda }_{\min },{\lambda }_{\max }$ 为矩阵A的最小特征值与最大特征值。

下面，我们给出Jacobi预条件迭代的收敛性定理。

定理4.2 若矩阵$A,B$ 为由全正基生成的配置矩阵，即全正矩阵，且$\omega$ 满足$0<\omega <\frac{2}{{\lambda }_{\max }}$ 时，Jacobi预

条件迭代(3.8)收敛且可求得最优的$\omega$ 使得迭代矩阵的谱半径达到最小，其中${\lambda }_{\max }$ 为矩阵$T=B_1^{\text{T}}\otimes A_1$ 的最大特征值。

证明. 若矩阵$A,B$ 为全正矩阵，则$P_1=D_1^{-1},P_2=D_2^{-1}$ 为全正矩阵，故$B_1=B{P}_1,A_1=P_{2}A$ 也为全正矩阵，则矩阵$T=B_1^{\text{T}}\otimes A_1$ 的特征值都为正实数，故由定理4.1可知Jacobi预处理迭代(3.8)收敛，且都可求得最优的$\omega$ 使得该迭代的迭代矩阵的谱半径达到最小。

4.2. Gauss-seidel预处理迭代收敛性分析

由于(3.3)还可以看成残差迭代

$x_{k+1}=x_k+\left(c-B^{\text{T}}\otimes A{x}_k\right)$ ，

与外推法

$x_{k+1}=\omega x_{k+1}+\left(1-\omega \right)x_k$ ，

相结合得到的，故迭代(3.3)中的$\omega$ 可以看作外推参数，我们可以类比外推法讨论迭代(3.3)的收敛性定理。下面，我们给出外推法在一般条件下的收敛性定理[11]。

定理4.3 设$S=\left\{{\nu }_j={\gamma }_j+i{\eta }_j,j=1,\cdots ,n\right\}$ 为定常迭代格式：

$X_{k+1}=G{X}_k+C,k=0,1,\cdots ,n$ (4.1)

的迭代矩阵G的谱，${\gamma }_M,{\gamma }_m$ 分别为迭代矩阵G的特征值的最大和最小实部，并且满足$1>{\gamma }_M\ge {\gamma }_j\ge {\gamma }_m$ ，${\eta }_M$ 是G的模最大的虚部，则有：

(1) 设$G_{\omega }=\left(1-\omega \right)I+\omega G$ ，迭代格式(4.1)的外推格式：

$X_{k+1}=G_{\omega }{X}_k+\omega C,k=0,1,\cdots ,n$

收敛当且仅当$0<\omega <\zeta$ ，这里$\zeta =\underset{j}{\min }\left\{\frac{2\left(1-{\gamma }_j\right)}{{\left(1-{\gamma }_j\right)}^2+{\eta }_j^2}\right\}$ 。

设

${\omega }_m=\{\begin{array}{l}{\omega }_1\vartheta \le \psi \\ {\omega }_{*}\vartheta \ge \psi \end{array}$ ，

其中，${\omega }_1=\frac{1-{\gamma }_M}{{\left({\gamma }_M-1\right)}^2+{\gamma }_j^2}$ ，${\omega }^{\ast }=\frac{2}{2-\left({\gamma }_m+{\gamma }_M\right)}$ ，$\vartheta =\left(1-{\gamma }_M\right)\left({\gamma }_M-{\gamma }_m\right)$ ，$\psi =2{\eta }_M^2$ ，我们$\underset{\omega }{\min }\rho \left(G_{\omega }\right)\le \rho \left(G_{{\omega }_m}\right)=\{\begin{array}{ll}\frac{{\eta }_M}{{\left({\left({\gamma }_m-1\right)}^2+{\eta }_M^2\right)}^{1/2}}<1,\hfill & {\omega }_m={\omega }_1\hfill \\ \frac{{\left({\left({\gamma }_M-{\gamma }_m\right)}^2+4{\eta }_M^2\right)}^{1/2}}{2-{\gamma }_M-{\gamma }_m}<1,\hfill & {\omega }_m={\omega }_{\ast }\hfill \end{array}$

当且仅当满足以下条件时等式成立：1) 当$\vartheta \le \psi$ 时，${\gamma }_M+i{\eta }_M\in S$ ；2) 当$\vartheta \ge \psi$ 时，${\gamma }_m+i{\eta }_M,{\gamma }_M+i{\eta }_M\in S$ .

下面，我们给出迭代(3.3)的收敛性定理。

定理4.4 设${\mu }_j={\alpha }_j+i{\beta }_j,j=1,\cdots ,n$ 为$T=B^{\text{T}}\otimes A$ 的所有特征值。${\alpha }_M,{\alpha }_m$ 分别为矩阵T的特征值的最大和最小实部，${\beta }_M$ 为矩阵T的特征值的最大虚部，$\zeta =\underset{j}{\min }\left\{\frac{2{\alpha }_j}{{\alpha }_j^2+{\beta }_j^2}\right\}$ ，$\vartheta ={\alpha }_m\left({\alpha }_M-{\alpha }_m\right)$ ，$\psi =2{\beta }_M^2$ ，${\omega }_1=\frac{{\alpha }_m}{{\alpha }_m^2+{\beta }_M^2}$ ，${\omega }^{\ast }=\frac{2}{{\alpha }_M+{\alpha }_m}$ 。若${\alpha }_j>0,j=1,\cdots ,n$ ，则有：

1) 当$0<\omega <\zeta$ ，迭代(3.3)收敛；

2) 当

${\omega }_m=\{\begin{array}{l}{\omega }_1\vartheta \le \psi \\ {\omega }_{*}\vartheta \ge \psi \end{array}$ ，

时，迭代(3.3)的迭代矩阵$G_{\omega }=I-\omega B^{\text{T}}\otimes A$ 的谱半径上界达到最小值，并且有

$\underset{\omega }{\min }\rho \left(G_{\omega }\right)\le \rho \left(G_{{\omega }_m}\right)=\{\begin{array}{ll}\frac{{\beta }_M}{{\left({\alpha }_m^2+{\beta }_M^2\right)}^{1/2}}<1,\hfill & {\omega }_m={\omega }_1\hfill \\ \frac{{\left({\left({\alpha }_M-{\alpha }_m\right)}^2+4{\beta }_M^2\right)}^{1/2}}{{\alpha }_M+{\alpha }_m}<1,\hfill & {\omega }_m={\omega }_{\ast }\hfill \end{array}$

证明. 设$G=I-B^{\text{T}}\otimes A$ ，${\nu }_j={\gamma }_j+i{\eta }_j,j=1,\cdots ,n$ 矩阵G的特征值，则有${\eta }_j={\beta }_j,{\gamma }_j=1-{\alpha }_j$ 。由于${\alpha }_j>0,j=1,\cdots ,n$ ，故有$1>{\gamma }_M\ge {\gamma }_j\ge {\gamma }_m$ ，其中${\gamma }_M,{\gamma }_m$ 分别为迭代矩阵G的特征值的最大和最小实部，因此，根据定理4.3，我们有(1)和(2)成立。

下面，我们给出Gauss-Seidel预条件迭代的收敛性定理。

定理4.5 若矩阵B(A)为全正矩阵，矩阵A(B)为对称正定矩阵或H矩阵，设${\mu }_j={\alpha }_j+i{\beta }_j,j=1,\cdots ,n$ 为

$T=B^{\text{T}}\otimes A$ 的所有特征值，当$0<\omega <\underset{j}{\min }\left\{\frac{2{\alpha }_j}{{\alpha }_j^2+{\beta }_j^2}\right\}$ 时，单边Gauss-Seidel预处理迭代(3.11) ((3.10))收敛。

证明. 设${\gamma }_l={\sigma }_l+i{\tau }_l,l=1,\cdots ,n$ 为(3.11)中矩阵$A_2={\left(D_2-L_2\right)}^{-1}A$ 的特征值，由于矩阵A为对称正定矩阵或H矩阵，由引理4可知：

$\rho \left(I-{\left(D_2-L_2\right)}^{-1}A\right)=\rho \left({\left(D_2-L_2\right)}^{-1}{U}_2\right)<1$ .

故对于$j=1,\cdots ,n$ 都有

$\left|1-\left({\sigma }_l+i{\tau }_l\right)\right|<1$ ，

则有

${\left(1-{\sigma }_l\right)}^2+{\tau }_l^2<1$ ，

从而可得${\sigma }_l>0$ ，即矩阵$A_2$ 的特征值的实部都大于0。又矩阵B为全正矩阵，即矩阵B的特征值都为正实数，从而矩阵$T=B^{\text{T}}\otimes A$ 的所有特征值的实部都大于0。故由定理4.4可知迭代(3.11)收敛。

定理4.6若矩阵$A,B$ 都为对称正定矩阵或H矩阵，设${\gamma }_l={\sigma }_l+i{\tau }_l,l=1,\cdots ,n$ 为(3.12)中矩阵$A_2={\left(D_2-L_2\right)}^{-1}A$ 的特征值，${\eta }_k={\zeta }_k+i{\upsilon }_k,k=1,\cdots ,n$ 为(3.12)中矩阵$B_2=B{\left(D_1-L_1\right)}^{-1}$ 的特征值，

${\mu }_j={\alpha }_j+i{\beta }_j,j=1,\cdots ,n$ 为$T=B^{\text{T}}\otimes A$ 的所有特征值，当$0<\omega <\underset{j}{\min }\left\{\frac{2{\alpha }_j}{{\alpha }_j^2+{\beta }_j^2}\right\}$ 时，且满足以下(1) (2)中的一点：

1) 对于任意$k,l$ 都有$\left|{\tau }_l\right|<\left|{\sigma }_l\right|,\left|{\upsilon }_k\right|<\left|{\zeta }_k\right|$ .

2) 对于任意$k,l$ 都有${\tau }_l\ast {\upsilon }_k<0$ .

则双边Gauss-Seidel预处理迭代(3.12)收敛。

证明. 由于矩阵$A,B$ 都为对称正定矩阵或H矩阵，同理定理4.5的证明可知矩阵$A_2,B_2$ 的特征值的实部${\sigma }_l,{\zeta }_k$ 都大于0.若(1)或者(2)成立，可得矩阵T的特征值的实部${\alpha }_j={\sigma }_l\ast {\zeta }_k-{\tau }_l\ast {\upsilon }_k>0$ 恒成立，此时由定理4.4可得双边Gauss-Seidel预处理迭代(3.12)收敛。

5. 数值实验

下面，我们给出两个数值示例来说明上述提出的求解矩阵方程(1.1)的迭代方法的性能。在Intel(R)Core(TM) i7-8550U CPU@1.80GHz1.99GHz处理器上，用MATLAB R2016进行了数值实验。我们使用迭代步长(记为IT)、计算时间(记为CPU)和残差(记为RES)三个迭代参数来测试上述的迭代方法，本节所有迭代均从零矩阵开始，当残差范数满足${\Vert R_k\Vert }_F\le {10}^{-8}$ 时停止，其中$R_k=C-A{X}_{k}B$ ，${\Vert G\Vert }_F$ 表示G的Frobenius范数。

例1：考虑矩阵方程(1.1)，其中矩阵A定义为

$A_{10\times 10}=\left(\begin{array}{cccccc}1& & & & & \\ \frac{5}{24}& \frac{7}{12}& \frac{5}{24}& & & \\ & \frac{5}{24}& \frac{7}{12}& \frac{5}{24}& & \\ & & \ddots & \ddots & \ddots & \\ & & & \frac{5}{24}& \frac{7}{12}& \frac{5}{24}\\ & & & & & 1\end{array}\right)$ ，

$B=A^{\text{T}}$ ，C为10阶随机矩阵。

在本例中，我们首先对ω取不同值利用Jacobi预处理迭代(3.9)求解矩阵方程(1.1)，由定理4.2可求出最优参数${\omega }_{\text{opt}}$ 。如表1和图1所示迭代次数和CPU时间方面，最优参数${\omega }_{\text{opt}}$ 表现最好。

其次，我们分别利用Richardson迭代(3.4)，Jacobi预处理迭代(3.9)以及Gauss-Seidel预处理迭代(3.13)求解矩阵方程(1.1)，其中Richardson迭代中的ω取其最优的ω为1.9129，Jacobi预处理迭代中的ω取其最优的ω为0.6800，Gauss-Seidel预处理迭代中的ω取定理4.6中对应的${\omega }_m$ 为1.0000，此时上述三种迭代方法都收敛。我们将三种迭代方法进行比较，可得表2和图2。由表2和图2可知Jacobi预处理迭代方法比广义Richardson迭代方法迭代次数少，收敛速度更快，Gauss-Seidel预处理迭代方法的迭代次数最少，收敛速度最快。

Table 1. Numerical results of preprocessed Jacobi iteration with different values of ω

表1. ω取不同值时Jacobi预处理迭代的数值结果

	ω = 0.2	ω = 0.3	ω = 0.4	ω = 0.5	ω = 0.6	ω = 0.7	ω = ωopt
IT	753	512	389	314	263	446	233
CPU	0.015496	0.011436	0.009033	0.006586	0.006594	0.008983	0.005802
RES	4.81 × 10^−8	3.16 × 10^−8	2.31 × 10^−8	1.79 × 10^−8	1.48 × 10^−8	1.32 × 10^−8	1.27 × 10^−8

Table 2. The number of iteration steps (CPU time) for the three iteration methods under a given error

表2. 三种迭代方法在给定误差下的迭代步数(CPU时间)

算法	ωopt/ωm	IT(CPU/s)
		RES < 10^−^6	RES < 10^−^7	RES < 10^−^8	RES < 10^−^9	RES < 10^−^10	RES < 10^−^11
Richardson	1.929	230 (0.006949)	262 (0.009199)	293 (0.009521)	324 (0.010317)	355 (0.010921)	387 (0.013777)
J-Richardson	0.68	178 (0.003935)	208 (0.005769)	233 (0.005802)	268 (0.0081642)	299 (0.0095736)	329 (0.010421)
GS-Richardson	1	46 (0.001644)	52 (0.004854)	59 (0.005640)	65 (0.006953)	71 (0.007642)	78 (0.008916)

Figure 1. The convergence curve of preprocessed Jacobi iteration for different values of ω

图1. ω取不同值时Jacobi预处理迭代的收敛曲线

Figure 2. Convergence curves of three iterative methods

图2. 三种迭代方法的收敛曲线

类似利用PIA进行曲面拟合，我们可以将本文的算法运用到曲面拟合中，下面给出示例。

例2：考虑曲面

$z=\sin \left(x\right)\ast \cos \left(y\right)$ ，

其中$x\in \left(-2\pi ,2\pi \right),y\in \left(-2\pi ,2\pi \right)$ 。

这里我们考虑矩阵A为由均匀三次b样条基生成的配置矩阵，

$A_{n\times n}=\left(\begin{array}{cccccc}1& & & & & \\ \frac{5}{24}& \frac{7}{12}& \frac{5}{24}& & & \\ & \frac{5}{24}& \frac{5}{24}& \frac{5}{24}& & \\ & & \ddots & \ddots & \ddots & \\ & & & \frac{5}{24}& \frac{5}{24}& \frac{5}{24}\\ & & & & & 1\end{array}\right)$ ，

$B=A^{\text{T}}$ ，选取均匀的插值点，则我们可以得到$n=18$ 时，运用Jacobi预处理迭代(3.9)曲面拟合的结果如图3所示。

Figure 3. Interpolation surface image in Example 2

图3. 例2插值曲面图像

6. 总结

本文首先将矩阵方程AXB = C转化为线性系统Tx = c，通过引入可调参数ω，提出了求解矩阵方程AXB = C的广义Richardson迭代法及其Jacobi和Gauss-Seidel-type预条件迭代法，并分别分析了它们的收敛性。此外，讨论了在这些迭代方法中如何选择最优参数ω，并在一些特殊情况下得到了最优参数。最后，通过数值实验验证了所提算法的有效性。