大型连续Sylvester 方程外推的CSCS 迭代

刘仲云 张芳

(长沙理工大学数学与统计学院,长沙,410114)

1 引言

大型连续Sylvester 方程[1]是最常用的线性矩阵方程之一.矩阵方程形如

这里A ∈Rn×n,B ∈Rm×m,E ∈Rn×m是给定的矩阵.众所周知,当系数矩阵A和−B没有公共特征值时,方程(1.1)有唯一解[2,3].

Sylvester 方程是设计Luenberger 观测的经典方法, 它广泛应用于信号处理、控制和系统理论[4–7],也常出现在计算不变量子空间的Riccati方程的线性和广义特征值问题中[8–10],也可用于设计常微分方程数值解的隐式龙格 库塔积分公式和分块多步公式中[11].矩形域上分离椭圆边值问题的有限差分离散产生的方程也可以写成Sylvester 方程[12,13],优化理论、结构动力学、固体力学、地质分子光学等领域都涉及到Sylvester 方程的求解[14].

求解连续Sylvester 方程(1.1)主要有两种方法.第一种方法是将矩阵方程(1.1)转换成一个线性方程Φx=c,这里Φ=Im ⊗A+BT ⊗In是矩阵A和BT的Kronecker 积的和,向量x和c是矩阵X和C的列向量,直接法和迭代法均可用于求解该线性方程. 第二种方法是对方程(1.1)进行迭代.

当系数矩阵阶数比较大时,线性方程Φx=c中的系数矩阵Φ 的阶数会相当大,一般会出现数据存储和计算时间方面的困难.这说明第一种方法主要用于中小维问题.例如: 文献[15]中提出的Bartel stewart 方法是基于使用特征值QR 算法将矩阵A和B化为实Schur 形式,然后使用直接法来求解几个线性方程.文献[16]提出的Hessenberg Schur 法是将矩阵A化为Hessenberg 形式,将矩阵B分解为拟三角Schur 形式,其速度比Bartels Stewart 方法快.这些方法都属于直接法,并在MATLAB 中进行了应用.

对于大型稀疏Sylvester 方程,一般用迭代法求解. 比较常见的有古典迭代法,Smith 迭代法[17],ADI 迭代法[18]等. 最近比较流行的迭代法, 是利用系数矩阵A和B的分裂来建立类似ADI 迭代格式的迭代法,例如,白中治在2011 年提出的HSS 迭代法[19].自该方法提出后,几个更为有效的类似的分裂迭代法相继出现,例如: 2013 年王湘提出的PSS 迭代法[20], 2014 年郑青青提出的NSS 迭代法[21]等.这些方法都无条件收敛到方程的精确解.

文献[1]考虑了(1.1)中系数矩阵A和B都是正定Toeplitz 矩阵的情况,提出了CSCS 迭代法.通过这种CSCS 迭代方法,将求解一般连续Sylvester 方程的问题转化为包含两个循环矩阵和反循环矩阵的连续Sylvester 方程的求解问题.通过这种转化可以利用傅里叶算法快速求解.在对流扩散方程的离散化中就出现了这种结构的矩阵[13,22].

上述这些情况激发了我们对连续Sylvester 方程进一步的研究兴趣. 本文提出一个外推的CSCS(ECSCS)迭代.文章的结构如下: 第二节介绍基本定义,第三节介绍CSCS 迭代,第四节给出ECSCS 迭代法,并分析其收敛性,第五节通过数值实验验证该方法的有效性.

2 基本定义

对于任意的矩阵K ∈Rn×n,K∗,ρ(K)和λ(K)分别表示它的共轭转置矩阵、谱半径和特征值.Ki,j是矩阵K的第i行j列元素.若x ∈R,则x的实部记为Re(x),虚部记为Im(x).如果矩阵(K+K∗)是正定的,则矩阵K也是正定的.该正定条件也等价于: 存在一个非零向量z ∈R,使得Re(z∗Kz)>0,即对于矩阵K的任意特征值λ,有Re(λ)>0[23].

任意一个Toeplitz 矩阵T都有一个循环反循环分裂,即:T=CT+ST,这里CT是循环矩阵,ST是反循环矩阵,定义如下:

众所周知,循环矩阵C和反循环矩阵S可由傅里叶矩阵F和酉矩阵ˆF=FD进行如下对角化:

3 CSCS 迭代

若方程(1.1)中的系数矩阵A ∈Rn×n,B ∈Rm×m是正定Toeplitz 矩阵,则矩阵A,B可以进行如下循环反循环分裂:

给定常数α,β及适当维数的单位矩阵I,则有

由此可知,连续Sylvester 方程(1.1)可以等价地写成不动点矩阵方程

文献[1]中提出的求解(1.1)的CSCS 迭代定义如下.

CSCS 迭代给定一个初始矩阵X(0)∈Rn×m,使用以下迭代格式计算X(k+1)∈Rn×m(k=0,1,2,···):

对于上面的CSCS 迭代,文献[1]证明了如下收敛性定理.

因而,对任意γ>0,有ρ(T(γ))≤σ(γ)<1.

4 ECSCS 迭代

为了进一步加快CSCS 迭代(3.1)的收敛速度,本节给出ECSCS 迭代.

ECSCS 迭代给定一个初始矩阵X(0)∈Rn×m,使用以下迭代格式计算X(k+1)∈Rn×m(k=0,1,2,···):

用矩阵向量形式,ECSCS 迭代可以等价为

进而变成一步迭代格式

关于ECSCS 迭代,有如下的收敛性定理:

定理2在定理1 的假设下,(4.1)中产生的ECSCS 迭代序列X(k)收敛到(1.1)的精确解X∗,其中ω ∈(0,).

5 数值实验

例1 考虑连续Sylvester 方程(1.1),其中矩阵M,N ∈Rn×n都是三对角矩阵:M= tridiag(−1,2,−1),N= tridiag(−0.5,0,0.5). 数值结果见表1和表2.

表1 SSOR,HSS,CSCS 和ECSCS 的IT 与CPU

表2 SSOR,HSS,CSCS 和ECSCS 的IT 与CPU

例2 考虑Sylvester 方程(1.1),其中矩阵

数值结果见表3.

表3 SSOR,HSS,CSCS 和ECSCS 的IT 与CPU

从例1 和例2 中可以看到: CSCS 迭代在迭代次数和计算效率上都优于HSS 迭代和SSOR 迭代,而ECSCS 比CSCS 收敛的效果更好,并且系数矩阵的阶数越大,效果越明显,因此外推的CSCS迭代很好地提高了CSCS 迭代法的收敛速度.