Sylvester矩阵方程在很多领域中都有涉及, 例如矩形域上的椭圆边值问题离散后会得到一个Sylvester形式的矩形方程；在常微分方程定性理论研究及数值求解的隐式Runge-Kwutta方法与块方法中常会碰到这类方程的求解问题；在线性统计领域中也会有Sylvester矩阵方程；在控制理论及应用中, 极点配置、观测器及构造Sylvester函数中都涉及上述方程的求解.除此之外, 特殊形式的Sylvester方程, 即Lyapunov方程, 在实际应用中, 特别是在连续和离散稳定性分析中具有重要意义, 在系统与控制领域中经常会遇到.与此同时, Sylvester矩阵方程在图像恢复、模型降阶、神经网络、信号处理等领域中也有涉及.Sylvester矩阵方程的求解问题通常包括精确解和数值解的研究.尽管在系统稳定性分析等众多应用中, 精确解的研究非常重要, 但是, 当系统的参数未知时, 不可能得到精确解, 此时数值解的研究就非常必要.近些年, 国内外很多学者对各种类型的Sylvester矩阵方程的数值解进行研究, 得到了一些有效的数值求解算法^[1-15].

本文从Sylvester矩阵方程AXB=F的迭代解法入手, 先给出了它的迭代解的收敛条件, 然后介绍了计算Sylvester矩阵方程AXB+CXD=F数值解的分层迭代法，以及解决Sylvester矩阵方程AX+XB=C的松弛迭代法.

1 矩阵方程AXB+CXD=F的松弛迭代法 1.1 相关引理

对于任意的两个整数i < j, I[i, j]表示集合{i, i+1, …, j}, 对任意方阵A, , tr (A), ‖A‖分别表示矩阵共轭, 矩阵的迹和矩阵的F范数.对X=[x₁x₂…x_n]∈C^m×n, 设vec (X)=[x₁^Tx₂^T…x_n^T]^T, ‖A‖=‖vec (A)‖.对于任意矩阵M和N, M∈C^r×m, N∈C^n×s, MXN表示克罗内克积, 有vec(MXN)=(N^T⊗M)vec (X).

引理1^[2]   对于AXB=F, 其中A∈C^p×m, B∈C^n×q, F∈C^p×q是已知矩阵, X∈C^m×n是未知矩阵.若满足或, 且求解方程AXB=F的迭代算法表示为X(k+1)=X(k)+μA^T(F-AX(k)B)B^T.若方程AXB=F有唯一解X_*, 则迭代解X(k)收敛到唯一解, 即.

引理2^[2]   若A∈C^p×m是列满秩矩阵, B∈C^n×q是行满秩矩阵(p≥m, n≤q), 则方程AXB=F有唯一解, 且X=(A^TA)^-1A^TFB^T(BB^T)^-1, 所对应的齐次方程AXB=O有唯一解X=O.

1.2 迭代格式的建立

考虑Sylvester矩阵方程

(1)

其中:A, C∈C^m×r；B, D∈C^s×n；F∈C^m×n是已知的常数矩阵；X∈C^r×s是待求矩阵.

假设X_*表示方程(1) 的精确解, X(k)表示X_*的第k步迭代的近似解.应用分层辨识原理^{[2-5, 7-8]}, 定义矩阵

(2)

(3)

引入一个松弛因子, 建立如下松弛递归格式:

(4)

(5)

其中:μ是收敛因子；ω是一个松弛因子, 满足0 < ω < 1.

将方程(2) 和(3) 代入到方程(4) 和(5) 中, 得到如下格式:

(6)

(7)

由于上述两个等式的右端包括未知矩阵X, 算法无法实现.根据分层辨识原理^{[2-5, 7-8]}, 方程(6) 和(7) 中X可以被它的迭代近似值X(k-1) 代替, 得到

(8)

(9)

由于求的是近似解X(k)而不是X₁(k)和X₂(k), 遵循平衡策略形成第k步迭代的近似解:

(10)

(11)

(12)

于是迭代算法改成如下格式:

(13)

1.3 收敛性分析

定理1   设方程AXB+CXD=F有唯一解X_*, 其中A, C∈C^m×r, B, D∈C^s×n, F∈C^m×n, X∈C^r×s.设λ_max(AA^T)=λ₁, λ_max(BB^T)=λ₂, λ_max(CC^T)=λ₃, λ_max(DD^T)=λ₄.若收敛因子μ满足, 则对于任意初值X(0), 算法(13) 给出的迭代解收敛于唯一解X_*, 即.

证明  定义误差矩阵

(14)

(15)

(16)

设

(17)

将式(15) 代入到式(10), (11)：

将式(16), (17) 代入到上式, 可得X₁(k)=X(k-1)-ωμA^T(ξ(k-1)+η(k-1))B^T,

X₂(k)=X(k-1)-(1-ω)μC^T(ξ(k-1)+η(k-1))D^T.

由于|X|²=tr (X^TX), tr (AB)=tr (BA), tr (A)=tr (A^T), 得

.设λ_max(AA^T)=λ₁, λ_max(BB^T)=λ₂, λ_max(CC^T)=λ₃, λ_max(DD^T)=λ₄.可得‖A^T(ξ(k-1)+η(k-1))B^T‖²≤λ₁λ₂‖(k-1)+η(k-1)², 则 tr [(ξ(k-1)+η(k-1))ξ^T(k-1)+(ξ(k-1)+η(k-1))^Tξ(k-1)].则+ω²μ²λ₁λ₂‖ξ(k-1)+η(k-1)‖²-μωtr [(ξ(k-1)+η(k-1))ξ^T(k-1)+(ξ(k-1)+η(k-1))^Tξ(k-1)].

同样(1-ω)²μ²λ₃λ₄‖ξ(k-1)+η(k-1)‖²-μ(1-ω)tr [(ξ(k-1)+η(k-1))η^T(k-1)+(ξ(k-1)+η(k-1))^Tη(k-1)].

将等式(12) 两边同时取范数, 有

可得.若μ满足, 有, 则‖ξ(k)+η(k)‖²→0.进而得ξ(k)+η(k)→0, 结合引理1, .

2 数值模拟

例  考虑方程AXB+CXD=F, 其中: 设方程AXB+CXD=F的精确解为

利用算法(10)~(12), 取初始矩阵, 定义相对误差.根据定理1计算得到, λ_max(AA^T)=λ₁=2, λ_max(BB^T)=λ₂=2, λ_max(CC^T)=λ₃=5, λ_max(DD^T)=λ₄=2, 故有

通过大量实验, 当ω=0.49及μ==0.283 29时, 有更小的相对误差和较快的迭代速度, 如图 1~图 3所示.

当ω=0.49, μ=100/353时, 表 1给出用算法(13) 和文献[2]中算法求得解与相对误差δ₁.

通过图表可以发现, 算法(13) 明显优于文献[2]中的算法.

3 结论

本文讨论了实数域下Sylvester矩阵方程的迭代解法, 其主要思想是基于梯度迭代法与分层迭代原理, 通过引入松弛因子, 建立新的迭代格式, 从而获得Sylvester矩阵方程的数值解.还给出了相应的收敛性分析, 进行了理论证明, 得出如下结论:在特定条件下, 对于任意初始值, 方程的迭代解均收敛于精确解.最后通过数值算例说明了所给新算法的有效性、可行性以及优越性.从讨论可以看出, 同等情况下, 带有松弛因子的迭代算法与之前的方法相比, 可以带来较小的迭代误差与较快的迭代速度.另外, 本文给出的迭代算法更具有普遍性, 可以解多种Sylvester矩阵方程.