基于SOR预条件的全局GMRES方法求解Sylvester矩阵方程

SOR Preconditioned Global GMRES Method for Solving Sylvester Matrix Equations

摘要: Sylvester矩阵方程在控制理论、图像处理及微分方程数值解等领域应用广泛。对于大规模问题，直接求解方法因计算量和存储需求巨大而不再适用。全局GMRES (Generalized Minimal Residual Method)方法是求解此类大规模方程的有效迭代方法。为提升其收敛速度，本文引入了基于SOR (Successive Over-Relaxation)迭代的预条件技术，构建了SOR预条件全局GMRES方法。数值实验表明，与传统的全局GMRES方法相比，新方法能显著减少迭代步数与计算时间，验证了其有效性。

Abstract: The Sylvester matrix equation is widely used in fields such as control theory, image processing, and numerical solutions to differential equations. For large-scale problems, direct solution methods become impractical due to their high computational cost and storage requirements. The global Generalized Minimal Residual Method (GMRES) is an effective iterative approach for solving such large-scale equations. To enhance its convergence speed, this paper introduces a preconditioning technique based on the Successive Over-Relaxation (SOR) iteration, establishing the SOR-preconditioned global GMRES method. Numerical experiments demonstrate that, compared to the traditional global GMRES method, the proposed approach significantly reduces the number of iteration steps and computational time, confirming its effectiveness.

1. 引言

Sylvester矩阵方程是科学计算中的一类重要方程，其形式为：

$AX+XB=C$

其中$A\in {ℝ}^{m\times m},B\in {ℝ}^{n\times n},C\in {ℝ}^{m\times n}$ 为已知矩阵，$X\in {ℝ}^{m\times n}$ 为待求的未知矩阵。该方程在控制和通信理论、模型降阶、图像恢复等诸多领域都有着重要的应用[1] [2]。

求解Sylvester方程的方法主要分为直接法和迭代法两大类。直接法(如Bartels-Stewart算法[3])适合求解中小规模问题。对于由偏微分方程离散化等产生的大规模稀疏问题，直接法往往因计算复杂度和存储需求过高而变得不再适用。因此，发展高效、低存储的迭代法成为研究的重点。

在众多迭代法中，基于Krylov子空间的方法，如GMRES [4]，因其适用于非对称矩阵而备受关注。为了更高效地处理矩阵形式的未知量，有学者提出了全局Krylov子空间方法的概念[5]。这类方法将整个矩阵视为迭代变量，在扩展的Krylov子空间中寻找解，通常比将其向量化的方法具有更高的计算效率和更低的存储需求。

尽管全局GMRES方法在处理大规模问题上较为有效，但其收敛速度严重依赖于系数矩阵A和B的谱性质。当矩阵谱分布不佳时，收敛可能非常缓慢。预条件技术是改善迭代法收敛性的关键手段[6]。其核心思想是将原系统转化为一个等价的、具有更优谱分布的系统，从而加速迭代收敛。

本文旨在研究如何将SOR预条件技术与全局GMRES方法结合，以高效求解Sylvester方程。本文结构如下：第二节介绍预备知识和全局GMRES方法；第三节讨论SOR预条件技术及其实现；第四节通过数值实验验证算法性能；第五节总结全文。

2. 全局GMRES方法

考虑Sylvester方程$AX+XB=C$ 。利用Kronecker积和向量化算子$vec(\cdot )$ ，可将其等价地转化为线性方程组：

$Ax=c$

其中$A=I_n\otimes A+B^{\text{T}}\otimes I_m,x=vec\left(X\right),c=vec\left(C\right)$ 。此时如果应用广义最小残差法(GMRES)等Krylov子空间方法求解会发现：当$m,n$ 很大时，$mn\times mn$ 的矩阵$A$ 将导致巨大的存储和计算成本。因此需要采用全局迭代法来避免这种显式的向量化[5]。

定义线性算子：$L\left(X\right)=AX+XB$ ，则原方程等价于$L\left(X\right)=C$ 。与传统Krylov子空间针对向量不同，全局Krylov子空间是针对矩阵定义的。对于初始残差矩阵$R_0=C-L\left(X_0\right)$ (其中$X_0$ 为初始猜测，通常取零矩阵)，其对应的$k$ 维全局Krylov子空间定义为：

$K_k\left(L,R_0\right)=span\left\{R_0,L\left(R_0\right),L^2\left(R_0\right),\cdots ,L^{k-1}\left(R_0\right)\right\}$

全局GMRES方法的目标是在该矩阵Krylov子空间中寻找近似解。通过矩阵的Frobenius内积${\langle U,V\rangle }_F=trace\left(U^{\text{T}}V\right)$ 及其诱导的范数${\Vert U\Vert }_F=\sqrt{{\langle U,U\rangle }_F}$ ，生成一组F-正交的矩阵序列$\left\{V_1,V_2,\cdots ,V_k\right\}$ ，满足：

${\langle V_i,V_j\rangle }_F=trace\left(V_i^{\text{T}}{V}_j\right)={\delta }_{ij}=\{\begin{array}{cc}1& i=j\\ 0& i\ne j\end{array}$

全局Arnoldi过程的算法[5]如下：

算法1全局Arnoldi过程

1. 令 $V_1=R_0/{\Vert R_0\Vert }_F$ .

2. For$j=1,2,\cdots ,k$

$W=L\left(V_j\right)=A{V}_j+V_{j}B$

For $i=1,2,\cdots ,j$

$h_{ij}={\langle W,V_j\rangle }_F$

$W=W-h_{ij}{V}_i$

End For

$h_{j+1,j}={\Vert W\Vert }_F$ ，若$h_{j+1,j}=0$ 则停止，否则$V_{j+1}=W/h_{j+1,j}$

End For

此过程生成了一个F-正交基${\text{V}}_k=\left[V_1,V_2,\cdots ,V_k\right]$ 和一个上Hessenberg矩阵${\overline{H}}_k\in {ℝ}^{\left(k+1\right)\times k}$ 。它们满足全局Arnoldi关系：

$L\left({\text{V}}_k\right)={\text{V}}_{k+1}{\overline{H}}_k$

全局GMRES方法的目标是在k维全局Krylov子空间$K_k$ 中寻找近似解$X_k=X_0+{\displaystyle \sum _{i=1}^k{y}_i{V}_i}$ 。令向量$y={\left[y_1,y_2,\cdots ,y_k\right]}^{\text{T}}\in {ℝ}^k$ ，该近似解可以通过最小化残差范数来确定：

$\underset{y\in {ℝ}^k}{\min }{\Vert C-L\left(X_k\right)\Vert }_F=\underset{y\in {ℝ}^k}{\min }{\Vert C-L\left(X_0+{\displaystyle \sum _{i=1}^k{y}_i{V}_i}\right)\Vert }_F$

利用全局Arnoldi关系，上式可转化为：

$\underset{y\in {ℝ}^k}{\min }{\Vert {\Vert R_0\Vert }_F{e}_1-{\overline{H}}_{k}y\Vert }_2$

其中$e_1={\left[1,0,\cdots ,0\right]}^{\text{T}}\in {ℝ}^{k+1}$ 。这是一个小型的最小二乘问题，可以通过Givens旋转高效求解得到向量y，从而得到近似解

$X_k=X_0+{\displaystyle \sum _{i=1}^k{y}_i{V}_i}$

3. SOR预条件全局GMRES方法

3.1. SOR预条件子

全局GMRES方法是求解Sylvester方程的有效方法。但是对于病态或困难问题，全局GMRES的收敛速度可能会较为缓慢。因此，引入高效的预条件技术至关重要。

预条件的基本思想是将原系统$Ax=c$ 转化为一个具有更优谱性质的新系统，从而加速收敛。对于Sylvester方程，我们寻求一个预条件子$M$ ，使得$M^{-1}Ax=M^{-1}c$ 更容易求解。

逐次超松弛(SOR: Successive Over-Relaxation)方法作为一种经典的分裂迭代法，其迭代矩阵可以自然地导出一种分裂预条件子。SOR预条件子具有形式简单、易于实现等优点。我们首先回顾对于普通线性系统$Ax=b$ 的SOR迭代过程[7]：

设$A=D-L-U$ ，其中$D,-L,-U$ 分别为矩阵A的对角、严格下三角和严格上三角部分。SOR迭代格式为：

$x^{\left(k+1\right)}={\left(D-\omega L\right)}^{-1}\left[\left(1-\omega \right)D+\omega U\right]x^{\left(k\right)}+\omega {\left(D-\omega L\right)}^{-1}b$

其中$\omega \in \left(0,2\right)$ 为松弛因子。相应的SOR预条件子为

$M_{SOR}=\frac{1}{\omega }\left(D-\omega L\right)$

我们推广到Sylvester方程$AX+XB=C$ 。SOR预条件子$M_{SOR}$ 定义为使得一次SOR迭代等价于求解$M_{SOR}^{-1}\left(AX+XB\right)=M_{SOR}^{-1}\left(C\right)$ 的算子。

我们对矩阵A和B进行如下分裂：$A=D_A-L_A-U_A,B=D_B-L_B-U_B$ ，其中$D_A,D_B,-L_A,-L_B,-U_A,-U_B$ 分别是矩阵A和B的对角、严格下三角和严格上三角部分。

此时$AX+XB=C$ 对应的线性方程组$Ax=c$ 在 Kronecker 乘积形式下可以分裂为：

$A=\left(I\otimes D_A+D_B^{\text{T}}\otimes I\right)-\left(I\otimes L_A+U_B^{\text{T}}\otimes I\right)-\left(I\otimes U_A+L_B^{\text{T}}\otimes I\right)$

因此相应的SOR预条件子为

$M_{SOR}=\frac{1}{\omega }\left(I\otimes \left(D_A-\omega L_A\right)+{\left(D_B-\omega U_B\right)}^{\text{T}}\otimes I\right)$

从结构上可以看出，$M_{SOR}$ 主要由系数矩阵的对角矩阵和下三角矩阵构成，通过引入松弛因子$\omega$ 进行缩放，可以调整预条件算子的性质。当矩阵的对角占优特性不明显时，适当选择松弛因子$\omega$ 可以增强对角部分的作用，改善矩阵的条件数，从而加速迭代收敛。

而且由于$M_{SOR}$ 的构成核心是两个下三角矩阵$\left(D_A-\omega L_A\right)$ 和${\left(D_B-\omega U_B\right)}^{\text{T}}$ ，因此其对应的线性系统可以高效求解。具体来说，作用于残差$R=C-AX-XB$ 的SOR预条件步$Z=M_{SOR}^{-1}\left(R\right)$ 可以通过求解以下方程来实现：

$\left(D_A-\omega L_A\right)Z+Z\left(D_B-\omega U_B\right)=\omega R$

因为该Sylvester方程的系数矩阵均为三角矩阵，因此可以利用逐列或逐行的方式快速求解，计算成本远低于直接处理原方程。我们给出求解系数矩阵为上、下三角矩阵的Sylvester方程的快速算法如下：

算法2 求解Sylvester方程$LX+XU=C$ 算法

输入：$L\in {ℝ}^{m\times m}$ 是下三角矩阵，$U\in {ℝ}^{n\times n}$ 是上三角矩阵，$C\in {ℝ}^{m\times n}$

输出：$X\in {ℝ}^{m\times n}$

For $j=1,2,\cdots ,n$

For $i=1,2,\cdots ,m$

$X_{ij}=\frac{1}{L_{ii}+U_{jj}}\left(C_{ij}-L\left(i,1:i-1\right)X\left(1:i-1,j\right)-X\left(i,1:j-1\right)U\left(1:j-1,j\right)\right)$

End For

End For

基于上述SOR预条件技术，我们构建出求解Sylvester方程的SOR预条件全局GMRES算法：

算法3 SOR预条件全局GMRES算法(SOR-P-GL-GMRES)

1．初始化：选择初始解$X_0$ ，计算残差$R_0=C-A{X}_0-X_{0}B$

2．$V_1=M_{SOR}^{-1}\left(R_0\right)/{\Vert M_{SOR}^{-1}\left(R_0\right)\Vert }_F$

3．For $j=1,2,\cdots ,k$

$W=M_{SOR}^{-1}\left(A{V}_j+V_{j}B\right)$

For $i=1,2,\cdots ,j$

$h_{ij}={\langle W,V_j\rangle }_F$

$W=W-h_{ij}{V}_i$

End For

$h_{j+1,j}={\Vert W\Vert }_F$ ，若$h_{j+1,j}=0$ 则停止，否则$V_{j+1}=W/h_{j+1,j}$

End For

4．最小化问题：寻找y使得${\Vert {\Vert R_0\Vert }_F{e}_1-{\overline{H}}_{k}y\Vert }_2$ 最小

5．更新解：$X_k=X_0+{\displaystyle \sum _{i=1}^k{y}_i{V}_i}$

3.2. 松弛因子的选取策略

SOR预条件全局GMRES算法方法的一个难点在于如何确定最优松弛因子${\omega }_{opt}$ 。理论上，最优松弛因子与系数矩阵$A$ 和$B$ 的谱特性密切相关。对于对称正定矩阵等特殊矩阵，可以找出较为合适的松弛因子，但对于系数矩阵非对称的Sylvester方程，难以精确计算其最优松弛因子。此时可以采用数值搜索方法来选取松弛因子，即在一定范围内对松弛因子$\omega$ 进行搜索，以找到使迭代算法收敛速度最快的$\omega$ 值。

其中试算法是一种较为直接的数值搜索方法。其基本思路是在一个合理的范围内，循环尝试多个不同的$\omega$ 值，针对每个$\omega$ 值运行迭代算法，并记录其迭代次数或收敛时间，最后比较不同$\omega$ 值下的收敛速度，选择使收敛速度最快的$\omega$ 作为最优松弛因子。

如果系数矩阵$A$ 和$B$ 的维数很大，用试算法可能成本太高。此时可以考虑一种替代方案：从原问题中提取一个小模型，其矩阵$A_{small}$ 和$B_{small}$ 保持原矩阵的主要性质但规模小得多。在这个小问题上执行对$\omega$ 的搜索，然后将找到的最佳值${\omega }_{opt}$ 用于原始的大规模问题。通过这些策略，可以使得SOR预条件全局GMRES方法适应多样化的应用场景。

4. 数值实验

本节通过数值算例验证SOR预条件全局GMRES (SOR-P-GL-GMRES)方法的有效性。我们将其与未预条件的全局GMRES (GL-GMRES)方法进行比较。所有实验均在MATLAB R2024a环境中进行，机器配置为Macbook M4，16 GB RAM。在下面的例子中，初始解均取为$X_0=O$ 。停机准则是

$\frac{{\Vert R_k\Vert }_F}{{\Vert R_0\Vert }_F}\le {10}^{-11}$

例1 我们考虑二维对流–扩散方程$-\Delta u+v\cdot \nabla u=f$ 离散化后产生的Sylvester方程。令$A=tridiag\left(-1-\alpha ,4,-1+\alpha \right),B=tridiag\left(-1-\beta ,4,-1+\beta \right)$ 为三对角矩阵。右端项C随机生成。$A\in {ℝ}^{m\times m},B\in {ℝ}^{n\times n},C\in {ℝ}^{m\times n}$ 。取$m=160,n=180,\alpha =0.2,\beta =1.6$ ，松弛因子$\omega =1.1$ 。表1和图1展示了两种方法的收敛性能。从中可以看出，SOR预条件技术显著加速了收敛。SOR-P-GL-GMRES需要26步迭代，而GL-GMRES需要58步迭代。

Table 1. Performance comparison of global GMRES and SOR-preconditioned global GMRES (Example 1)

表1. 全局GMRES与SOR预条件全局GMRES的性能比较(例1)

方法	迭代步数	耗时(秒)	相对残差
GL-GMRES	58	0.35	8.92 × 10^−^12
SOR-P-GL-GMRES	26	0.22	3.65 × 10^−^12

Figure 1. Plot of relative residual norm versus iteration steps (Example 1)

图1. 相对残差范数随迭代步数下降图(例1)

例2 这个例子的结构和上个例子相同，我们取$m=500,n=300,\alpha =0.1,\beta =1.2$ ，松弛因子$\omega =1.2$ 。

两种方法的收敛性能展示在表2和图2中。从结果可以清晰地看出SOR-P-GL-GMRES的性能也要优于GL-GMRES。SOR-P-GL-GMRES需要24步迭代，而GL-GMRES需要49步迭代。

Figure 2. Plot of relative residual norm versus iteration steps (Example 2)

图2. 相对残差范数随迭代步数下降图(例2)

Table 2. Performance comparison of global GMRES and SOR-preconditioned global GMRES (Example 2)

表2. 全局GMRES与SOR预条件全局GMRES的性能比较(例2)

方法	迭代步数	耗时(秒)	相对残差
GL-GMRES	49	5.26	7.57 × 10^−^12
SOR-P-GL-GMRES	24	3.56	3.89 × 10^−^12

5. 结论

本文给出了求解大规模Sylvester矩阵方程的SOR预条件全局GMRES方法。数值实验结果表明，SOR预条件技术能有效改善算法的收敛性，减少迭代步数和计算时间。未来的工作将集中于研究更高效的预条件子，例如基于低秩近似或多重网格方法的预条件子，以应对更复杂和病态的问题。

基金项目

由2025年上海市级大学生创新训练计划项目S202510251170资助。