关于求解矩阵方程AXB = C的广义Richardson迭代及其收敛性

The Generalized Richardson Iteration for Solving the Matrix Equation AXB = C and Its Convergence Analysis

摘要: 本文研究了对于方程AXB = C在传统的Richardson方法基础上，与外推法结合得到的广义Richardson迭代方法。首先，提出广义Richardson迭代方法，然后证明其收敛性。最后，通过数值实验，验证了该迭代方法比传统的渐进迭代逼近法方法(PIA)更有效。

Abstract: This paper investigates the generalized Richardson iterative method for solving the matrix equation AXB = C, which combines extrapolation with the traditional Richardson method. Initially, the generalized Richardson iterative method is proposed, followed by the proof of its convergence. Finally, through numerical experiments, it is verified that this iterative method is more effective than the traditional Progressive Iterative Approximation (PIA) approach.

1. 引言

矩阵方程AXB = C在科学计算和工程领域有着广泛的运用，是一个值得深入研究的问题，在实际问题中也发挥着重要作用，在电力系统潮流计算中，可以通过求解此方程组了解电力网络中各节点的电压以及功率信息；在图像处理技术中，将图像转变成矩阵的形式，通过求解AXB = C这个方程组，可以通过数据恢复原始图像；在计算机辅助几何设计中的曲面拟合中，也有着广泛的运用。而方程是AX = b方程AXB = C的一种特殊情况，当矩阵B为单位矩阵时，C为列向量时，矩阵方程AXB = C就成了Ax = c的形式。因此，如何求解矩阵方程AXB = C引起了国内外广泛学者的关注。

对于求解矩阵方程AXB = C，一般是使用直接法和迭代法进行求解，但是直接法数值不稳定，对于大规模矩阵需要较多的计算和存储资源，而迭代法[1]-[3]是通过逐步逼近方程AXB = C的解，在计算复杂度和存储需求方面较低。现在我们使用的迭代法大多是古典迭代法和投影迭代法，但是迭代法使得系数矩阵的条件数太大了。还有一种解决线性方程组的方法，特别是稀疏矩阵问题和大规模问题中，利用Kronecker积将线性方程组的求解问题转化为更高维度的矩阵方程的求解问题，这样会使得计算量变大。对于矩阵方程AXB = C其Kronecker的形式为：$\left(B^{\text{T}}\otimes A\right)vec\left(X\right)=vec\left(C\right)$ ，该格式在理论研究收敛性中比较实用，但在求解过程中不太符合实际，导致计算量太大。Tian [4]等人发展了一些古典迭代方法，所考虑的分裂是$B^{\text{T}}\otimes A$ 的Jacobi和Gauss-Siedel分裂。

本文通过外推法引入参数$\omega$ 给出了求解方程AXB = C的广义Richardson迭代法，并分析了这钟迭代方法的收敛性，最后利用数值实验证明了算法的有效性。

2. 相关引理

引理1 kronecker积中矩阵A与B的满足的基本运算法则：

数乘运算：$k\left(A\otimes B\right)=kA\otimes B=A\otimes kB$ ，$\forall k\in Ϝ$ ，$Ϝ$ 为数域；

分配率：$\left(A+B\right)\otimes C=A\otimes C+B\otimes C$ ，其中$A,B\in {ℂ}^{m\times n}$ ，$C\in {ℂ}^{p\times q}$ ；

结合率：$\left(A\otimes B\right)\otimes C=A\otimes \left(B\otimes C\right)$ ，其中$A,B\in {ℂ}^{m\times n}$ ，$C\in {ℂ}^{r\times l}$ 。

引理2 设$\mu \left(A\right)$ 和$\lambda \left(B\right)$ 分别是$A\in R_{n\times n}$ 和$B\in R_{m\times m}$ 的谱，则

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

引理3 设$A\in {ℂ}^{m\times n},B\in {ℂ}^{n\times p},C\in {ℂ}^{s\times r},D\in {ℂ}^{r\times l}$ ，则$\left(A\otimes B\right)\left(B\otimes D\right)=AB\otimes CD$ 。

引理4 设$A\in {ℂ}^{m\times m},B\in {ℂ}^{n\times n}$ ，则${\left(A\otimes B\right)}^{\prime }=\left(A^{\prime }\otimes B^{\prime }\right)$ ，如果$A,B$ 的逆矩阵存在，则还有${\left(A\otimes B\right)}^{-1}=A^{-1}\otimes B^{-1}$ 成立。

3. Richardson方法

当解系数矩阵为对称正定阵的大规模线性方程组问题时，我们一般采取Richardson迭代方法，这里我们考虑线性方程Ax = b的求解问题，这里的A是给定的n阶正定矩阵，b是一个n维向量，我们要求出x这个n维向量，这我们可以采取Richardson迭代，其对应的迭代公式为：

$x^{\left(n+1\right)}=x^{\left(n\right)}+\omega \left(b-A{x}^{\left(n\right)}\right),n=0,1,\cdots$

其中$\omega$ 为松弛因子，$x^{\left(0\right)}$ 为初始值，Richardson迭代算法还有一种形式为：

$x^{\left(n+1\right)}=\left(I-\omega A\right)x^{\left(n\right)}+\omega b,n=0,1,\cdots$

此时迭代矩阵$B\left(\omega \right)=I-\omega A$ 。

我们设矩阵A的特征值为$0\le {\lambda }_{\min }={\lambda }_1\le {\lambda }_2\le {\lambda }_3\le \cdots \le {\lambda }_n\le {\lambda }_{n+1}\le {\lambda }_{\max }$ ，则满足当$0<\omega <\frac{2}{{\lambda }_{\max }}$ 时，Richardson迭代收敛，且$\omega ={\omega }_{opt}=\frac{2}{{\lambda }_{\max }+{\lambda }_{\min }}$ 时，Richardson迭代的收敛速度最快。

4. 广义Richardson方法

我们把上述描述的求解线性方程AX = b的Richardson迭代推广到矩阵方程AXB = C上，对于矩阵方程AXB = C的求解问题，这里的$A\in {ℂ}^{n\times n},B\in {ℂ}^{m\times m},C\in {ℂ}^{n\times m}$ ，我们要求出X这个$n\times m$ 的矩阵，这我们提出广义Richardson迭代，其对应的迭代格式为：

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

其可表示成等价的向量形式：

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

其中$x=vec\left(X\right),c=vec\left(C\right)$ 。

5. 广义Richardson的收敛性分析

5.1. A、B特征值为正实数

定理1：当A、B特征值为正实数，$B^{\text{T}}\otimes A$ 的特征值为${\lambda }_i,i=1,\cdots ,n$ ，当$\omega =\frac{2}{{\lambda }_{\max }+{\lambda }_{\min }}$ 时，迭代$x^{k+1}=\left(I-\omega B^{\text{T}}\otimes A\right)x^k+\omega C$ 的谱半径达到最小$\rho =\frac{{\lambda }_{\max }-{\lambda }_{\min }}{{\lambda }_{\max }+{\lambda }_{\min }}$ 。

证明：对于广义的Richardson迭代：

$x^{\left(k+1\right)}=x^{\left(k\right)}+\omega \left(C-A{x}^{\left(k\right)}B\right)$ ，

其可表示成等价的向量形式：

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

其中$x=vec\left(X\right),c=vec\left(C\right)$ .该迭代的迭代矩阵$G_{\omega }=I-\omega B^{\text{T}}\otimes A$ ，收敛因子${\rho }_{\omega }=I-\omega B^{\text{T}}\otimes A$ ，假设$B^{\text{T}}\otimes A$ 的特征值为${\lambda }_i,i=1,\cdots ,n$ ，并且为正实数，${\mu }_i$ 满足

$1-\alpha {\lambda }_{\min }\le {\mu }_i\le 1-\alpha {\lambda }_{\max }$ ，

如果迭代矩阵所有的特征值都为正，即${\lambda }_{\text{min}}>0$ ，则为了使得该迭代收敛，则需要满足下述条件：

$1-\omega {\lambda }_{\max }>-1$ ，

$1-\omega {\lambda }_{\min }<1$ ，

第一个条件意味着$\omega <\frac{2}{{\lambda }_{\max }}$ ，第二个条件意味着$\omega >0$ ，两式结合则可以推出使得迭代收敛应满足：

$0<\omega <\frac{2}{{\lambda }_{\max }}$ ，

设该迭代最优的$\omega$ 为${\omega }_{opt}$ ，当迭代矩阵的谱半径取到最小值时，此时的$\omega$ 为${\omega }_{opt}$ ，其对应的谱半径为：

$\rho G_{\omega }=\max \left\{\left|1-\omega {\lambda }_{\min }\right|,\left|1-\omega {\lambda }_{\max }\right|\right\}$ ，

Figure 1. The function curve of $\rho G_{\omega }$ with respect to $\omega$

图1. $\rho G_{\omega }$ 关于$\omega$ 的函数曲线

如图1所示，关于$\omega$ 的函数当斜率为正的曲线$\left|1-\omega {\lambda }_{\max }\right|$ 与斜率为负的曲线$\left|1-\omega {\lambda }_{\min }\right|$ 相交时，这时取到最优的$\omega$ ，应满足

$1-\omega {\lambda }_{\min }=-1+\omega {\lambda }_{\max }$ ，

可以得到：

${\omega }_{opt}=\frac{2}{{\lambda }_{\min }+{\lambda }_{\max }}$ ，

将其替换为两条曲线中的一条，得到相应最优的谱半径：

${\rho }_{opt}=\frac{{\lambda }_{\text{max}}-{\lambda }_{\text{min}}}{{\lambda }_{\text{min}}+{\lambda }_{\text{max}}}$ ，

5.2. A、B存在复特征值

令$B^{\text{T}}\otimes A$ 的特征值为${\lambda }_i,i=1,\cdots ,n$ ，${\lambda }_i={\alpha }_j+i{\beta }_j,j=1,\cdots ,n$ ，且实部满足${\alpha }_m\le {\alpha }_j\le {\alpha }_M<1,j=1,\cdots ,n$ ，则迭代矩阵$G_{\omega }$ 的特征值为：

${\mu }_j=1-\omega \left(\alpha +\beta i\right)=1-\omega \alpha +\omega \beta i,j=1,\cdots ,n$ ，

迭代矩阵$G_{\omega }$ 的谱半径为：

$\rho \left(G_{\omega }\right)=\max \left|{\mu }_j\right|=\max \left|1-\omega \left(\alpha +\beta i\right)\right|,j=1,\cdots ,n$ ，

定理2：设${\lambda }_i={\alpha }_j+i{\beta }_j,j=1,\cdots ,n$ ，为矩阵$B^{\text{T}}\otimes A$ 的特征值，令

${\alpha }_m={\min }_{1\le j\le n}\left\{{\alpha }_j\right\}$ ，${\alpha }_M={\max }_{1\le j\le n}\left\{{\alpha }_j\right\}$ ，${\beta }_M={\max }_{1\le j\le n}\left|{\beta }_j\right|$ ，${\omega }_1=\frac{{\alpha }_m}{{\alpha }_m^2+{\beta }_M^2}$ ，${\omega }_2=\frac{{\alpha }_M}{{\alpha }_M^2+{\beta }_M^2}$ ，$A={\alpha }_m\left({\alpha }_M-{\alpha }_m\right)$ ，$B=2{\beta }_M^2$ ，$\xi ={\min }_j\left\{\frac{2{\alpha }_j}{{\alpha }_j^2+{\beta }_j^2}\right\}$ ，则有：

1) 当$0<\omega <\xi$ ，广义WPIA收敛。

2) 当${\omega }_{opt}\le \{\begin{array}{l}{\omega }_1,A\le B\\ {\omega }^{\text{*}},A\ge B\end{array}$ ，广义WPIA的迭代矩阵$G_{\omega }$ 的上界取最小值，即

$\min \rho \left(G_{\omega }\right)\le \{\begin{array}{l}\frac{{\beta }_M}{{\left({\alpha }_m^2+{\beta }_M^2\right)}^{\frac{1}{2}}},{\omega }_{opt}={\omega }_1,\\ \frac{{\left[{\left({\alpha }_m-{\alpha }_M\right)}^2+4{\beta }_M^2\right]}^{\frac{1}{2}}}{{\alpha }_m+{\alpha }_M},{\omega }_{opt}={\omega }^{*},\end{array}$

证明：

由${\mu }_j=1-\omega \left(\alpha +\beta i\right)=1-\omega \alpha +\omega \beta i,j=1,\cdots ,n$ ，要使迭代收敛，则对应的谱半径小于1，即$\left|{\mu }_j\right|<1$ ，满足${\left(1-\omega {\alpha }_i\right)}^2+{\left(\omega {\beta }_i\right)}^2<1$ ，解得$0<\omega <\xi$ ，$\xi ={\min }_j\left\{\frac{2{\alpha }_j}{{\alpha }_j^2+{\beta }_j^2}\right\}$ 时，该迭代收敛，(1)得证。

若要得到最优的$\omega$ ，则要找到谱半径上界的最小值，

${\rho }^2\left(G_{\omega }\right)=\max {\left|{\mu }_i\right|}^2=\max \left|{\left(1-\omega \alpha \right)}^2+{\left(\omega \beta \right)}^2\right|\le {\left(\max \left|1-\omega {\alpha }_i\right|\right)}^2+{\omega }^2{\beta }_M^2$ ，

我们只对$\omega >0$ 感兴趣，则设$\max \left|y_i\right|=y_M>0$ ，则可以的得到下述式子：

当${\alpha }_i<\frac{1}{\omega }$ ，则$1-\omega {\alpha }_i>0$ ，

$\max \left|1-\omega {\alpha }_i\right|=\max \left\{1-\omega {\alpha }_i\right\}=1-\omega {\alpha }_m$ ，

当${\alpha }_i>\frac{1}{\omega }$ ，则$1-\omega {\alpha }_i<0$ ，

$\max \left|1-\omega {\alpha }_i\right|=\max \left\{\omega {\alpha }_i-1\right\}=\omega {\alpha }_M-1$ ，

当${\alpha }_m<\frac{1}{\omega }<{\alpha }_M$ ，因此参数满足$\frac{1}{{\alpha }_M}<\omega <\frac{1}{{\alpha }_m}$ ，

$\begin{array}{c}\max \left|1-\omega {\alpha }_i\right|=\max \left\{\max \left\{\omega {\alpha }_i-1\right\},{\alpha }_i\langle \frac{1}{\omega },\max \left\{\omega {\alpha }_i-1\right\},{\alpha }_i\rangle \frac{1}{\omega }\right\}\\ =\max \left\{\omega {\alpha }_M-1,1-\omega {\alpha }_m\right\}\end{array}$

综上所述，可得：

$\max \left|1-\omega {\alpha }_i\right|=\{\begin{array}{l}1-\omega {\alpha }_m,\frac{1}{{\alpha }_M}<\omega \le \frac{2}{{\alpha }_M+{\alpha }_m},\\ \omega {\alpha }_M-1,\frac{2}{{\alpha }_M+{\alpha }_m}<\omega \le \frac{1}{{\alpha }_M},\end{array}$

下面我们记${\omega }^{*}=\frac{2}{{\alpha }_M+{\alpha }_m}$ ，则可以得到：

${\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{l}{\left(1-\omega {\alpha }_m\right)}^2+{\omega }^2{\beta }_M^2,0<\omega \le {\omega }^{*},\\ {\left(\omega {\alpha }_M-1\right)}^2+{\omega }^2{\beta }_M^2,{\omega }^{*}<\omega ,\end{array}$

现在求出使得${\rho }^2\left(G_{\omega }\right)$ 最小的$\omega$ ：

$\min {\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{l}\min g_1\left(\omega \right),0<\omega \le {\omega }^{*},\\ \min g_2\left(\omega \right),{\omega }^{*}<\omega ,\end{array}$

为了寻找${\min }_{0<\omega \le {\omega }^{*}}{g}_1\left(\omega \right)$ 和${\min }_{{\omega }^{*}<\omega }{g}_2\left(\omega \right)$ ，我们可以考虑对$g_1\left(\omega \right)$ ，$g_2\left(\omega \right)$ 进行求导可得${g^{\prime }}_1\left(\omega \right)\ge 0\left(\omega >{\omega }_1\right)$ 和${g^{\prime }}_2\left(\omega \right)\ge 0\left(\omega >{\omega }_2\right)$ ，其中：

${\omega }_1=\frac{{\alpha }_m}{{\alpha }_m^2+{\beta }_M^2},{\omega }_2=\frac{{\alpha }_M}{{\alpha }_M^2+{\beta }_M^2},$

现在我们可以得到结论，如果$\omega \le {\omega }_1$ ，则满足：

${\alpha }_m\left({\alpha }_M-{\alpha }_m\right)\le 2{\beta }_M^2$ .

当$0<{\omega }_1\le {\omega }^{*}$ ，$g_1\left(\omega \right)$ 在$0<{\omega }_1$ 是单调递减的，$g_1\left(\omega \right)$ 在${\omega }_1\le {\omega }^{*}$ 是单调递增的，所以在${\omega }_1$ 处取得最小，即当$0<{\omega }_1\le {\omega }^{*}$ 时，${\min }_{0<{\omega }_1\le {\omega }^{*}}{g}_1\left(\omega \right)=g_1\left({\omega }_1\right)$ 。

当${\omega }_1\ge {\omega }^{*}$ 时，$g_1\left(\omega \right)$ 是单调递减的，所以在${\omega }^{*}$ 取得最小，即当${\omega }_1\ge {\omega }^{*}$ 时，${\min }_{{\omega }_1\ge {\omega }^{*}}{g}_1\left(\omega \right)=g_1\left({\omega }^{*}\right)$ 成立。

当${\omega }_2\le {\omega }^{*}$ ，$g_2\left(\omega \right)$ 是一个单调递增函数，所以在${\omega }^{*}$ 取得最小，即当$\omega \ge {\omega }^{*}$ 时，${\min }_{\omega \ge {\omega }^{*}}{g}_2\left(\omega \right)=g_2\left({\omega }^{*}\right)=g_1\left({\omega }^{*}\right)$ 成立。

综上所述，可总结为：

$\min {\rho }^2\left(G_{\omega }\right)\le \{\begin{array}{l}{g}_1\left(\omega \right),A\le B,\\ g_1\left({\omega }^{*}\right),A\ge B,\end{array}$

其中$A={\alpha }_m\left({\alpha }_M-{\alpha }_m\right)$ ，$B=2{\beta }_M^2$ 则可以得到：

$\min \rho \left(G_{\omega }\right)\le \{\begin{array}{l}\frac{{\beta }_M}{{\left({\alpha }_m^2+{\beta }_M^2\right)}^{\frac{1}{2}}},A\le B,\\ \frac{{\left[{\left({\alpha }_m-{\alpha }_M\right)}^2+4{\beta }_M^2\right]}^{\frac{1}{2}}}{{\alpha }_m+{\alpha }_M},A\ge B,\end{array}$

对于上述的$\min \rho \left(G\right)$ 的上界，我们可以作下述说明：

1) 当$A\le B$ 时，显然上界总是小于1，且以下式子成立：

$\min \rho \left(G_{\omega }\right)=g_1\left({\omega }_1\right)=\frac{{\beta }_M}{{\left({\alpha }_m^2+{\beta }_M^2\right)}^{\frac{1}{2}}}<1$ ，

故在${\alpha }_m+i{\beta }_M$ 为G的特征值时：

$\min \rho \left(G_{\omega }\right)=g_1\left({\omega }_1\right)=\frac{{\beta }_M}{{\left({\alpha }_m^2+{\beta }_M^2\right)}^{\frac{1}{2}}}$ ，${\omega }_{opt}={\omega }_1$ ，

2) 当$A\ge B$ ，即${\alpha }_m{\alpha }_M>2{\beta }_M^2\ge {\beta }_M^2$ ，故可以得到：

${\left({\alpha }_m-{\alpha }_M\right)}^2+4{\beta }_M^2<\left({\alpha }_m^2+{\beta }_M^2\right)$ ，

那么可以得到：

$g_1\left({\omega }^{*}\right)=\frac{{\left[{\left({\alpha }_m-{\alpha }_M\right)}^2+4{\beta }_M^2\right]}^{\frac{1}{2}}}{{\alpha }_m+{\alpha }_M}<1$ ，

故在${\alpha }_M+i{\beta }_M$ 以及${\alpha }_m+i{\beta }_M$ 为G的特征值时：

$\min \rho \left(G_{\omega }\right)=g_1\left({\omega }^{*}\right)=\frac{{\left[{\left({\alpha }_m-{\alpha }_M\right)}^2+4{\beta }_M^2\right]}^{\frac{1}{2}}}{{\alpha }_m+{\alpha }_M}$ ，${\omega }_{opt}={\omega }^{*}$ .

5.3. 数值实验

下面我们考虑将广义Richardson迭代运用到曲面拟合[5]，首先我们回顾一下张量积曲面拟合的传统方法的细节，设两个标准化的正交积为$\left({\mu }_0,\cdots ,{\mu }_n\right)$ 和$\left({\nu }_0,\cdots ,{\nu }_n\right)$ 。给定一组控制点序列为$\left\{P_{ij}\right\},i=0,1,\cdots ,m;j=0,1,\cdots ,n$ ，则生成的张量基曲面可以表示为：

$S\left(\mu ,\nu \right)={\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^n{P}_{ij}{B}_i\left(\mu \right)B_j\left(\nu \right),$

且参数值满足$t_0<t_1<\cdots <t_m$ 和$s_0<s_1<\cdots <s_n$ 这两组实递增序列，将其对应的参数对$\left(t_i,s_j\right)$ 赋予控制点$\left\{P_{ij}\right\}$ ，从$i=0,1,\cdots ,m;j=0,1,\cdots ,n$ ，这个点列为初始控制点，则对应生成的初始张量积曲面为：

$S^{\left(0\right)}\left(t,s\right)={\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^n{P}_{ij}^{\left(0\right)}{\mu }_i\left(t\right){\nu }_j\left(s\right).$

则进行第K + 1次迭代后可生成的曲面为：

$S^{\left(k+1\right)}\left(t,s\right)={\displaystyle \sum }_{i=0}^m{\displaystyle \sum }_{j=0}^n{P}_{ij}^{\left(k+1\right)}{\mu }_i\left(t\right){\nu }_j\left(s\right).$

对于K + 1次迭代，计算其对应的校正向量：

$\{\begin{array}{l}{\Delta }_{uv}^{\left(k\right)}=P_{uv}-S^{\left(k\right)}\left(t_u,s_v\right),\\ P_{uv}^{\left(k+1\right)}=P_{uv}^{\left(k\right)}+{\Delta }_{uv}^{\left(k\right)},\end{array}$ (4.1)

其中，$u=0,1,\cdots ,n;v=0,1,\cdots ,m$ .则可以得到

$P_{uv}^{\left(k+1\right)}=P_{uv}^{\left(k\right)}+P_{uv}-S^{\left(k\right)}\left(t_u,s_v\right),$

该迭代称为渐进迭代逼近，简称为PIA。

将$P_{uv}^{\left(k\right)}$ 排列成矩阵形式可以得到：

$P^{\left(k\right)}=\left(P_{00}^{\left(k\right)},P_{01}^{\left(k\right)},\cdots ,P_{0m}^{\left(k\right)},P_{10}^{\left(k\right)},P_{11}^{\left(k\right)},\cdots ,P_{1m}^{\left(k\right)},\cdots ,P_{n0}^{\left(k\right)},P_{n1}^{\left(k\right)},\cdots ,P_{nm}^{\left(k\right)}\right)$

如果$\underset{k\to \infty }{\lim }{S}^{\left(k\right)}\left(t_u,S_v\right)=P_{uv}$ .则初始的$S^{\left(0\right)}\left(t_u,S_v\right)$ 有PIA性质。

设A，B矩阵分别为正交积为$\left({\mu }_0,\cdots ,{\mu }_n\right)$ 和$\left({\nu }_0,\cdots ,{\nu }_n\right)$ 的配置矩阵，则(2.1)可以写成：

$P^{\left(k+1\right)}=P^{\left(k\right)}+\left[P-\left(B\otimes A\right)P^{\left(k\right)}\right]$ ，(4.2)

由上述的(2.2)式，我们也可以将其写成下面形式，

$\{\begin{array}{l}{P}_x^{\left(k+1\right)}=P_x^{\left(k\right)}+\left[P_x-\left(B\otimes A\right)P_x^{\left(k\right)}\right],\\ P_y^{\left(k+1\right)}=P_y^{\left(k\right)}+\left[P_y-\left(B\otimes A\right)P_y^{\left(k\right)}\right],\\ P_z^{\left(k+1\right)}=P_z^{\left(k\right)}+\left[P_z-\left(B\otimes A\right)P_z^{\left(k\right)}\right],\end{array}$ (4.3)

迭代(2.3)在数学上则等价于Richardson迭代，可以用来求解下列方程组：

${\left(B\otimes A\right)}_x=P_x$ ，

${\left(B\otimes A\right)}_y=P_y$ ，

${\left(B\otimes A\right)}_z=P_z$ ，

看成数值线性代数问题，最好就是求解下面三个方程：

$AX{B}^{\prime }=C_1$ ，

$AX{B}^{\prime }=C_2$ ，(4.4)

$AX{B}^{\prime }=C_3$ ，

其中，$vec\left(X\right)=x,vec\left(Y\right)=y,vec\left(Z\right)=z,vec\left(C_1\right)=P_x,vec\left(C_2\right)=P_y,vec\left(C_3\right)=P_z$ .

这里我们比较PIA算法和广义Richardson，刘成志[6]等人提出的三次均匀B样条扩展曲面的渐进迭代逼近法得出配置矩阵与形状参数λ有关，这里我们取A，B矩阵为H-矩阵的三次B样条基，具体形式为：

$A=\left(\begin{array}{ccccc}1& 0& & & \\ \frac{4-\lambda }{24}& \frac{8-\lambda }{12}& \frac{4-\lambda }{24}& \cdots & \\ & & \frac{4-\lambda }{24}& \frac{8-\lambda }{12}& \frac{4-\lambda }{24}\\ & & & 0& 1\end{array}\right),B=A^{\text{T}}$

其中$-2\le \lambda \le 1$ ，这里我们让$\lambda =-1$ ，求解了(4.4)对应的三个线性方程，要求的精度为10⁻⁸，结果如下表所示(表1)。

Table 1. CPU times (seconds) for PIA and generalized richardson

表1. PIA与广义Richardson的CPU时间(秒)

n	PIA	广义Richardson
n = 10	0.003	0.002
n = 20	0.057	0.01
n = 30	0.348	0.03
n = 40	1.668	0.036
n = 50	5.728	0.037
n = 60	17.49	0.07
n = 70	45.6	0.072
n = 80	104.68	0.107

从表中数据可以清楚显示，广义Richardson明显快于PIA。当n的取值越大，广义Richardson的收敛速度越好，当$n=80$ 时，PIA方法所以需要的时间达到了10,468秒，而广义Richardson仅仅需要0.107秒。下面我们展示$\lambda =-1,n=10$ 和$\lambda =-1,n=30$ 的两个插值曲面(图2)。

(a) (b)

Figure 2. Interpolated surface plots when $\lambda = -1,n=10$ (left) and $\lambda =-1,n=30$ (right)

图2. 当$\lambda = -1,n=10$ (左) $\lambda =-1,n=30$ (右)时插值曲面的图像

6. 结论

本文给出了解方程$AXB=C$ 的一种算法，该算法将外推法与Richardson方法相结合，通过寻找$\omega$ 的取值来提高收敛速度。并且对A、B矩阵分为特征值为实数和复数两种情况并且分别做了收敛性分析，最后利用曲面拟合来研究PIA算法和广义Richardson算法的速度。数值实验表明：广义Richardson方法明显快于PIA方法。