基于复西尔维斯特矩阵方程的改进双阶尺度分裂数值计算方法
Improved Numerical Calculation Method for Second-Order Scale Splitting Based on the Complex Sylvester Matrix Equation
DOI: 10.12677/airr.2024.132033, PDF, HTML, XML, 下载: 33  浏览: 84 
作者: 杨创勋:广东工业大学数学与统计学院,广东 广州
关键词: 西尔维斯特矩阵方程MDSS迭代法最优收敛因子Sylvester Matrix Equation MDSS Iterative Method Optimal Parameters
摘要: 连续时间的矩阵方程是矩阵方程中十分重要的一个类型,在矩形域上的椭圆边值问题的数值解法、线性统计、振动结构的共振控制、二次矩阵的特征值配置问题、受噪声影响的图像复原等问题中有重要地位。由于连续时间的矩阵方程广泛的应用背景,因此,对连续时间的Sylvester方程的数值解法的研究具有重要的理论和实际意义.改进的双步尺度分裂(MDSS)方法是研究一类大型复杂对称线性系统的一种有效方法。在本文中,我们将在双步尺度分裂(DSS)方法的基础上通过改进,得到MDSS方法,来求解复数域上的连续时间的西尔维斯特矩阵方程的近似解,证明了该迭代序列在任意初始条件的情况下都收敛于西尔维斯特矩阵方程的唯一解,并确定了其最优参数和相应的最优收敛因子。最后,给出了一个测试问题来说明该新技术的有效性。
Abstract: The continuous time matrix equation is a very important type of matrix equation, which plays an important role in numerical solutions of elliptic boundary value problems in rectangular domains, linear statistics, resonance control of vibration structures, eigenvalue assignment of quadratic matrices, image restoration affected by noise, and other problems. Due to the wide application background of continuous time matrix equations, the study of numerical solutions for the continuous time Sylvester equation has important theoretical and practical significance. The improved two-step scale splitting (MDSS) method is an effective method for studying a class of large complex symmetric linear systems. In this article, we will improve the two-step scale splitting (DSS) method and obtain the MDSS method to solve the approximate solution of the continuous time Sylvester matrix equation in the complex field. We prove that the iterative sequence converges to a unique solution of the Sylvester matrix equation under any initial condition, and determine its optimal parameters and corresponding optimal convergence factors. Finally, a testing question was presented to demonstrate the effectiveness of the new technology.
文章引用:杨创勋. 基于复西尔维斯特矩阵方程的改进双阶尺度分裂数值计算方法[J]. 人工智能与机器人研究, 2024, 13(2): 313-321. https://doi.org/10.12677/airr.2024.132033

1. 引言

在这里,我们重点关注西尔维斯特矩阵方程:

A X + X B = C , A m × m , B n × n , C m × n (1)

其中必须确定矩阵 X m × n 。式(1)在控制理论和许多其他工程分支中具有非常重要的应用。矩阵方程在科学计算和工程技术等诸多领域都有着广泛的应用,如控制论 [1] 、系统理论、信号处理 [2] 、图像恢复 [3] 、极点配置 [4] 、扰动分析、模型降阶 [5] 、神经网络 [6] 、矩阵逼近问题 [7] 、电力系统 [8] 、微分方程的数值解 [9] 等。因此,对连续时间的Sylvester方程数值解法的研究具有重要的理论和实际意义。对于Sylvester矩阵方程,传统的求解方法是将这个矩阵方程用Kronecker积转化为形如的等价线性方程组,然后再用求解线性方程组的方法去计算其数值解。然而在实际问题中,遇到的系数矩阵与规模一般都是较大的,用Kronecker积转化后,会使阶数也变得非常大,在很大程度上增加了工作量及计算机的储存空间,而要得到其解析式也十分困难,此时这种传统的解决方法变得不可行.近年来,国内外许多学者致力于研究线性方程组的各种求解方法,并推广到矩阵方程的求解,这些方法可以概括为直接法和迭代法两大类。直接法是一种精确的求解方法,常见的直接法有Bartels-stewart法 [10] 和Hessenberg-Schur法 [11] ,其原理是利用正交相似变换将矩阵转化为三角矩阵或者上Hessenberg矩阵,然后再利用回代法求解线性方程。当求解一些小型线性矩阵方程时,我们通常采用直接法。但对于大型稀疏Sylvester矩阵方程,由于系数矩阵的规模非常大,直接法在求解方程时可能会产生许多非零元,将消耗大量的存储空间和计算时间,所以对于大型稀疏矩阵,我们一般采用迭代法来求解。古典迭代法是基于矩阵分裂的迭代算法,其基本思想是将线性方程组问题转化为求解一个不动点方程的问题,进而构造迭代格式。因此,分裂法是常见的有效迭代法,包括Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法、交替方向法 [12] 、HSS迭代法 [13] 和Smith方法 [14] 等。在本文中,我们将在双步尺度分裂(DSS)方法的基础上通过改进,得到MDSS方法,来求解复数域上的连续时间的西尔维斯特矩阵方程的近似解,证明了该迭代序列在任意初始条件的情况下都收敛于西尔维斯特矩阵方程的唯一解,并确定了其最优参数和相应的最优收敛因子。最后,给出了一个测试问题来说明该新技术的有效性。

2. 复西尔维斯特矩阵方程的MDSS算法的研究

在本节中,先给出本文中使用的一些基本定义和符号。

符号 是所有复数的集合,维数为n的复向量的集合用 n 表示。符号 是所有实数的集合,并用 n 表示维数为n的实向量的集合。具有m行和n列的复数(或实数)矩阵的集合用 m × n (或 m × n )表示。用 表示克罗内克积,用I表示单位矩阵,用0表示零矩阵为0。此外,AT是A的转置,以及我们得到 A H = A ¯ T ,其中 A ¯ 表示矩阵A的共轭。

对于向量x,我们定义 x 2 为向量x的2范数,以及我们定义 . F 向量x的Frobenius范数,其中对于任意矩阵A,我们定义 A F = t r ( A A H ) 。在这期间,我们定义 t r ( A ) 为矩阵A的迹。同样的,我们定义Sp(A)为矩阵A的谱,则矩阵A的谱半径定义为 ρ ( A ) 。我们使得矩阵 A = ( a 1 , a 2 , , a n ) m × n ,其中 a k 是矩阵A的第k列。接着,我们可以定义 v e c ( A ) = ( a 1 , a 2 , , a n ) T ,其中 v e c ( A ) m × n 中的向量。

在研究本文提出的新算法前,我们将列出其他求解复西尔维斯特矩阵方程的算法,我们计划简要描述一下求解矩阵方程(1)的MHSS、PMHSS和GMHSS迭代方法 [15] [16] 。

2.1. 解决矩阵方程(1)的MHSS算法

通过计算求解 X ( k + 1 ) m × n ,再使用下列式子,使得 { X ( k ) } k = 0 收敛。

{ ( α I + W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( β I + U ) = ( α I i T ) X ( k ) + X ( k ) ( β I i V ) + C ( α I + T ) X ( k + 1 ) + X ( k + 1 ) ( β I + V ) = ( α I + i W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( β I + i U ) i C (2)

其中, α > 0 β > 0 X ( 0 ) m × n 是一个初始矩阵,I是单位矩阵。因此我们得到矩阵

M ( γ ) = ( γ I + D ) 1 ( γ I + i H ) ( γ I + H ) 1 ( γ I i D )

为MHSS方法的迭代矩阵,其中为 H = I W + U I D = I T + V I γ = α + β 。此外,还证明了 ρ ( M ( γ ) ) σ ( γ ) < 1 ,其中

σ ( γ ) : = max γ 2 + ( λ i W + λ j U ) 2 γ + λ i W + λ j U

因此,MHSS迭代收敛于西尔维斯特矩阵等式(1)的精确解 X m × n ,其中得到最优参数 γ ˜ = λ min λ max ,其中 λ min = λ min W + λ min U λ max = λ max W + λ max U ,同时 λ min W = min S p ( W ) λ min U = min S p ( U ) λ max U = max S p ( U ) 以及 λ max W = max S p ( W ) 。此外还有

σ ( γ ˜ ) = κ ( H ) + 1 κ ( H ) + 1

其中

κ ( H ) = λ min λ max

2.2. 解决矩阵方程(1)的PMHSS算法

通过计算求解 X ( k + 1 ) m × n ,再使用下列式子,使得 { X ( k ) } k = 0 收敛。

{ ( γ P 1 + W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( γ P 2 + U ) = ( γ P 1 i T ) X ( k ) + X ( k ) ( γ P 2 i V ) + C ( γ P 1 + T ) X ( k + 1 ) + X ( k + 1 ) ( γ P 2 + V ) = ( γ P 1 + i W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( γ P 2 + i U ) i C (3)

其中 γ > 0 X ( 0 ) m × n 是一个初始矩阵, P 1 P 2 是对称正定矩阵。而且

P ( γ ) = ( γ P + D ) 1 ( γ P + i H ) ( γ P + H ) 1 ( γ P i D )

为PMHSS方法的迭代矩阵,其中 P = I P 1 + P 2 T I

2.3. 解决矩阵方程(1)的GMHSS算法

通过计算求解 X ( k + 1 ) m × n ,再使用下列式子,使得 { X ( k ) } k = 0 收敛。

{ ( α I + W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( β I + U ) = ( α I i T ) X ( k ) + X ( k ) ( β I i V ) + C ( τ I + T ) X ( k + 1 ) + X ( k + 1 ) ( μ I + V ) = ( τ I + i W ) X ( k + 1 2 ) + X ( k + 1 2 ) ( μ I + i U ) i C (4)

其中, α β τ μ 为实正常数。GMHSS方法的迭代矩阵为

M ( γ , θ ) = ( θ I + D ) 1 ( θ I + i H ) ( γ I + H ) 1 ( γ I i D )

因此,对于不同的参数 α β τ μ ,提出了不同的方法。

接着,便是本文提出算法的研究:

2.4. 解决矩阵方程(1)的MDSS算法

在本节中,我们尝试应用MDSS迭代方法来求解矩阵方程(1)。首先,设 A = W + i T B = U + i V ,其中WTUV是实对称正定矩阵,根据上述等式以及变换,我们可以将矩阵方程(1)写成

( W + i T ) X + X ( U + i V ) = C . (5)

现在假设 α β 是一个正的实常数。通过在等式(6)中分别乘以式子 ( α β i ) ( β α i ) ,我们得到

( α W + i α T i β W + β T ) X + X ( α U + i α V i β U + β V ) = ( α i β ) C (6)

以及式子

( β W + i β T i α W + α T ) X + X ( β U + i β V i α U + α V ) = ( β i α ) C (7)

在这里,我们使得 C 1 = ( α i β ) C C 2 = ( β i α ) C ,然后我们得到

( α W + β T ) X + X ( α U + β V ) = i ( β W α T ) X + i X ( β U α V ) + C 1 (8)

以及

( α T + β W ) X + X ( α V + β U ) = i ( α W β T ) X + i X ( α U β V ) + C 2 (9)

然后通过上述式子,我们得到了求解矩阵方程(1)的算法。

通过计算求解 X ( k + 1 ) m × n ,再使用下列式子,使得 { X ( k ) } k = 0 收敛:

{ ( α W + β T ) X ( k + 1 2 ) + X ( k + 1 2 ) ( α U + β V ) = i ( β W α T ) X ( k ) + i X ( k ) ( β U β V ) + C 1 ( α T + β W ) X ( k + 1 ) + X ( k + 1 ) ( α V + β U ) = i ( α W β T ) X ( k + 1 2 ) + i X ( k + 1 2 ) ( α U β V ) + C 2 (10)

其中, α > 0 β > 0 X ( 0 ) m × n 是一个初始矩阵,I是单位矩阵。很明显,矩阵 α W + β T α U + β V α T + β W α V + β U 是对称正定的。

对此,我们可以得到以下定理。

定理2.1 首先,我们定义 E = D + i H ,得到

D = I W + U I , H = I T + V I (11)

其中,I是单位矩阵。根据上述式子,便得到了MDSS方法的迭代矩阵,如下式所示:

K ( α , β ) = ( α H + β D ) 1 ( α D β H ) ( α D + β H ) 1 ( β D α H ) (12)

对此,我们得到迭代矩阵 K ( α , β ) 的谱半径为

ρ ( K ( α , β ) ) = max | ( α λ β ) ( α β λ ) ( α λ + β ) ( α + β λ ) | < 1 , α , β > 0 (13)

其中 S p ( D H 1 ) 表示矩阵 D H 1 的谱,则MDSS迭代(11)在任意初始条件下都收敛于矩阵方程(1)的唯一精确解 X m × n

证明:通过使用Kronecker积,我们可以按照以下的计算结构编写算法(3.1):

{ ( I ( α W + β T ) + ( α U + β V ) T I ) X ( k + 1 2 ) = ( I i ( β W α T ) + i ( β U α V ) I ) X ( k ) + C 1 ( I ( α T + β W ) + ( α V + β U ) I ) X ( k + 1 ) = ( I i ( α W β T ) + i ( α U β V ) I ) X ( k + 1 2 ) + C 2 (14)

接着,根据算法演练,上述两个等式,可以分别重写为

( α D + β H ) x ( k + 1 2 ) = i ( β D α H ) x ( k + 1 ) + C 1 (15)

( α H + β D ) x ( k + 1 ) = i ( α D β H ) x ( k + 1 2 ) + C 2 (16)

其中,我们设定 C 1 = v e c ( C 1 ) C 2 = v e c ( C 2 ) x ( k ) = v e c ( x ( k ) ) 。根据上述等式转换,我们得到

Q 1 = ( α D + β H ) 1 ( α H β D ) (17)

以及

Q 2 = ( α H + β D ) 1 ( α D β H ) (18)

因此,算法MDSS过程的迭代矩阵为

K ( α , β ) = Q 2 Q 1 = ( α H + β D ) 1 ( α D β H ) ( α D + β H ) 1 ( α H β D )

接着,我们假设 λ p , q H λ p T λ q V 分别表示H,T和V的特征值( p = 1 , , m q = 1 , , n )。可以得到

λ p , q H = λ p T + λ q V > 0 , p = 1 , , m ; q = 1 , , n

其中H是非奇异的,因此,对于算法MDSS的迭代矩阵 K ( α , β ) 的谱半径 ρ ( K ( α , β ) ) ,我们可以写出来,得到以下结果

ρ ( K ( α , β ) ) = ρ ( ( α H + β D ) 1 ( α D β H ) ( α D + β H ) 1 ( α H β D ) ) = ρ ( H 1 ( α I + β D H 1 ) 1 ( α D H 1 β I ) H H 1 ( α D H 1 + β I ) 1 ( α I β D H 1 ) H ) = ρ ( ( α I + β D H 1 ) 1 ( α D H 1 β I ) ( α D H 1 + β I ) 1 ( α I β D H 1 ) )

= max λ S p ( D H 1 ) | α β λ α + β λ × α λ β α λ + β | = max λ S p ( D H 1 ) | α 2 λ + β 2 λ α β α β λ 2 α 2 λ + β 2 λ + α β + α β λ 2 | = max λ S p ( D H 1 ) | α λ [ ( α + β 2 α ) β ( λ + 1 λ ) ] α λ [ ( α + β 2 α ) + β ( λ + 1 λ ) ] | = max λ S p ( D H 1 ) | ( α + β α ) β ( λ + 1 λ ) ( α + β α ) + β ( λ + 1 λ ) | (19)

根据上述一系列等式,其中等式中任何一个 α , β , λ > 0 ,因此我们可以得到下列不等式

( α + β α ) β ( λ + 1 λ ) < ( α + β α ) + β ( λ + 1 λ ) (20)

因此,我们得到,迭代矩阵 K ( α , β ) 的谱半径 ρ ( K ( α , β ) ) 有以下关系式:

ρ ( K ( α , β ) ) = max λ S p ( D H 1 ) | ( α λ β ) ( α β λ ) ( α λ + β ) ( α + β λ ) | < 1 , α , β > 0 (21)

根据上述不等式(22),表明了MDSS迭代(11)在任意初始条件下都收敛于矩阵方程(1)的唯一精确解 X * m × n 。定理得证。

根据上述定理,我们得到该算法可以在任意初始条件下都收敛于唯一的精确解。对此,我们还需得到参数,使得算法的谱半径 ρ ( K ( α , β ) ) 取得最小值。因此,我们需要得到定理,并证明定理的可行性。

定理2.2 此外,使函数 ρ ( K ( α , β ) ) 最小化的参数 α , β

α β = u v + u v 4 2 (22)

其中

u = min λ min ζ λ max { ζ + ζ 1 } v = max λ min ζ λ max { ζ + ζ 1 } (23)

并设为

λ min = min λ j S p ( D H 1 ) { λ j } , λ max = max λ j S p ( D H 1 ) { λ j } (24)

接着,我们得到收敛因子 ρ ( K ( α , β ) )

ρ ( K ( α , β ) ) = κ 1 κ + 1 (25)

其中,

κ = u v (26)

证明:从上述所得到的等式与结果上可以看出

ρ ( K ( α , β ) ) = max { α + β 2 α 1 β u α + β 2 α 1 + β u , α + β 2 α 1 β v α + β 2 α 1 + β v } (27)

接着,我们可以根据下列等式得到其最小值

α + β 2 α 1 β u α + β 2 α 1 + β u = α + β 2 α 1 β v α + β 2 α 1 + β v (28)

当上述等式成立时,我们便可以得到谱半径 ρ ( K ( α , β ) ) 的最小值。根据上述关系,我们得到

( α + β 2 α 1 ) 2 = β 2 u v (29)

接着得到下列等式

α 2 α β u v + β 2 = 0 (30)

因此,通过求解上述等式,对于谱半径 ρ ( K ( α , β ) ) ,我们得到

ρ ( K ( α , β ) ) = α + β 2 α 1 β v α + β 2 α 1 + β v (31)

根据上述结果以及等式,我们最终可以得到最小化的谱半径的值

ρ ( K ( α , β ) ) = u v u u v + u = v u 1 v u + 1 = κ 1 κ + 1 (32)

定理得证。

3. 数值例子

为了测试文中提出的新方法的有效性,我们给出了一个求解西尔维斯特矩阵等式的测试问题(1.1)的方法。所有的计算和绘图都是在MATLAB软件上完成的。我们使用X(0) = 0 (零矩阵)来进行初始迭代,如果当前的迭代满足以下条件,则终止迭代

R ( X ( K ) ) F R ( X ( 0 ) ) F 10 10

其中 R ( X ( K ) ) = C A X ( K ) X ( K ) B 。在测试过程中,我们会报告迭代次数(IT)和计算时间(以秒为单位)。对于西尔维斯特矩阵方程 A X + X B = C ,与 A = W + i T B = U + i V ,其中,我们假设

W = U = ( K + 3 3 θ I n ) , T = V = ( K + 3 + 3 θ I n )

其中,矩阵K的形式为 K = I n V m + V m I m ,其中 V m = h 2 T r i ( 1 , 2 , 1 ) m × m 。因此,K是一个大小为n × n的对角矩阵,具有 n = m 2 I n I m 分别是维数为n和m的单位矩阵。在这里,我们设置了

θ = h = 1 m + 1

此外,对于矩阵C,我们设置C = K。同时,我们通过把矩阵方程(1)乘以 h 2 来将它归一化。最后,我们用 α β 表示了MDSS方法中的最优参数。这些参数通过实验被指定,使每种方法的迭代次数最少。

在数值实验中,我们设置 m = 2 。通过输入上述参数,在MATLAB中得到如下图形:

图1可以看出,本文提出的MDSS算法在任意的初始条件下,可以收敛于唯一的数值解,说明算法在实际应用中的可行性也可以得到验证。

Figure 1. The number of iterations of MDSS algorithm under logarithmic conditions

图1. MDSS算法在对数条件下的迭代次数

4. 结论

在本文中,我们在双步尺度分裂(DSS)方法的基础上通过改进,得到MDSS方法,来求解复数域上的连续时间的西尔维斯特矩阵方程的近似解,证明了该迭代序列在任意初始条件的情况下都收敛于西尔维斯特矩阵方程的唯一解,并确定了其最优参数和相应的最优收敛因子。最后,通过数值实验可以证明该算法的可行性。

参考文献

[1] Zhou, K., Doyle, J.C. and Glover, K. (1996) Robust and Optimal Control. Prentice-Hall, Upper Saddle River, New Jersey.
[2] Anderson, B.D.O., Agathoklis, P., Jury, E.I. and Mansour, M. (1986) Stability and the Matrix Lyapunov Equation for Discrete 2-Dimensional Systems. IEEE Transactions on Circuits and Systems, 33, 261-267.
https://doi.org/10.1109/TCS.1986.1085912
[3] Calvetti, D. and Reichel, L. (1996) Application of ADI Iterative Methods to the Restoration of Noisy Images. SIAM Journal on Matrix Analysis and Applications, 17, 165-186.
https://doi.org/10.1137/S0895479894273687
[4] Bru, R., Mas, J. and Urbano, A.M. (1994) An Algorithm for the Single-Input Pole Assignment Problem. SIAM Journal on Matrix Analysis and Applications, 15, 393-407.
https://doi.org/10.1137/S0895479889161125
[5] Baur, U. and Benner, P. (2008) Cross-Gramian Based Model Reduction for Data-Sparse Systems. Electronic Transactions on Numerical Analysis, 31, 256-270.
https://doi.org/10.1137/070711578
[6] Zhang, Y.N., Jiang, D.C. and Wang, J. (2002) A Recurrent Neural Network for Solving Sylvester Equation with Time-Varying Coefficients. IEEE Transactions on Neural Networks, 13, 1053-1063.
https://doi.org/10.1109/TNN.2002.1031938
[7] Liao, A.P., Bai, Z.Z. and Lei, Y. (2005) Best Approximate Solution of Matrix Equation AXB CYD = E. SIAM Journal on Matrix Analysis and Applications, 27, 675-688.
https://doi.org/10.1137/040615791
[8] Ilic, M. (1989) New Approaches to Voltage Monitoring and Control. IEEE Control Systems Magazine, 9, 5-11.
https://doi.org/10.1109/37.16743
[9] Bai, Z.Z., Huang, Y.M. and Ng, M.K. (2007) On Preconditioned Iterative Methods for Burgers Equations. SIAM Journal on Scientific Computing, 29, 415-439.
https://doi.org/10.1137/060649124
[10] Bartels, R.H. and Stewart, G.W. (1972) Solution of the Matrix Equation AX XB = C. Communications of the ACM, 15, 820-826.
https://doi.org/10.1145/361573.361582
[11] Golub, G.H., Nash, S. and Loan, C.V. (1979) A Hessenberg-Schur Method for the Problem AX XB = C. IEEE Transactions on Automatic Control, 24, 909-913.
https://doi.org/10.1109/TAC.1979.1102170
[12] Benner, P., Li, R.C. and Truhar, N. (2009) On the ADI Method for Sylvester Equations. Journal Computational and Applied Mathematics, 233, 1035-1045.
https://doi.org/10.1016/j.cam.2009.08.108
[13] Bai, Z.-Z. (2011) On Hermitian and Skew-Hermitian Spliting Iteration Methods for Continuous Sylvester Equations. Journal of Computation and Mathematics, 29, 185-198.
https://doi.org/10.4208/jcm.1009-m3152
[14] Smith, R.A. (1968) Matrix Equation XA BX = C. SIAM Journal on Applied Mathematics, 16, 198-201.
https://doi.org/10.1137/0116017
[15] Bai, Z.-Z., Benzi, M. and Chen, F. (2010) Modified HSS Iteration Methods for a Class of Complex Symmetric Linear Systems. Computing, 87, 93-111.
https://doi.org/10.1007/s00607-010-0077-0
[16] Hezari, D., Salkuyeh, D.K. and Edalatpour, V. (2016) A New Iterative Method for Solving a Class of Complex Symmetric System of Linear Equations. Numerical Algorithms, 73, 927-955.
https://doi.org/10.1007/s11075-016-0123-x