位移全局GPBiCG方法求解Stein矩阵方程X + AXB = C

Shifted Global GPBiCG Method for Solving the Stein Matrix Equation X + AXB = C

摘要: 广义点积型双共轭梯度法(Generalized Product-type Biconjugate Gradient, GPBiCG)作为共轭梯度平方法(Conjugate Gradient Squared, CGS)和稳定双共轭梯度法(Biconjugate Gradient Stabilized, BiCGStab)的推广，收敛速度快于BiCGStab，收敛曲线比CGS更光滑。因此全局型GPBiCG方法被用于求解具有多个右端项的线性系统。基于全局GPBiCG方法，本文提出一种位移全局GPBiCG方法(Shifted Global GPBiCG, SGl-GPBiCG)求解Stein矩阵方程。该方法充分利用了Stein矩阵方程的位移结构。最后，我们给出数值算例验证了新方法的有效性。

Abstract: Generalized Product-type Biconjugate Gradient (GPBiCG) method can be regarded as generalizations of Conjugate Gradient Squared (CGS) and Biconjugate Gradient Stabilized (BiCGStab) methods which is faster than BiCGStab and its convergence is smoother than CGS. Then the global variant of GPBiCG method is proposed for linear systems with multiple right-hand sides. In this paper, we present a shifted global GPBiCG (SGl-GPBiCG) method for solving the Stein matrix equation based on the global GPBiCG method. The new method makes full use of the shifted structure of the Stein matrix equation. Finally, numerical examples are given to illustrate the effectiveness of the method.

1. 引言

Stein矩阵方程在图像复原[1] [2]，控制系统设计[3] [4]，神经网络[5]，模型降阶[6]等领域有着广泛的应用。Stein矩阵方程又称为离散型的Sylvester方程，其标准形式为

$X+AXB=C,$ (1.1)

其中，$A\in R^{n\times n},B\in R^{s\times s},C\in R^{n\times s}$ ，其存在唯一解的充要条件为$\lambda \left(A\right)\lambda \left(B\right)\ne 1$ ，$\lambda \left(A\right),\lambda \left(B\right)$ 分别表示矩阵$A$ 和$B$ 的所有特征值的集合。本文只考虑Stein矩阵方程有唯一解的情况。

Krylov子空间方法[7]是求解大型线性方程组$Ax=b$ 的有效方法，相应的块Krylov子空间方法[8]-[11]和全局Krylov子空间方法[12]-[14]是求解具有多右端项的线性系统$AX=B$ 的有效方法。近年来，块方法和全局方法被推广用于求解Stein矩阵方程，其中系数矩阵$A$ 和$B$ 的维数不同是构造不同方法的关键。当矩阵$A$ 和$B$ 的维数都小时，直接法是最有效的选择，比如Hessenberg-Schur方法[15]和Bartels-Stewart 方法[16]。当矩阵$A$ 很大，$B$ 很小时，黄等人[17]提出块Arnoldi Stein方法和非对称的块Lanczos Stein方法求解Stein矩阵方程。文献[18] [19]提出全局Arnoldi方法求Stein矩阵方程的低秩解。文献[20]提出块Arnoldi方法求其低秩解。但是，上述这些方法都没有利用到Stein矩阵方程的特殊位移结构。当矩阵$A$ 和$B$ 都很大时，Smith方法[21]被提出用于Stein矩阵方程的求解。然而，Smith方法要求矩阵$A$ 和$B$ 的谱半径满足$\rho \left(A\right)\rho \left(B\right)<1$ 。本文主要研究大规模Stein矩阵方程的求解。

在文献[22]中，张邵良提出了一种广义点积型双共轭梯度法(Generalized Product-type Biconjugate Gradient, GPBiCG)求解线性方程组。在文献[23]中，Dehghan等人提出了一种位移GPBiCG方法求解位移线性系统$Ax+\delta x=b$ 。这种位移技巧来源于Frommer的工作[24]。受到上述方法的启发，本文利用Stein矩阵方程的位移结构，基于全局GPBiCG方法提出了一种位移全局GPBiCG方法(Shifted Global GPBiCG, SGl-GPBiCG)。该方法适合任意维数的系数矩阵$A$ 和$B$ ，不受矩阵维数的影响。

为了推导新方法，首先将矩阵方程(1.1)命名为位移系统，将矩阵方程

$AXB=C,$ (1.2)

命名为种子系统。为了方便起见，定义线性算子${\rm M}$ 如下

${\rm M}:{\rm M}\left(X\right)=AXB,$

相应地，种子系统(1.2)可以改写成

${\rm M}\left(X\right)=C,$ (1.3)

而位移系统(1.1)可写成

$\left({\rm M}+I\right)X=C.$ (1.4)

对给定的矩阵$V\in R^{n\times s}$ 和算子${\rm M}$ ，定义矩阵Krylov子空间如下

$K_m\left(M,V\right)=span\left\{V,M\left(V\right),\cdots ,M^{m-1}\left(V\right)\right\},$ (1.5)

其中

$M^i\left(V\right)=M\left(M^{i-1}\left(V\right)\right).$

由于Krylov子空间的位移不变性，我们有

$K_m\left(M,V\right)=K_m\left(\left(M+I\right),V\right).$ (1.6)

因此，在初始近似解为零矩阵的条件下，种子系统与位移系统在第$m$ 次迭代时得到的近似解在同一个矩阵Krylov子空间中，这为后续方法的推导提供了关键的理论支撑。

2. 全局GPBiCG算法

在这节，我们借鉴文献[25]中求解$AX=B$ 的全局GPBiCG方法，直接以算子形式给出求解种子系统(1.3)的全局GPBiCG算法(Gl-GPBiCG)，具体迭代过程见表1。

Table 1. Global GPBiCG method for solving seed system (1.3)

表1. 全局GPBiCG算法求解种子系统(1.3)

算法2.1：全局GPBiCG算法求解种子系统(1.3)

给定初值$X_0\in R^{m\times n}$ ，计算$R_0=C-{\rm M}\left(X_0\right),$ 选择$R_0^{*}=R_0$ ，满足$\langle R_0^{*},R_0\rangle \ne 0,$ 令$T_{-1}=W_{-1}=0_{N\times S},{\beta }_{-1}=0,$ 对于$n=0,1,\cdots$ 计算 $P_n=R_n+{\beta }_{n-1}\left(P_{n-1}-U_{n-1}\right),$ ${\alpha }_n=\frac{{\langle R_0^{*},R_n\rangle }_F}{{\langle R_0^{*},{\rm M}{P}_n\rangle }_F},$ $T_n=R_n-{\alpha }_n{\rm M}{P}_n,$ $Y_n=T_{n-1}-R_n-{\alpha }_n{W}_{n-1}+{\alpha }_n{\rm M}{P}_n,$ ${\zeta }_n==\frac{{\langle Y_n,Y_n\rangle }_F{\langle {\rm M}{T}_n,T_n\rangle }_F-{\langle Y_n,T_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F}{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,Y_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F},$ ${\eta }_n=\frac{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,T_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,T_n\rangle }_F}{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,Y_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F},$ 如果$n=0$ ，则${\zeta }_n=\frac{\langle {\rm M}{T}_n,T_n\rangle }{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle },{\eta }_n=0.$ $U_n={\zeta }_n{\rm M}{P}_n+{\eta }_n\left(T_{n-1}-R_n+{\beta }_{n-1}{U}_{n-1}\right),$ $Z_n={\zeta }_n{R}_n+{\eta }_n{Z}_{n-1}-{\alpha }_n{U}_n,$ $R_{n+1}=T_n-{\eta }_n{Y}_n-{\zeta }_n{\rm M}{T}_n,$ ${\beta }_n=\frac{{\alpha }_n}{{\zeta }_n}\cdot \frac{{\langle R_0^{*},R_{n+1}\rangle }_F}{{\langle R_0^{*},R_n\rangle }_F},$ $X_{n+1}=X_n+{\alpha }_n{P}_n+Z_n,$ $W_n={\rm M}{T}_n+{\beta }_n{\rm M}{P}_n,$

结束.

在上述算法中，全局GPBiCG方法是通过残差矩阵优化近似解，残差矩阵表达式为

$R_n=H_n\left({\rm M}\right){\varphi }_n\left({\rm M}\right)R_0,$ (2.1)

其关键在于残差多项式${\varphi }_n$ 与辅助矩阵$H_n$ 的构建。回顾张绍良[22]提出的GPBiCG方法，${\varphi }_n$ 是BiCG的残差多项式，残差多项式${\varphi }_n$ 和辅助多项式${\phi }_n$ 之间的基本递归关系如下

${\varphi }_0\left(t\right)=1,{\phi }_0\left(t\right)=1,$ (2.2)

${\varphi }_{n+1}\left(t\right)={\varphi }_n\left(t\right)-{\alpha }_{n}t{\phi }_n\left(t\right).$ (2.3)

消去${\phi }_n$ ，得到${\varphi }_n$ 的单独的三项递归关系

${\varphi }_0\left(t\right)=1,{\varphi }_1\left(t\right)=\left(1-{\alpha }_{0}t\right){\varphi }_0\left(t\right),$ (2.4)

${\varphi }_{n+1}\left(t\right)=-\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}{\varphi }_{n-1}\left(t\right)+\left(1+\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}\right){\varphi }_n\left(t\right)-{\alpha }_{n}t{\varphi }_n\left(t\right).$ (2.5)

且辅助矩阵$H_n$ 满足如下三项递归关系

$H_0\left(t\right)=1,H_1\left(t+1\right)=\left(1-{\zeta }_0\left(t\right)\right)H_0\left(t\right),$ (2.6)

$H_{n+1}\left(t\right)=\left(1+{\eta }_n-{\zeta }_n\left(t\right)\right)H_n\left(t\right)-{\eta }_n{H}_{n-1}\left(t\right).$ (2.7)

3. 位移全局GPBiCG算法

在本节中，我们在全局GPBiCG方法求解种子系统(1.3)的基础上，提出求解位移系统(1.4)的位移全局GPBiCG方法。

3.1. 位移系统的残差

众所周知，GPBiCG方法是在BiCG方法的残差多项式的基础上建立新的三项递推关系得到的。基于以上思想，我们推导求解位移系统的位移全局GPBiCG方法。

考虑种子系统的残差多项式的递归关系，定义位移系统的残差多项式为

${\varphi }_0^s\left(t+1\right)=1,{\varphi }_1^s\left(t+1\right)=\left(1-{\alpha }_0\left(t+1\right)\right){\varphi }_0^s\left(t+1\right),$ (3.1)

${\varphi }_{n+1}^s\left(t+1\right)=-\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}{\varphi }_{n-1}^s\left(t+1\right)+\left(1+\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}\right){\varphi }_n^s\left(t+1\right)-{\alpha }_n^s\left(t+1\right){\varphi }_n^s\left(t+1\right).$ (3.2)

定义辅助多项式

${\phi }_0^s\left(t+1\right)=1,$ (3.3)

${\phi }_n^s\left(t+1\right)=\frac{{\varphi }_n^s\left(t+1\right)-{\varphi }_{n+1}^s\left(t+1\right)}{{\alpha }_n^s\left(t+1\right)}.$ (3.4)

并且通过对(3.2)和(3.4)的一些简单运算，得到以下两个递归关系

${\varphi }_0^s\left(t+1\right)=1,{\phi }_0\left(t\right)=1,$ (3.5)

${\varphi }_{n+1}^s\left(t+1\right)={\varphi }_n^s\left(t+1\right)-{\alpha }_n^s\left(t+1\right){\phi }_n^s\left(t+1\right),$ (3.6)

${\phi }_{n+1}^s\left(t+1\right)={\varphi }_{n+1}^s\left(t+1\right)+{\beta }_n^s{\phi }_n^s\left(t+1\right).$ (3.7)

考虑前面辅助矩阵$H_n$ 的三项递归关系，定义位移系统的辅助矩阵

$H_0^s\left(t+1\right)=1,H_1^s\left(t+1\right)=\left(1-{\zeta }_0^s\left(t+1\right)\right)H_0^s\left(t+1\right),$ (3.8)

$H_{n+1}^s\left(t+1\right)=\left(1+{\eta }_n^s-{\zeta }_n^s\left(t+1\right)\right)H_n^s\left(t+1\right)-{\eta }_n^s{H}_{n-1}^s\left(t+1\right).$ (3.9)

定义

$G_0^s\left(t+1\right)=1,$ (3.10)

$G_n^s\left(t+1\right)=\frac{H_n^s\left(t+1\right)-H_{n+1}^s\left(t+1\right)}{{\zeta }_n^s\left(t+1\right)}.$ (3.11)

然后用这个定义对(3.9)进行一些简单计算，得到辅助矩阵$H_n^s$ 的两项递归关系

$H_{n+1}^s\left(t+1\right)=H_n^s\left(t+1\right)-{\zeta }_n^s\left(t+1\right)G_n^s\left(t+1\right),$ (3.12)

$\begin{array}{c}{G}_{n+1}^s\left(t+1\right)=H_{n+1}^s\left(t+1\right)+{\zeta }_n^s\frac{{\eta }_{n+1}^s}{{\zeta }_{n+1}^s}{G}_n^s\left(t+1\right)G_{n+1}^s\left(t+1\right)\\ =H_{n+1}^s\left(t+1\right)+{\zeta }_n^s\frac{{\eta }_{n+1}^s}{{\zeta }_{n+1}^s}{G}_n^s\left(t+1\right).\end{array}$ (3.13)

现在计算位移系统残差$R_n^s=H_n^s\left({\rm M}\right){\varphi }_n^s\left({\rm M}\right)R_0$ ，使用与全局GPBiCG相同的推导方法，并定义新的辅助矩阵如下

$\{\begin{array}{l}{T}_n^s=H_n^s{\varphi }_{n+1}^s,Y_n^s=\left({\rm M}+I\right)G_{n-1}^s{\varphi }_{n+1}^s,P_n^s=H_n^s{\phi }_n^s\hfill \\ W_n^s=\left({\rm M}+I\right)H_n^s{\phi }_{n+1},U_n^s=\left({\rm M}+I\right)G_n^s{\phi }_n^s,Z_n^s=G_n^s{\varphi }_n^s.\hfill \end{array}$ (3.14)

得到中间矩阵的递归关系如下

$P_n^s=R_n^s+{\beta }_n^s\left(P_{n-1}^s-U_{n-1}^s\right),$

$Y_n^s=T_{n-1}^s-R_n^s-{\alpha }_n^s{W}_{n-1}^s+{\alpha }_n^s\left({\rm M}+I\right)P_n^s,$

$T_n^s=R_n^s-{\alpha }_n^s\left({\rm M}+I\right)P_n^s,$

$U_n^s={\zeta }_n^s\left({\rm M}+I\right)P_n^s+{\eta }_n^s\left(T_{n-1}^s-R_n^s+{\beta }_{n-1}^s{U}_{n-1}^s\right),$

$Z_n^s={\zeta }_n^s{R}_n^s+{\eta }_n^s{Z}_{n-1}^s-{\alpha }_n^s{U}_n^s,$

$X_{n+1}^s=X_n^s+{\alpha }_n^s{P}_n^s-Z_n^s,$

$R_{n+1}^s=T_n^s-{\eta }_n^s{Y}_n^s-{\zeta }_n^s\left({\rm M}+I\right)T_n^s,$

$W_n^s=\left({\rm M}+I\right)T_n^s+{\beta }_n^s\left({\rm M}+I\right)P_n^s.$

比较算法2.1中对应式子不难发现，全局CPBiCG直接求解位移系统需要额外计算$\left(M+I\right)$ 与中间矩阵的乘积，相较于求解种子系统，增加了一次矩阵乘矩阵运算。

3.2. 残差的共线性

除了Krylov子空间的位移不变性，残差的共线性也是实现位移方法的一个关键因素。根据文献[24]中的定理1知，在种子系统和位移系统的初始近似解为0的前提下，两系统的残差共线。接下来，我们借鉴文献[23]的推导方式，分别为残差多项式${\varphi }_n\left(t\right)$ ，${\varphi }_n^s\left(t+1\right)$ 和辅助多项式$H_n\left(t\right)$ ，$H_n^s\left(t+1\right)$ 建立共线性。

首先，引入系数${\delta }_n$ ，令

${\varphi }_n\left(t\right)={\delta }_n{\varphi }_n^s\left(t+1\right).$ (3.15)

从(2.5)有

${\varphi }_{n+1}^s\left(t+1\right)=-\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}{\varphi }_{n-1}^s\left(t+1\right)+\left(1+\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}-{\alpha }_n^s\right){\varphi }_n^s\left(t+1\right)-{\alpha }_n^s\left(t+1\right){\varphi }_n^s\left(t+1\right),$

利用(3.15)式，将${\varphi }_{n+1}^s\left(t+1\right),{\varphi }_n^s\left(t+1\right),{\varphi }_{n-1}^s\left(t+1\right)$ 分别用

$\left(\frac{1}{{\delta }_{n+1}}\right){\varphi }_{n+1}\left(t\right),\left(\frac{1}{{\delta }_n}\right){\varphi }_n\left(t\right),\left(\frac{1}{{\delta }_{n-1}}\right){\varphi }_{n-1}\left(t\right),$

代替，进而得到

$\left(\frac{1}{{\delta }_{n+1}}\right){\varphi }_{n+1}\left(t\right)=-\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}\left(\frac{1}{{\delta }_{n-1}}\right){\varphi }_{n-1}\left(t\right)+\left(1+\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}-{\alpha }_n^s\right)\left(\frac{1}{{\delta }_n}\right){\varphi }_n\left(t\right)-{\alpha }_n^s\left(t+1\right)\left(\frac{1}{{\delta }_n}\right){\varphi }_n\left(t\right),$

整理一下得到

${\varphi }_{n+1}\left(t\right)=-\frac{{\delta }_{n+1}}{{\delta }_{n-1}}\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}{\varphi }_{n-1}\left(t\right)+\frac{{\delta }_{n+1}}{{\delta }_n}\left(1+\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}-{\alpha }_n^s\right){\varphi }_n\left(t\right)-\frac{{\delta }_{n+1}}{{\delta }_n}{\alpha }_n^{s}t{\varphi }_n\left(t\right).$ (3.16)

比较(2.5)和(3.16)，有

$\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}=\frac{{\delta }_{n+1}}{{\delta }_{n-1}}\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s},$ (3.17)

$1+\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}=\frac{{\delta }_{n+1}}{{\delta }_n}\left(1+\frac{{\alpha }_n^s{\beta }_{n-1}^s}{{\alpha }_{n-1}^s}-{\alpha }_n^s\right),$ (3.18)

${\alpha }_n^s=\frac{{\delta }_n}{{\delta }_{n+1}}{\alpha }_n.$ (3.19)

将(3.19)代入(3.18)得到系数${\beta }_n^s$ 的转化公式

${\beta }_n^s={\left(\frac{{\delta }_n}{{\delta }_{n+1}}\right)}^2{\beta }_n,$ (3.20)

将(3.20)和(3.19)代入(3.18)并化简，得到辅助系数${\delta }_n$ 一个三项递归关系

${\delta }_{n+1}=\left(1+{\alpha }_n\right){\delta }_n+\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}\left({\delta }_n-{\delta }_{n-1}\right).$ (3.21)

同理引入系数${\xi }_n$ ，设

$H_n={\xi }_n{H}_n^s.$ (3.22)

将(3.22)代入(3.10)化简得到

$H_{n+1}\left(t\right)=-\frac{{\xi }_{n+1}}{{\xi }_{n-1}}{\eta }_n^s{H}_{n-1}\left(t\right)+\frac{{\xi }_{n+1}}{{\xi }_n}\left(1+{\eta }_n^s-{\zeta }_n^s\right)H_n\left(t\right)-\frac{{\xi }_{n+1}}{{\xi }_n}{\zeta }_n^{s}t{H}_n\left(t\right),$ (3.23)

同理比较(3.23)和(2.7)得到系数${\eta }_n^s$ 和${\zeta }_n^s$ 的转换公式

${\eta }_n^s=\frac{{\xi }_{n-1}}{{\xi }_{n+1}}{\eta }_n,$ (3.24)

${\zeta }_n^s=\frac{{\xi }_n}{{\xi }_{n+1}}{\zeta }_n,$ (3.25)

其中辅助系数${\xi }_n$ 的三项递归关系为

${\xi }_{n+1}=\left(1+{\eta }_n+{\zeta }_n\right){\xi }_n-{\eta }_n{\xi }_{n-1}.$ (3.26)

考虑共线系数${\delta }_n$ 和${\xi }_n$ 都是三项递归关系，迭代时需要两个初值，因此我们需要计算${\delta }_0,{\xi }_0,{\delta }_1,{\xi }_1$ 。根据残差的共线性，结合上式，当$n\ge 1$ ，我们可以得到

$R_n={\xi }_n{\delta }_n{R}_n^s.$

当$n=0$ 时，$X_0^s=X_0=0$ ，那么$R_0^s=R_0=C$ ，所以${\xi }_0{\delta }_0=1$ 。

结合(3.1)式和(3.15)，可得

${\varphi }_1\left(t\right)={\delta }_1{\varphi }_1^s={\delta }_1\left(1-{\alpha }_0^s\left(t+1\right)\right){\varphi }_0^s\left(t+1\right),$

又因为${\varphi }_0\left(t\right)={\delta }_0{\varphi }_0^s\left(t+1\right)$ ，所以

${\varphi }_1\left(t\right)={\delta }_1{\varphi }_1^s={\delta }_1\left(1-{\alpha }_0^s\left(t+1\right)\right){\varphi }_0^s\left(t+1\right)=\frac{{\delta }_1}{{\delta }_0}\left(1-{\alpha }_0^s\left(t+1\right)\right){\varphi }_0\left(t\right),$

当$t=-1$ 时，${\varphi }_1\left(-1\right)=\frac{{\delta }_1}{{\delta }_0}{\varphi }_0\left(-1\right)$ ，对比(2.4)可以得到

${\delta }_1=\left(1+{\alpha }_0\right){\delta }_0,$

同理有

${\xi }_1=\left(1+{\zeta }_0\right){\xi }_0.$

现在，重新考虑位移全局GPBiCG递归关系中的$R_n^s$ 和$T_n^s$ ，并使用残差共线多项式

$R_n^s=\frac{1}{{\xi }_n{\delta }_n}{R}_n,T_n^s=\frac{1}{{\xi }_n{\delta }_{n+1}}{T}_n,$ (3.27)

来降低它们的运行成本，得到

$\left({\rm M}+I\right)T_n^s=\left({\rm M}+I\right)\frac{1}{{\xi }_n{\delta }_{n+1}}{T}_n=\frac{1}{{\xi }_n{\delta }_{n+1}}{\rm M}{T}_n+\frac{1}{{\xi }_n{\delta }_{n+1}}{T}_n.$ (3.28)

因此，在实际计算时，只需要计算种子系统中的$M{T}_n^s$ ，而不需要计算位移系统中的$\left({\rm M}+I\right)T_n^s$ ，减少了算法的计算量。

在上述推导的基础上，我们给出求解位移系统的一种位移全局GPBiCG算法见表2。

Table 2. The shifted global GPBiCG method for solving shifted system (1.4)

表2. 位移全局GPBiCG算法求解位移系统(1.4)

算法3.1：位移全局GPBiCG算法求解位移系统(1.4)

给定初值$X_0=X_0^s=0_{n\times s}$ ，计算$R_0=C-\left({\rm M}+I\right)\left(X_0\right),$ 选择$R_0^{*}=R_0$ ，满足$\langle R_0^{*},R_0\rangle \ne 0,$ 令$T_{-1}=W_{-1}=0_{N\times S},{\beta }_{-1}=0,{\delta }_0={\xi }_0={\delta }_1={\xi }_1=1$ 对于$n=0,1,\cdots$ 计算 $P_n=R_n+{\beta }_{n-1}\left(P_{n-1}-U_{n-1}\right),$

${\alpha }_n=\frac{{\langle R_0^{*},R_n\rangle }_F}{{\langle R_0^{*},{\rm M}{P}_n\rangle }_F},$ $T_n=R_n-{\alpha }_n{\rm M}{P}_n,$ $Y_n=T_{n-1}-R_n-{\alpha }_n{W}_{n-1}+{\alpha }_n{\rm M}{P}_n,$ ${\zeta }_n==\frac{{\langle Y_n,Y_n\rangle }_F{\langle {\rm M}{T}_n,T_n\rangle }_F-{\langle Y_n,T_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F}{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,Y_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F},$ ${\eta }_n=\frac{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,T_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,T_n\rangle }_F}{{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle }_F{\langle Y_n,Y_n\rangle }_F-{\langle Y_n,{\rm M}{T}_n\rangle }_F{\langle {\rm M}{T}_n,Y_n\rangle }_F},$ 如果$n=0$ ，则${\zeta }_n=\frac{\langle {\rm M}{T}_n,T_n\rangle }{\langle {\rm M}{T}_n,{\rm M}{T}_n\rangle },{\eta }_n=0.$ $U_n={\zeta }_n{\rm M}{P}_n+{\eta }_n\left(T_{n-1}-R_n+{\beta }_{n-1}{U}_{n-1}\right),$ $Z_n={\zeta }_n{R}_n+{\eta }_n{Z}_{n-1}-{\alpha }_n{U}_n,$ $R_{n+1}=T_n-{\eta }_n{Y}_n-{\zeta }_n{\rm M}{T}_n,$ ${\beta }_n=\frac{{\alpha }_n}{{\zeta }_n}\cdot \frac{{\langle R_0^{*},R_{n+1}\rangle }_F}{{\langle R_0^{*},R_n\rangle }_F},$ $X_{n+1}=X_n+{\alpha }_n{P}_n+Z_n,$ $W_n={\rm M}{T}_n+{\beta }_n{\rm M}{P}_n,$ ${\delta }_{n+1}=\left(1+{\alpha }_n\right){\delta }_n+\frac{{\alpha }_n{\beta }_{n-1}}{{\alpha }_{n-1}}\left({\delta }_n-{\delta }_{n-1}\right),$ ${\xi }_{n+1}=\left(1+{\eta }_n+{\zeta }_n\right){\xi }_n-{\eta }_n{\xi }_{n-1},$ $P_n^s=R_n^s+{\beta }_{n-1}^s\left(P_{n-1}^s-U_{n-1}^s\right),$ ${\alpha }_n^s=\frac{{\delta }_n}{{\delta }_{n+1}}{\alpha }_n,$ ${\beta }_n^s={\left(\frac{{\delta }_n}{{\delta }_{n+1}}\right)}^2{\beta }_n,$ $T_n^s=\frac{1}{{\xi }_n{\delta }_{n+1}}{T}_n,$ ${\eta }_n^s=\frac{{\xi }_{n-1}}{{\xi }_{n+1}}{\eta }_n,$ ${\zeta }_n^s=\frac{{\xi }_n}{{\xi }_{n+1}}{\zeta }_n,$ $U_n^s={\zeta }_n^s\left({\rm M}+I\right)P_n^s+{\eta }_n^s\left(T_{n-1}^s-R_n^s+{\beta }_{n-1}^s{U}_{n-1}^s\right),$ $Z_n^s={\zeta }_n^s{R}_n^s+{\eta }_n^s{Z}_{n-1}^s-{\alpha }_n^s{U}_n^s,$ $X_{n+1}^s=X_n^s+{\alpha }_n^s{P}_n^s+Z_n^s,$ $R_{n+1}^s=\frac{1}{{\xi }_{n+1}{\delta }_{n+1}}{R}_{n+1},$

结束.

□

4. 数值算例

在这节，我们给出两个数值实验来说明位移全局GPBiCG方法(SGl-GPBiCG)的有效性。我们从迭代次数(Its)和计算时间(CPU，单位秒)两方面比较位移型方法SGl-GPBiCG与非位移型方法Gl-GPBiCG。在所有算例中，初始近似解都为零矩阵，最大迭代次数为2500次，停止迭代的条件为${\Vert R_{n+1}\Vert }_F/{\Vert R_0\Vert }_F<{10}^{-6}$ 。

例4.1 给定如下Stein矩阵方程

$X+AXB=C,$

其中

$A={\left(\begin{array}{ccccc}{D}_A& 0& 0& \cdots & 0\\ I_u& D_A& 0& \ddots & 0\\ 0& I_u& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & D_A& 0\\ 0& 0& \cdots & I_u& D_A\end{array}\right)}_{n\times n}\in R^{u^2\times u^2}$ ，$D_A={\left(\begin{array}{ccccc}14& 7& 0& \cdots & 0\\ 1& 14& 7& \cdots & 0\\ 0& 1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 14& 7\\ 0& 0& \cdots & 1& 14\end{array}\right)}_{u\times u}\in R^{u\times u}$ ，$B={\left(\begin{array}{ccccc}8& 3& 0& \cdots & 0\\ 3& 8& 3& \cdots & 0\\ 0& 3& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & 8& 3\\ 0& 0& \cdots & 3& 8\end{array}\right)}_{s\times s}\in R^{s\times s}$ ，$C=eye\left(n,s\right)$ 。

数值结果见表3和图1。从表3可以看出SGl-GPBiCG方法比Gl-GPBiCG的迭代次数更少。此外，从图1也可以看出，SGl-GPBiCG方法的收敛速度优于Gl-GPBiCG方法。

Table 3. Numerical results for Example 4.1

表3. 例4.1的数值结果

Matrix order	Methods	Its	CPU
$n=u^2=225,s=200$	Gl-GPBICG	38	2.0023
	SGl-GPBICG	21	3.0354
$n=u^2=625,s=200$	Gl-GPBICG	47	10.0729
	SGl-GPBICG	25	11.2680
$n=u^2=1296,s=200$	Gl-GPBICG	52	25.0586
	SGl-GPBICG	31	37.8359

(a) $n=u^2=225,s=200$

(b) $n=u^2=625,s=200$

(c) $n=u^2=1296,s=200$

Figure 1. Convergence curves for Example 4.1

图1. 例4.1的收敛曲线

例4.2 给定Stein矩阵方程为

$X+AXB=C,$

其中

$A=\text{triu}\left(rand\left(n\right),1\right)+diag\left(5+diag\left(rand\left(n\right)\right)\right),$

$B=\text{tril}\left(rand\left(n\right),1\right)+diag\left(8+diag\left(rand\left(n\right)\right)\right),C=rand\left(n\right).$

数值结果见表4和图2。从表4可以看出，SGl-GPBiCG方法比Gl-GPBiCG方法的迭代次数更少，消耗的CPU时间也更少。另外，从图2也可以看出，SGl-GPBiCG方法的收敛速度明显优于Gl-GPBiCG方法，并且随着矩阵阶数$n$ 的增加，收敛效果更明显。

Table 4. Numerical results for Example 4.2

表4. 例4.2的数值结果

Matrix order	Methods	Its	CPU
n = 100	Gl-GPBICG	25	0.6745
	SGl-GPBICG	21	0.4743
n = 200	Gl-GPBICG	45	4.3436
	SGl-GPBICG	33	3.6176
n = 300	Gl-GPBICG	75	20.7008
	SGl-GPBICG	53	12.9074
n = 500	Gl-GPBICG	178	77.1783
	SGl-GPBICG	50	12.9074

(a) $n=100$

(b) $n=200$

(c) $n=300$

(d) $n=500$

Figure 2. Convergence curves for Example 4.2

图2. 例4.2的收敛曲线

5. 总结

在本文中，我们提出一种位移全局GPBiCG方法求解大型Stein矩阵方程。该方法充分利用了Stein矩阵方程的特殊位移结构。同时，该方法去除了系数矩阵$A$ 和$B$ 间的维数限制。最后，实验结果验证了位移全局GPBiCG方法求解Stein矩阵方程的有效性，尤其在矩阵$A$ 和$B$ 均为大规模矩阵的情况下优势明显。

NOTES

^*通讯作者。