Advances in Applied Mathematics
Vol. 11  No. 03 ( 2022 ), Article ID: 49893 , 9 pages
10.12677/AAM.2022.113162

Hermitian正定Toeplitz线性方程组的外推CSCS方法

傅毛里,刘仲云

长沙理工大学,数学与统计学院,湖南 长沙

收稿日期:2022年2月28日;录用日期:2022年3月23日;发布日期:2022年3月30日

摘要

本文基于Toeplitz矩阵有循环与反循环分裂(CSCS)的事实,提出求解Hermitian正定Toeplitz线性方程组的外推CSCS方法,并给出其最优双参数 α , β ,以及最优外推参数 ω 的选取策略。并通过数值实验验证我们方法的有效性。

关键词

外推法,Toeplitz矩阵,循环与反循环分裂,双参数,Hermitian正定矩阵

Extrapolated Circulant and Skew Circulant Splitting Method for Hermitian Positive Definite Toeplitz Systems

Maoli Fu, Zhongyun Liu

School of Mathematics and Statistics, Changsha University of Science and Technology, Changsha Hunan

Received: Feb. 28th, 2022; accepted: Mar. 23rd, 2022; published: Mar. 30th, 2022

ABSTRACT

Based on the fact that a Toeplitz matrix admits a circulant and skew circulant splitting (CSCS), we propose the extrapolated CSCS method to solve Hermitian definite Toeplitz systems and discuss the strategy to select the optimal two parameters α , β , and the extrapolated parameter ω . The effectiveness of our method is verified by numerical experiments.

Keywords:Extrapolated Method, Toeplitz Matrices, Circulant and Skew Circulant Splitting, Two Parameters, Hermitian Definite Matrices

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

本文考虑用迭代法求解线性方程组

T x = b , (1.1)

其中 x , b n T n 是Hermitian正定Toeplitz矩阵。

在信号处理和控制论 [1] [2] [3] 等许多应用中经常遇到这类Toeplitz线性方程组。关于这类线性方程组的解法,主要分为直接法和迭代法两大类 [4] [5] [6] [7]。由于直接法通常存在数值不稳定的缺点 [8],这使得迭代法更加受到人们的关注。

设系数矩阵 T = M N ,若M非奇异,则可推导出如下求解线性方程组(1.1)的一个迭代序列

x ( k + 1 ) = G x ( k ) + M 1 b, k = 0 , 1 , , (1.2)

其中 G = M 1 N 为迭代矩阵。

对于非奇异线性方程组(1.1),迭代(1.2)收敛当且仅当

ρ ( G ) max { | μ | | μ σ ( G ) } < 1 ,

其中 σ ( G ) 是矩阵G的谱集, max { } 是取最大元素。

对迭代(1.2),常用的加速技巧是外推法 [9],其迭代格式为

x ( k + 1 ) = G ω x ( k ) + ω M 1 b , (1.3)

其中 G ω = ( 1 + ω ) I + ω M 1 N ,I是n阶单位矩阵。在(1.3)中可以看到,若 ω = 1 ,则(1.3)就退化成了(1.2)的迭代方法。

对于任意Toeplitz矩阵 T n 皆有以下分裂

T n = C n + S n , (1.4)

其中

C n = 1 2 [ t 0 t 1 + t n 1 t 2 n + t 2 t 1 n + t 1 t 1 + t 1 n t 0 t 3 n + t 3 t 2 n + t 2 t n 2 + t 2 t n 3 + t 3 t 0 t 1 + t n 1 t n 1 + t 1 t n 2 + t 2 t 1 + t 1 n t 0 ] ,

S n = 1 2 [ t 0 t 1 t n 1 t 2 n t 2 t 1 n t 1 t 1 t 1 n t 0 t 3 n t 3 t 2 n t 2 t n 2 t 2 t n 3 t 3 t 0 t 1 t n 1 t n 1 t 1 t n 2 t 2 t 1 t 1 n t 0 ] ,

显然, C n 是一个Hermitian循环矩阵, S n 是一个Hermitian反循环矩阵,分裂(1.4)为 T n 的一个CSCS。

基于CSCS(1.4),Ng提出了如下的CSCS [10]

{ ( α I + C ) x ( k + 1 / 2 ) = ( α I S ) x ( k ) + b , ( α I + S ) x ( k + 1 ) = ( α I C ) x ( k + 1 / 2 ) + b, (1.5)

其中 α 是一个给定的正常数。作者还证明了对于任意正定Toeplitz线性方程组CSCS迭代(1.5)都收敛。此外,作者还给出了CSCS迭代收缩因子的上界,它仅依赖于C和S的最大特征值与最小特征值。2012年Akhondi等人在CSCS迭代的基础上提出了双参数的加速版本,简称ACSCS [11] 迭代,其格式为

{ ( α I + C ) x ( k + 1 / 2 ) = ( α I S ) x ( k ) + b, ( β I + S ) x ( k + 1 ) = ( β I C ) x ( k + 1 / 2 ) + b . (1.6)

(1.6)可写成如下的单步迭代

x ( k + 1 ) = R x ( k ) + b ˜ , (1.7)

其中 R = ( β I + S ) 1 ( β I C ) ( α I + C ) 1 ( α I S ) b ˜ = ( β I + S ) 1 b + ( β I + S ) 1 ( β I C ) ( α I + C ) 1 b

ACSCS迭代的计算主要是循环矩阵和反循环矩阵与向量的乘积,这可以使用快速Fourier变换 [12] (FFT)实现,其原理是循环矩阵利用Fourier矩阵进行对角化,反循环矩阵利用Fourier矩阵与一个对角矩阵的乘积实现对角化,

F ˜ = F × Ω ,

其中F是Fourier矩阵, Ω 是对角矩阵,其主对角的元素为 ( 1 , e π i / n , , e ( n 2 ) π i / n , e ( n 1 ) π i / n ) F ˜ 为酉矩阵,则有

C = F * Λ 1 F , S = F ˜ * Λ 2 F ˜ , (1.8)

其中 Λ 1 Λ 2 都是对角矩阵,对角元分别为循环矩阵C和反循环矩阵S的特征值,并且循环矩阵与反循环矩阵的特征值只需一次FFT就可以精确算出来。

2. 预备知识

接下来给出文中用到的基本定义及引理。

定义2.1 [13] 若矩阵 T n n × n 满足

S n = 1 2 [ t 0 t 1 t 2 n t 1 n t 1 t 0 t 3 n t 2 n t n 2 t n 3 t 0 t 1 t n 1 t n 2 t 1 t 0 ] ,

则称 T n 为Toeplitz矩阵,若还满足 t k = t ¯ k ( k = 1 , 2 , , n 1 ) ,则称 T n 为Hermitian Toeplitz矩阵。

定义2.2 [14] 设 C 2 π 表示定义在区间 [ π , π ] 内所有以 2 π 为周期的连续实值函数,对所有 f ( x ) C 2 π ,若

t k = 1 2 π π π f ( x ) e i k x d x ( k = 0 , ± 1 , ± 2 , )

f ( x ) 的Fourier系数,则称 f ( x ) 为Toeplitz矩阵 T n 的生成函数。

定义2.3 [14] 若 C n n × n 是Toeplitz矩阵,且满足 c k = c n k ( k = 1 , 2 , , n 1 ) ,则称 C n 为循环矩阵。

定义2.4 [14] 若 S n n × n 是Toeplitz矩阵,且满足 s k = s n k ( k = 1 , 2 , , n 1 ) ,则称 S n 为反循环矩阵。

定义2.5 [14] 若 F n n × n 满足

( F n ) j , k = 1 n e 2 π i j k n ,

则称 F n 为Fourier矩阵。其中 i = 1 为虚数单位,e为自然对数的底, π 为圆周率。

引理2.1 [11] 设 T n 是Hermitian正定Toeplitz矩阵, C n , S n 为(1.4)所定义的矩阵,只要n足够大,则 C n S n 的特征值均是正实的。

3. 外推CSCS迭代方法

为了进一步加速ACSCS迭代(1.7),我们提出了以下求解线性方程组(1.1)的迭代

{ x ˜ ( k + 1 ) = R x ( k ) + b ˜ , x ( k + 1 ) = ω x ˜ ( k + 1 ) + ( 1 ω ) x ( k ) , (3.1)

我们称(3.1)为外推循环与反循环分裂迭代方法,简称EACSCS,其中 α , β , ω 是给定的正常数。由(1.8)可知,循环矩阵和反循环矩阵与向量乘积可以利用快速Fourier变换实现,故EACSCS每次迭代的计算量为 O ( n log n )

迭代(3.1)可以理解为直接分裂迭代,主要是便于收敛性分析,但从实际计算效果来看,该迭代不如下面的基于残差校正的迭代 [15]

{ x ( k + 1 / 2 ) = x ( k ) + ( α I + S ) 1 ( b T x ( k ) ) , x ˜ ( k + 1 ) = x ( k + 1 / 2 ) + ( β I + S ) 1 ( b T x ( k + 1 / 2 ) ) , x ( k + 1 ) = ω x ˜ ( k + 1 ) + ( 1 ω ) x ( k ) (3.2)

接下来我们分析迭代(3.1)的收敛性。在不引起歧义的情况下,我们将矩阵维数的下标省略。设T是一个Hermitian正定Toeplitz矩阵,由引理2.1可知C和S的特征值均是正实的,分别记为 0 < λ 1 < λ 2 < < λ n 0 < μ 1 < μ 2 < < μ n

下面我们讨论 α , β 的选取,为此先引入下面引理。

引理3.1 [16] 设T是一个Hermitian正定Toeplitz矩阵,R为(1.7)中所定义的迭代矩阵,则有

ρ ( R ) max 1 j n | α μ j β + μ j | max 1 j n | β λ j α + λ j | = ϕ ( α , β ) (3.3)

引理3.2设T是一个Hermitian正定Toeplitz矩阵,则有

max 1 j n | α μ j β + μ j | = { μ n α β + μ n , α α ^ , α μ 1 β + μ 1 , α > α ^ , μ 1 α ^ μ n , (3.4)

max 1 j n | β λ j α + λ j | = { λ n β α + λ n , β β ^ , β λ 1 α + λ 1 , β > β ^ . λ 1 β ^ λ n , (3.5)

证明 现在定义

f 1 ( μ ) = α μ β + μ , μ [ μ 1 , α ] , f 2 ( μ ) = μ α β + μ , μ [ α , μ n ] ,

可得

f 1 ( μ ) = β α ( β + μ ) 2 < 0 , f 2 ( μ ) = β + α ( β + μ ) 2 > 0 ,

f 1 ( μ ) 单调递减, f 2 ( μ ) 单调递增,所以 f 1 ( μ ) 的最大值为 f 1 ( μ 1 ) f 2 ( μ ) 的最大值 f 2 ( μ n ) ,即(3.4)得证,同理可证得(3.5)。

定理3.2设T是一个Hermitian正定Toeplitz矩阵,若 α , β 取值为

α ^ = μ ^ p λ ^ p + Δ μ ^ s + λ ^ s , (3.6)

β ^ = λ ^ p μ ^ p + Δ μ ^ s + λ ^ s , (3.7)

ϕ ( α , β ) 达到最小,其值为

ϕ ( α ^ , β ^ ) = θ 1 θ + 1 , (3.8)

其中 μ ^ p = μ 1 μ n λ ^ p = λ 1 λ n μ ^ s = μ 1 + μ n λ ^ s = λ 1 + λ n Δ = ( μ ^ p λ ^ p ) 2 + ( μ ^ s + λ ^ s ) ( μ ^ s λ ^ p λ ^ s μ ^ p ) θ = ( λ n + μ 1 ) ( λ 1 + μ n ) / [ ( λ n + μ n ) ( λ 1 + μ 1 ) ]

证明 有引理3.2可知,存在 α ^ [ μ 1 , μ n ] β ^ [ λ 1 , λ n ] 有(3.4),(3.4),为了求得 α ^ , β ^ 使得(3.3)中 ϕ ( α , β ) 达到最小,那么需要 α ^ , β ^ 让(3.4),(3.5),同时达到最小,即

{ μ n α β + μ n = α μ 1 β + μ 1 , λ n β α + λ n = β λ 1 α + λ 1 , (3.9)

利用 μ ^ p , λ ^ p , μ ^ s , λ ^ s ,则(3.9)可化简为

{ 2 ( α β μ ^ p ) + ( α β ) μ ^ s = 0 , 2 ( α β λ ^ p ) ( α β ) λ ^ s = 0. (3.10)

经过简单计算就可得出(3.6)和(3.7)。将 α ^ , β ^ 代入 ϕ ( α , β )

ϕ ( α ^ , β ^ ) = μ n α ^ β ^ + μ n × λ n β ^ α ^ + λ n = μ n ( μ ^ s + λ ^ s ) μ ^ p + λ ^ p Δ μ n ( μ ^ s + λ ^ s ) μ ^ p + λ ^ p + Δ × λ n ( μ ^ s + λ ^ s ) λ ^ p + μ ^ p Δ λ n ( μ ^ s + λ ^ s ) λ ^ p + μ ^ p + Δ = ( λ n + μ n ) ( λ 1 + μ n ) Δ ( λ n + μ n ) ( λ 1 + μ n ) + Δ × ( λ n + μ n ) ( λ n + μ 1 ) Δ ( λ n + μ n ) ( λ n + μ 1 ) + Δ = ( λ n + μ 1 ) ( λ 1 + μ n ) Δ ( λ n + μ 1 ) ( λ 1 + μ n ) + Δ = ( λ n + μ 1 ) ( λ 1 + μ n ) / [ ( λ n + μ n ) ( λ 1 + μ 1 ) ] 1 ( λ n + μ 1 ) ( λ 1 + μ n ) / [ ( λ n + μ n ) ( λ 1 + μ 1 ) ] + 1 .

虽然通过定理3.2可以找到 α ^ , β ^ 使得收缩上界达到最小,但不是谱半径 ρ ( R ) 达到最小的参数,数值实验表明在 α ^ , β ^ 的附近可以搜索找到使 ρ ( R ) 达到较小的数值最优 α * , β *

一旦 α * , β * 确定,接下来可以考虑外推参数的选取,关于外推参数的选取建议读者参阅文献 [9],文献中给出了外推法的上界最优参数,其选取方式如下

ω ^ = { 1 η n ( 1 η n ) 2 + τ n 2 , δ 1 δ 2 , 2 2 ( η 1 + η n ) , δ 1 > δ 2 , (3.11)

其中 η 1 , η 2 为R的特征值实部最小值和最大值, τ n 为R的特征值虚部最大值, δ 1 = ( η n η 1 ) ( 1 η n ) δ 1 = 2 τ n 2

由(3.11)可知,选取最优外推参数 ω 时,只与矩阵R的特征值的虚部最大值以及实部最小与最大值有关。我们在外推参数 ω 的理论最优值附近进行搜索以此得到数值最优的 ω *

4. 数值实验

我们使用EACSCS迭代求解三个例子,并与CSCS迭代和ACSCS迭代作比较。测试矩阵的阶数 n = 64 , 128 , 256 , 512 , 1024 ,迭代的初始向量 x ( 0 ) 是零向量,方程(1.1)右端向量b是元素均为1的向量,迭代的停机准则为:

ε = r ( k ) 2 r ( 0 ) 2 10 7 ,

其中 r ( k ) = b T x ( k ) 是第k步迭代的残差向量。所有数值实验均在Matlab R2018b上进行测试。

以下测试例子均来自文献 [14]。

例子4.1进行Hermitian正定Toeplitz矩阵T的生成函数为 f ( x ) = x 4 + 1 , x [ π , π ] ,其元素为

t k = { 4 π 2 ( 1 ) k k 2 24 ( 1 ) k k 2 , k 0 , π 4 5 + 1 , k = 0.

例子4.2进行Hermitian正定Toeplitz矩阵T的生成函数为 f ( x ) = k = 0 [ e i k ( log k + x ) / k + e i k ( log k + x ) / k ] + 4.2 x [ π , π ] ,其元素为

t k = { e 1 log k k , k > 0 , 4.2 , k = 0 , t ¯ k , k < 0.

例子4.3进行Hermitian正定Toeplitz矩阵T的生成函数为 f ( x ) = 2 k = 0 [ ( sin ( k x ) + cos ( k x ) ) / ( 1 + k ) 1.1 ] x [ π , π ] ,其元素为

t k = { 1 + 1 ( 1 + k ) 1.1 , k > 0 , 2 , k = 0 , t ¯ k , k < 0.

下表中只列出了数值最优参数 α * , β * , ω * ,其中计算 ω * 的过程中使用到了幂迭代,其迭代次数不超过3次。算法的迭代次数记为IT,算法的运行时间记为CPU (单位为秒)。

Table 1. Numerical results of example 4.1

表1. 例子4.1的数值结果

Table 2. Numerical results of example 4.2

表2. 例子4.2的数值结果

Table 3. Numerical results of example 4.3

表3. 例子4.3的数值结果

表1~3可知,EACSCS方法的迭代次数和运算时间均比其它两种方法少,达到了较好的加速目的。

如上,我们绘制了EACSCS,ACSCS,CSCS的迭代步数与相对残差的变化图,并且标出了停机准则 ε 线,分别用 Δ Δ Δ 曲线表示。

图1可知,EACSCS迭代的相对残差下降速度明显快于其它两种方法,有较好的加速效果。

Figure 1. When n = 1024, the number of iteration steps and relative residual changes of the three iterative methods in Example 4.1~4.3

图1. 当n = 1024时,例4.1~4.3中三种迭代方法的迭代步数与相对残差变化图

5. 结束语

基于Toeplitz矩阵有CSCS的事实,本文提出了外推带双参数的CSCS方法,通过理论分析,我们给出了方法中最优外推参数 α , β 以及外推参数 ω 的选取策略。利用数值实验,证明了我们的方法要明显优于ACSCS和CSCS。

文章引用

傅毛里,刘仲云. Hermitian正定Toeplitz线性方程组的外推CSCS方法
Extrapolated Circulant and Skew Circulant Splitting Method for Hermitian Positive Definite Toeplitz Systems[J]. 应用数学进展, 2022, 11(03): 1484-1492. https://doi.org/10.12677/AAM.2022.113162

参考文献

  1. 1. Chui, C. and Chan, A. (1982) Application of Approximation Theory Methods to Recursive Digital Filter Design. IEEE Transactions on Acoustics Speech & Signal Processing, 30, 18-24. https://doi.org/10.1109/TASSP.1982.1163842

  2. 2. Gilliam, D., Martin, C. and Lund, J. (1987) Analytic and Numerical Aspects of the Observation of the Heat Equation. 26th IEEE Conference on Decision & Control, Los Angeles, 9-11 December 1987, 975-976. https://doi.org/10.1109/CDC.1987.272541

  3. 3. King, R., Ahmadi, M., Gorgui, R., et al. (1989) Digital Filtering in One and Two Dimensions. Plenum Press, New York. https://doi.org/10.1007/978-1-4899-0918-3

  4. 4. Axelsson, O. (2003) Iterative Solution Methods. Cambridge University Press, New York.

  5. 5. Cao, Z.H. (1995) Convergence of Two-Stage Iterative Methods for the Solution of Linear Systems. Mathematica Numerica Sinica, 17, 98-109.

  6. 6. Ruggiero, V. and Galligani, E. (1990) An Iterative Method for Large Sparse Linear Systems on a Vector Computer. Computers& Mathematics with Applications, 20, 25-28. https://doi.org/10.1016/0898-1221(90)90065-R

  7. 7. Saad, Y. (2003) Iterative Methods for Sparse Linear Systems. 2nd Edition, SIAM, Philadelphia. https://doi.org/10.1137/1.9780898718003

  8. 8. Golub, G.H. and Loan, A. (2013) Matrix Computations. 4th Edition, Johns Hopkins University Press, Baltimore.

  9. 9. Yeyios, A. (1984) On the Optimization of an Extrapolation Method. Linear Algebra and Its Applications, 57, 191-203. https://doi.org/10.1016/0024-3795(84)90187-3

  10. 10. Ng, M.K. (2003) Circulant and Skew-Circulant Splitting Methods for Toeplitz Systems. Journal of Computational and Applied Mathematics, 159, 101-108. https://doi.org/10.1016/S0377-0427(03)00562-4

  11. 11. Akhondi, N. and Toutounian, F. (2012) Accelerated Circulant and Skew Circulant Splitting Methods for Hermitian Positive Definite Toeplitz Systems. Advances in Numerical Analysis, 2012, Article ID: 973407. https://doi.org/10.1155/2012/973407

  12. 12. Chan, R. and Ng, M.K. (1996) Conjugate Gradient Methods for Toeplitz Systems. SIAM Review, 38, 427-482. https://doi.org/10.1137/S0036144594276474

  13. 13. Horn, R.A. and Johnson, C.R. (2013) Matrix Analysis. 2nd Edition, Cambridge University Press, New York.

  14. 14. Chan, R. and Jin, X. (2007) An Introduction to Iterative Toeplitz Solvers. SIAM, Philadelphia. https://doi.org/10.1137/1.9780898718850

  15. 15. Bai, Z.Z. and Rozloznik, M. (2015) On the Numerical Behavior of Matrix Splitting Iteration Methods for Solving Linear Systems. SIAM Journal on Numerical Analysis, 53, 1716-1737. https://doi.org/10.1137/140987936

  16. 16. Liu, Z.Y., Zhou, Y. and Zhang, Y. (2020) On Inexact Alternating Direction Implicit Iteration for Continuous Sylvester Equations. Numerical Linear Algebra with Applications, 27, e2320. https://doi.org/10.1002/nla.2320

期刊菜单