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时间(秒)

从表中数据可以清楚显示，广义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方法。

参考文献