第33卷第2期南 京 航 空 航 天 大 学 学 报V o l . 33N o. 2
2001年4月 A p r . 2001Jou rnal of N an jing U n iversity of A eronau tics &A stronau tics
文章编号:100522615(2001) 0220139207
求解大规模矩阵问题的K rylov 子空间方法
戴 华
(南京航空航天大学理学院 南京, 210016)
摘要 。最近几年, 其研究工作取得了许多重大进展。K rylov 子空间方法若干进展的一个概述, 。、SY MM LQ 算法、M I N R ES 算法、G M s 、QM R 算法以及这些算法的块格式; 求解大型对称特征值问题的L s L s 算法; 求解大型非对称特征值问题的L anczo s 算法、A rno ldi 算法以及。提出了一些有待进一步研究的问题。关键词:矩阵; 线性方程组; 特征值; K rylov 子空间方法中图分类号:O 241. 6 文献标识码:A
引 言
大规模矩阵问题包括线性方程组、矩阵特征值问题等数值解法的理论研究、算法设计和软件研制是当今计算数学和科学工程计算研究的重大课题, 是大规模科学与工程计算的基础和重要组成部分, 其研究具有重要的理论意义和广泛的应用价值。尤其是其中大型非对称矩阵问题的数值计算, 因其理论上的复杂性, 导致在算法的理论研究和软件设计中产生了许多本质的困难。问题的挑战性吸引了国内外许多数值分析专家从事该领域的研究, 并取得了许多重大进展。
求解线性方程组和矩阵特征值问题的K rylov 子空间方法可以追溯到50年代初L anczo s
[1, 2]
理论分析和数值试验, 充分认识到K rylov 子空间
方法是求解大型线性方程组和大型矩阵特征值问题的一类最有效的方法。本文概述K rylov 子空间方法的若干进展, 其中包括作者对这些问题的研究成果, 并提出一些有待进一步研究的问题。
1 求解大型线性方程组的K rylov 子
空间方法
对线性方程组
A x =b
(1)
其中, A ∈R n ×n , b ∈R n , 设x 0∈R n 是式(1) 解的初始估计, r 0=b -A x 0是初始残量。称K k (A , r 0) =
k -span {r 0, A r 0, …, A
1
r 0}为由A 和r 0产生的
,
K rylov 子空间。K rylov 子空间方法求式(1) 的具有
[3][4]
A rno ldi , H estenes 和Stiefel 的工作。但在其后
20年, 人们一直认为这类方法是数值不稳定的, 很
如下形式的近似解
x k ∈x 0+K k (A , r 0) 并且使残量r k =b -A x k 满足Galerk in 条件
r k ⊥L k
(2) (3)
少用于实际计算。直到70年代初Paige 在其著名的博士论文[5]中重新研究了L anczo s 方法以及80年代初Saad [6]重新研究了A rno ldi 方法, 才使人们重新认识了这两个方法。之后, 人们又做了大量的
其中L k 是R n 中的一个子空间, 或使残量r k 的范数极小化, 即
基金项目:国家自然科学基金(编号:19671043) 、江苏省自然科学基金(编号:BK 97059) 、江苏省“333工程”基金和江苏省“青蓝工程”基金资助项目。
收稿日期:2000204211; 修订日期:2000211201
作者简介:戴 华, 男, 教授, 博士导师, 1959年6月生。
140南 京 航 空 航 天 大 学 学 报第33卷
‖r k ‖=
x ∈x 0+K k (A , r 0)
m in ‖b -A x ‖(4)
T k =
Α1Β20
Β2Α2ω
ω…Β2Α2ω0…
ωωωΒk …ωωωΒk 0
T
用式(2, 3) 定义近似解x k 的K rylov 子空间方法称为正交残量(OR ) 方法。用V k 和W k 分别表示子空间K k (A , r 0) 和L k 的基为列作成的矩阵, OR 方法可行的充要条件是W T 。OR 方法产生k AV k 非奇异的近似解为
x k =x 0+V k (W k AV k ) W k r 0
-1
T
ωω
ω00
Βk Αk 0
Α1Β2
T k =
(e )
(5)
ωω0
用式(2, 4) 定义近似解x k 的K rylov 子空间方法称为极小残量(M R ) 方法。如果最小二乘问题
m in k ‖r 0-AV k y ‖
y ∈R
Βk Αk Βk +(e )
(6)
的解为y , 则M R 方法产生的近似解为
x k =x 0+V k y
k =[v 1, …, v k ],则
AV k =V k T k +Βk +1v k +1e k =V k +1T k
T
如果在式(3) A , r , 方[7]。
在式(3) L k 或在式(4) 中选取不同的范数可导出不同的算法。
1. 1 对称正定线性方程组的共轭梯度法
1952年H estenes 和Stiefel [4]提出了求解对称
(11)
SY MM LQ 算法是在OR 方法中取L k =
K k (A , r 0) , W k =V k , 则W k AV k =T k , x k =x 0+
-1
ΒV k T k e 1; M I N R ES 算法是在式(10) 意义下的M R
方法, 由式(11) , 式(6) 化为m in k ‖r 0-AV k y ‖2=
y ∈R
m in k ‖Βe 1-T k y ‖2, 则容易求得y , 并由式(7) 即
y ∈R
(e )
正定线性方程组的共轭梯度法(CG ) 。该方法可以充分利用矩阵A 的稀疏性, 每次迭代的主要计算量是向量之间的运算。共轭梯度法是一个K rylov 子空间方法。实际上, 共轭梯度法是在向量范数
2
‖x ‖A =(A x , x ) 1 下的M R 方法[7], 近似解x k 满足如下极小化性质‖b -A x k ‖A
-1
得近似解x k 。
SY MM LQ 算法和M I N R ES 算法均以三项递
推公式进行计算, 具有计算量小、存贮量少的优点, 并且M I N R ES 算法的残量范数单调下降。这两个
算法至今仍是求解大型对称不定线性方程组的最有效方法。
对具有相同系数矩阵而具有不同右端项的大型线性方程组集A x i =b i (i =1, …, m ) , 对所有方程组AX =B (其中X =[x 1, …, x m ],B =[b 1, …, b m ]) 采用块方法同时进行迭代比对各个方程组A x i =b i
[10]
使用通常的迭代法更有效。O ′L eary 提出了求解
=
x ∈x 0+K k (A , r 0)
m in ‖b -A x ‖A
-1
(8)
并且[8]
‖x -x k ‖A ≤2‖x -x 0‖A
+k
(9)
大型对称线性方程组集的块CG 方法, 并利用块L anczo s 过程
[33]
其中ϑ=K (A ) =Κm ax (A ) m in (A ) 表示A 的谱条件Κ数。
1. 2 对称不定线性方程组的SY MM LQ 算法和
M I N R ES 算法
将SY MM LQ 算法推广到求解AX
=B , 给出了块L anczo s 方法。但是对所有不同的右
端项由块L anczo s 方法产生的序列收敛速度一般是不同的, 并且残量范数未必单调下降。一个有效的块方法应该能判断、自适应地收缩收敛的方程组, 但块L anczo s 方法的收缩问题一直未解决。作者[11]利用收敛性判据识别和收缩收敛的方程组, 提出了求解大型对称线性方程组集的自适应块L anczo s 方法。为了解决块L anczo s 方法中残量范
对对称不定或非对称矩阵A , 人们通常用如下极小化性质代替式(8)
‖b -A x k ‖2=
x ∈x 0+K k (A , r 0)
m in
‖b -A x ‖2(10)
1975年Paige 和Saunders [9]提出了求解对称不定线性方程组的数值稳定的SY MM LQ 算法和M I N R ES 算法。这两个算法利用对称L anczo s 过
数的振荡, 文[11]提出了M I N R ES 方法的一个块
版本——块M I N R ES 方法, 描述了块M I N R ES 方法的收缩技术, 建立了块L anczo s 方法与块M I N R ES 方法之间的关系。
程[1]产生K k (A , r 0) 的一组标准正交基v 1, …, v k 和
矩阵
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法141
1. 3 G M R ES 及相关算法
对非对称矩阵A , A rno ldi [3]提出了产生K rylov 子空间K k (A , r 0) 的一组标准正交基的方法, 即A rno ldi 算法。该算法构造K k (A , r 0) 的一组标准正交基v 1, …, v k 和矩阵
h 11h 12h 13…h 1k
h 21
H k =
h 22h 32
h 23h 33
并将文[11]中块L anczo s 方法和块M I N R ES 方法
应用于求解对称线性方程组集(13) , 提出了求解非对称线性方程组集的块双对角化L anczo s 方法和块双对角化M I N R ES 方法。这两个方法具有计算量和存贮量小的优点。大量的数值试验[21]显示块双对角化L anczo s 方法和块双对角化M I N R ES 方法优于块G M R ES 方法。
1. 4 L anczo s M Q R 方法
[1]T
s K A , r 0) 和K k (A , r 0)
…ωω
h k , k -1
h 2k
h k -1,
k
h
11h 21
H k =
(e )
ω…
h 12h 22ω0
h 13h 23h h kk h 1k ……h k , k -1
, 即非对称L anczo s 算法。K k (A , r 0) 和K k (A T , r
0) 的一组双正交v 1, …, v k 和w 1, …, w k , 并产生矩阵
Α1
Β2Α2ω
ω…Β2Α2ωω……
0h k -1, k
…ωω
ω∆k 0
00
∆2
T k =
ωωω0
h kk 0
…00h k +1, 记V k =[v 1, …, v k ],则
(e ) T
AV k =V k H k +h k +1, k v k +1e k =V k +1H k (12) 取L k =K k (A , r 0) , W k =V k , 则W T k AV k =H k , 由OR 方法可产生完全正交化方法[12], 该方法所得近似
-1
解为x k =x 0+ΒV k H k e 1。
在式(4) 中取22范数, 则由M R 方法可产生
[13]
A xelsson 2L east 2Squares 算法、OR THOD I R 算
[15][16]
法[14]、GCR 算法和G M R ES 算法。理论上这4个算法是等价的。从数值稳定性、计算量和存贮量方面考虑, G M R ES 算法是值得推荐的。该方法利用式(12) 将最小二乘问题式(4) 化为‖r k ‖2=
(e )
m in k ‖Βe 1-H k y ‖2, 并由式(7) 确定近似解。理论
y ∈R
Β
k Αk
Α1∆2
T k =
(e )
…ωωω∆k 0
ωωω00
00
Βk Αk ∆k + 记V k =[v 1, …, v k ],W k =[w 1, …, w k ],则
(e ) δAV k =V k T k +[0, …, 0, v k +1]=V k +1T k
δ]T T
A W =W T +[0, …, 0, w
k
k
k
k +1
(14)
T
V k 和W k 满足关系W k V k =I , 并且sp an (V k ) =
上已证明对正实矩阵, G M R ES 算法收敛[16], 但其收敛率与谱分布情况密切相关。
由于A rno ldi 过程具有长递推关系, 其计算量和存贮量较大。目前已有G M R ES 算法的多个实用变形, 其中应用最广泛的是重新开始G M R ES 算法[16], 但重新开始G M R ES 算法可能收敛很慢。为了克服这些缺点, 文[17]提出了混合G M R ES 算法; 而文[18]提出了利用特征向量的重新开始
[19]
讨论了混合G M R ES 算法。最近, 钟宝江
G M R ES 算法的有效实现。
对非对称线性方程组集AX =B , V ital [20]利用块A rno ldi 过程, 提出了块G M R ES 方法。但是块
作者在文G M R ES 方法的计算量和存贮量都很大。
[21]=B 对称化为
0A Y B
(13) =T
0A X
K k (A , r 0) , span (W k ) =K k (A , r 0) 。
T
取L k =K k (A T , r 0) , 利用非对称
L anczo s 算法, 由OR 方法导出了求解式(1) 的L anczo s 双正交化算法。但仍存在两方面的问题:
Saad
[22]
一是如何有效地计算收敛性判据; 二是如何将非对称L anczo s 过程与求解非对称三对角方程组有机结合起来。文[23]给出了计算L anczo s 双正交化算法收敛性判据的一个有效方法, 并且利用三对角矩阵的LQ 分解和导出的收敛性判据将非对称L anczo s 过程与求解三对角方程组巧妙地结合起来, 提出了求解方程组(1) 的UN SY MM LQ 方法。
M R 方法求y 使得
()
‖r k ‖2=m in k ‖V k +1(Βe 1-T k e y ) ‖2(15)
y ∈R
从而得到近似解x k =x 0+V k y 。因为V k +1不是列正交规范矩阵, 求解最小二乘问题(15) 的计算量和存
142南 京 航 空 航 天 大 学 学 报第33卷
[24]
贮量都较大。F reund 和N ach tigal 建议求y 使得
()
(16) m in ‖Βe 1-T k e y ‖2
y ∈R
k
给出了求解(1) 的QM R 算法。
δj +1, w δj +1) =0, 非对称L anczo s 算法发如果(v 生中断。为了解决这个问题, Parlett 等[25]和[26]
F reund 等分别提出了两个前瞻策略。但文[25]的计算工作量较大, 并且不易推广到块格式。
[27]
F reund 和M alho tra 将QM R 算法推广到求解大型非对称线性方程组集AX =B , 给出了块QM R 算法。
似, 并用相应的R itz 向量作为A 的特征向量的近似。这就是求解对称矩阵特征值问题的L anczo s 方法。
给出了L anczo s 方法中R itz 值的一
种收敛性估计, 不过其证明有错误; Paige [5]重新研
Kan iel
[28]
究了这一问题, 不仅给出了R itz 值的收敛性估计, 而且给出了R itz 向量的收敛性估计; Saad [29]在Paige 工作的基础上对R Paige 更简单的估计。
, v v k 是相互正交的, 并且至多n 。但是在出现舍入误差的情况下, 计v 1, …, v k 会失去正交性。因此, 长期以来人们一直认为L anczo s 方法是数值不稳定的。直到1971年, Paige [5]通过仔细的舍入误差分析, 发现了失去正交性恰与近似特征值的精度提高有关。之后, 经过大量的理论分析与数值试验, 人们充分认识到L anczo s 方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一。然而, 对于A 的某一指定的特征值Κ, k 究竟取多大时T k 才有与其非常靠近的特征值? 一般说来, 这取决于A 的特征值的分布、Κ与A 的其他特征值的分离程度以及Κ在
通常对位于谱区间两端并且A 的谱区间内的位置。
分离较好的特征值Κ, 当k νn 时T k 的特征值中就含有Κ的很好近似值。
在实际计算中, 若k 较大, 因{v i }的正交性逐渐失去, 常会发生A 的同一特征值被反复计算多次的现象。为了避免这一现象的产生, 文[30]提出了带选择正交化的L anczo s 方法。当矩阵的阶数非常高时, 常取一适当的k , 迭代使用L anczo s 方法(即重新开始L anczo s 方法) 。但其收敛速度与初始向量的选取有关, 文[31]利用Chebyshev 多项式给出了加速重新开始L anczo s 方法的一个技巧, 提出了迭代Chebyshev 2L anczo s 方法; 文[32]提出了隐式重新开始L anczo s 方法。
当要求的特征值密集或其中有重特征值时, L anczo s 方法及其变形的有效性与可靠性会下降。为了克服这个缺点, U nder w ood [33]分别提出了求解对称矩阵特征值问题的块L anczo s 方法, 并且文[33, 29]研究了块L anczo s 方法中R itz 对的收敛性, 其收敛率与初始向量有关。文[34]给出了选取初始向量的一种策略, 提出了块Chebyshev 2L anczo s 方法以加速块L anczo s 方法的收敛性。目前, 该方法是求解大型稀疏对称矩阵部分极端特征对的最有效方法。
2 矩阵特征值问题的K 方法
对矩阵特征值问题
A x =Κx
(17)
其中A 是n 阶矩阵, 记K 是C n 中的一个k 维子空
间。对子空间K 的正交投影方法就是求问题(17)
～υ使得的近似特征向量
k υC 。Κ和y 是如下k 维特征值问题的特征对
H υy (19) V k AV k y =Κ
～
～
υ为R itz 值,
对子空间K 的正交投影方法就是用R itz 对υ, ～(Κ, x ) 。随着子空间K 的维
～υ～
数增大, 残量r =A
υ, ～
上, 如果残量范数充分小, 通常取R itz 对(Κ
子空间K 的最佳选择之一是由A 和v 1产生的K rylov 子空间, 即K =K k (A , v 1) =sp an (v 1,
k -1
A v 1, …, A v 1) , 由此导出的方法称为K rylov 子空间方法。从稳定性和有效性考虑, 需要构造K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 并且使k 阶矩阵V H k AV k 的特征值问题(19) 容易求解。
2. 1 对称矩阵的L anczo s 方法
对对称(或H er m ite ) 矩阵A , 利用对称L anczo s 算法可产生K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 并且V H k AV k =T k 为对称三对角矩阵。计算T k 的特征值作为A 的特征值的近
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法143
2. 2 非对称矩阵的A rno ldi 方法
对非对称矩阵A , 用A rno ldi 算法可产生K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 且V H 。求解k AV k =H k 为上H essenberg 矩阵
特征值问题(17) 的A rno ldi 方法[3]就是求H k 的特征值作为A 的近似特征值, 并用相应的R itz 向量作为A 的近似特征向量。
用A rno ldi 算法构造K k (A , v 1) 的标准正交基
v 1, …, v k 时, 需要把所有基向量都存在主存中, 当k 较大时, A rno ldi 方法的存贮量很大甚至无法实
接近于A -1的M -1并且容易计算M -1。如果选取好的预条件矩阵, 各种K rylov 子空间算法的计算效率没有多大差别, 但如何识别并构造好的预条件矩阵是一个重要的研究课题。在这个领域已取得了一些进展[9], 但对高度非正规矩阵和对称不定矩阵, 如何选取合适的预条件矩阵仍有待进一步研究。
用K rylov 子空间方法求大型线性方程组近似。但除, 、正实矩阵, 给出理论误差, 有待进一步研究。
[40]
K rylov 子空间方法。D avidson 方法和广义
[41]
D avidson 方法都可以看作求解对称矩阵特征值
) =(M -问题的预条件方法, 这些方法用算子N (Θ-1ΘI ) (A -ΘI ) 产生投影子空间, 其中Θ是要求特征值的近似, M 是对称矩阵A 的近似。M -ΘI 可看作
) 变换为对对A -ΘA 的预条件矩阵。文[42]将N (Θ
-T T
) =L -1(A -Θ称矩阵H (ΘI ) L , 其中LL 是预条
) 应用件矩阵M -ΘI 的Cho lesky 分解, 并对H (ΘL anczo s 方法, 从而给出了预条件L anczo s 方法。
现。为了保证快速收敛, 理论上k 取越大越好; 上从存贮量考虑又不能这样做。力取最大可能的k , ldi 方法计算的R itz 。R itz 向量或几个R itz 向量的线性组合代替v 1重复使用A rno ldi 方法, 这个过程称为迭代A rno ldi 方法或重新开始A rno ldi 方法。重新开始A rno ldi 方法的收敛速度可能很慢, 文[35]提出利用Chebyshev 多项式加速重新开始A rno ldi 方法, 但至今尚无有效方法确定Chebyshev 迭代中的参数。因此, 如何利用已积累的信息逐次重新构造更好的投影子空间是提高重新开始A rno ldi 方法有效性的关键。文[36]提出了一种隐式重新开始A rno ldi 方法。
对大型非对称矩阵A , 当所要求的特征值密集时, 为了使A rno ldi 方法收敛, K k (A , v 1) 的维数k 会大到存贮量和计算量令人无法接受的程度; 另一方面, 当A 有重特征值时, 如何确定其重数及其对应的特征空间是数值计算中的一大难题, A rno ldi 方法在理论和算法上不能解决这一难题。为此, 作者提出了块A rno ldi 方法[37], 该方法不仅可以求密集特征值, 而且较A rno ldi 方法具有更高的可靠性和有效性。之后, 文[38]进一步研究了块A rno ldi 方法, 并且文[39]讨论了如何利用Chebyshev 多项式加速块A rno ldi 方法。
在此基础上, 作者提出了预条件块L anczo s 方法[43]。对非对称矩阵, 如何应用预条件技术提高K rylov 子空间方法的收敛速度有待探讨。
参
考
文
献
1 L anczo s C . A n iterati on m ethod fo r the so luti on of the
eigenvalue p roblem of linear differential and integral operato rs [J ]. J R es N at Bur Stand , 1950, 45:255～282
2 L anczo s C . So luti on of system s of linear equati ons by m ini m ized iterati ons [J ]. J R es N at Bur Stand , 1952, 49:33～53
3 A rno ldiW E . T he p rinci p le of m ini m ized iterati ons in
the so luti on of the m atrix eigenvalue p roblem [J ].
～29Q uart A pp lM ath , 1951, 9:17
4 H estenes M R , Stiefel E . M ethods of conjugate
gradients fo r so lving linear system s [J ]. J R es N at Bur
3 讨 论
K rylov 子空间方法是求解大型稀疏线性方程
～436Stand , 1952, 49:409
5 Paige C C .
T he computati on of eigenvalues and
eigenvecto rs of very large sparse m atrices [D ]:[dissertati on ]. L ondon :U niversity of L ondon , 19716 Saad Y . V ariati ons on A rno ldi ′s m ethod fo r computing
eigenelem ents of large unsymm etric m atrices [J ].
组的一类有效算法, 其收敛速度依赖于矩阵特征值
或奇异值的分布。改善特征值或奇异值分布的方法称为预条件技术, 即选取预条件矩阵M 使M -1A 尽可能接近于单位矩阵。通常的做法是令M 在某种意义下接近A 并且M -1的计算易于实现或选取
～295L inear A lgebra A pp l , 1980, 34:269
7 Saad Y ,
Schultz M
H .
Conjugate gradient 2like
144南 京 航 空 航 天 大 学 学 报第33卷
algo rithm s fo r so lving nonsymm etric linear system s [J ]. M ath Comp , 1985, 44:417～4248 Saad Y .
Iterative m ethods fo r sparse linear system s
[M ]. P W S Publish ing Company , 1996
9 Paige C C , SaundersM A . So luti on of sparse indefinite
system s of linear equati ons [J ]. S I AM J N um er A nal , 1975, 12:617～629
10 O ′L eary D P . T he block conjugate gradient algo rithm
and related m ethods [J ]. L inear A lgebra A pp l , 1980, 29:293～32211 D ai H ua .
Tw o algo rithm s fo r symm etric linear
system s w ith m ulti p le righ t 2hand sides [J ]. N um er
unsymm etric system s [J ]. S I AM J N um er A nal , 1982, 19:485～506
23 周树荃, 戴 华. 求解大型非对称线性方程组的
UN SY MM LQ 方法[J ]. 高等学校计算数学学报, 1989, 10(2) :118～130
24 Parlett B N , T aylo r D R , L iu Z A . A look 2ahead
L anczo s algo rithm fo r unsymm etric m atrices [J ].
～124M ath Comp , 1985, 44:105
25 F reund R W , Gutknech t , N ach tigal N M . A n
i m p lem on L anczo s algo rithm r 2H er m ].
S I AM J Sci Stat
158, R , N ach tigal N M . QM R :a quasi 2m ini m al
residual m ethod fo r non 2H er m itian linear system s [J ].
～M ath 2A J of Ch inese U niversities , 2000, 9:91
12 Saad Y . K rylov subspace m ethods r unsymm etric linear J M , 37:105～126
13 A xelsson O . gradient type m ethods fo r
unsymm etric and
inconsistent
system s of
linear
～339N um er M ath , 1991, 60:315
27 F reund R W , M alho tra M . A block QM R algo rithm
fo r non 2H er m itian linear system s w ith m ulti p le righ t 2
～hand sides [J ]. L inear A lgebra A pp l , 1997, 254:119
157
28 Kaniel S . E sti m ates fo r som e computati onal tech 2
niques in linear algebra [J ]. M ath Comp , 1966, 20:369
～16equati ons [J ]. L inear A lgebra A pp l , 1980, 29:1
14 Young D M , Jea K C . Generalized conjugate 2gradient
accelerati on of nonsymm etriable m ethods [J ]. L inear
～194A lgebra A pp l , 1980, 34:159
15 E isentat S C , E l m an H C , Schultz M H . V ariati onal
iterative m ethods fo r nonsymm etric system s of linear equati ons [J ]. S I AM J N um er A nal , 1983, 20:345～357
16 Saad Y , Schultz M H .
G M R ES :a generalized
m ini m al residual algo rithm fo r so lving nonsymm etric linear system s [J ]. S I AM J Sci Stat Comput , 1986, 7:856～869
17 N ach tigal N M , R eichel L , T refethen L N . A hybrid
G M R ES algo rithm fo r nonsymm etric linear system s [J ]. S I ～825AM J M atrix A nal A pp l , 1992, 13:79618 M o rgan R B . A restarted G M R ES m ethod augm ented
w ith eigenvecto rs [J ]. S I AM J M atrix A nal A pp l , 1995, 16:1154～1171
19 Zhong Bao jiang . I mp lem etati on of the hybrid G M R ES
algo rithm using Househo lder transfo r m ati ons [J ]. T ransacti ons of N anjing U niversity of A eronatuics &
～378
29 Saad Y . O n the rates of convergence of the L anczo s
and the block L anczo s m ethods [J ]. S I AM J N um er
～706A nal , 1980, 17:689
30 Parlett B N , Sco tt D S . T he L anczo s algo rithm w ith
selective o rthogonalizati on [J ]. M ath Comp , 1979, 33:217～238
31 戴 华, 周树荃. 求解大型对称特征值问题的迭代
Chebyshev 2L anczo s 方法[J ]. 南京航空航天大学学
报, 1986, 18(4) :25～34
32 Calvetti D , R eichel L , So rensen D C . A n i m p licitly
restarted L anczo s m ethod
fo r
large
symm etric
～21eigenvalue p roblem s [J ]. ETNA , 1994, 2:1
33 U nder w ood R R . A n iterative block L anczo s m ethod
fo r
the
so luti on
of
large
sparse
symm etric Califo rnia :
eigenp roblem s [D ]:
[dissertati on ].
Stanfo rd U niversity , 1975
34 周树荃, 戴 华. 求解大型对称特征值问题的块
Chebyshev 2L anczo s 方法[J ]. 南京航空航天大学学
～152A stronautics , 1997, 14:146
20 V ital B . E tude de quelques m ethodes de reso luti on de
p roblem es
lineaires
de
grande
taille
sur
m ulti p rocesseur [D ]:
[dissertati on ].
R ennes :
报, 1989, 21(4) :22～28
35 Ho D , Chatelin F , Bennani M . A rno ldi 2Chebyshev
m ethod fo r large scale nonsymm etric m atrices [J ]. RA I RO M ath M odelling and N um er A nal , 1990, 24:53～65
36 So rensen D C .
I mp licit app licati on of po lynom ial
filters in a k 2step A rno ldim ethod [J ]. S I AM J M atrix
U niversite de R ennes , 1990
21 D ai H ua . B lock bidiagonalizati on m ethods fo r m ulti p le
nonsymm etric linear system s [C ]. CER FA CS , TR 98 35, 1998PA
22 Saad Y . T he L anczo s bi o rthogonalizati on algo rithm
and o ther oblique p ro jecti on m ethods fo r so lving large
～385A nal A pp l , 1992, 13:357
37 戴 华. 求解大型非对称特征问题的块A rno ldi 方法
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法
1975, 17:87～94
145
[J ]. 南京航空航天大学学报, 1987, 19(1) :100～10838 Sadkane M . B lock 2A rno ldi and D avidson m ethods fo r
unsymm etric large eigenvalue p roblem s [J ]. N um er
41 M o rgan R B , Sco tt D S . Generalizati on of D avidson ′s
m ethod
fo r
computing
eigenvalues
of
sparse
symm etric m atrices [J ]. S I AM J Sci Stat Comput , 1986, 7:817～825
42 M o rgan R B , Sco tt D S . P reconditi oning the L anczo s
algo rithm fo r sparse symm etric eigenvalue p roblem s [J ]. S I ～593AM J Sci Computing , 1993, 14:58543 D ai H ua , L ancaster P . P oning block L anczo s
algo fo r so symm etric eigenvalue
～211M ath , 1993, 64:195
39 Sadkane M . A block A rno ldi 2Chebyshev m ethod fo r
computing the leading eigenpairs of large sparse unsymm etric m atrices [J ]. N um er M ath , 1993, 64:181
～193
40 D avidson E R . T he iterative calculati on of a few of
the low est eigenvalues and co rresponding eigenvecto rs of large real symm etric m atrices [J ]. J Comp Phys ,
～374s [J ]. , 18:365
Krylov Subspace for Large Sca le M a tr ix Problem s
D a i H ua
(Co llege of , N anjing U niversity of A eronautics &A stronautics N anjing 210016, P . R . Ch ina )
Abstract So lving large scale m atrix p rob lem s including large linear system s , eigenvalue p rob lem s , etc . , is a vital sub ject in com pu tati onal m athem atics and scien tific engineering com pu ting . T here have been i m po rtan t advances in the sub ject in recen t years . A com p rehen sive su rvey of som e developm en ts including the au tho rs ′w o rk s regarding K rylov sub sp ace m ethods fo r so lving large linear system s and eigenvalue p rob lem s is given .
Sp ecific top ics include :con jugate gradien t algo rithm ,
SY MM LQ
algo rithm , M I N R ES algo rithm , G M R ES algo rithm , L anczo s b i o rthogonalizati on algo rithm , QM R algo rithm and their b lock versi on s fo r so lving large linear system s , L anczo s algo rithm and its b lock versi on fo r so lving large symm etric eigenvalue p rob lem s , L anczo s algo rithm , A rno ldi algo rithm and their b lock versi on s fo r so lving large un symm etric eigenvalue p rob lem s . T he accelerati on techn iques and p reconditi on ing techn iques fo r large scale m atrix p rob lem s are discu ssed . Som e p rob lem s that need to be fu rther studied are p resen ted .
Key words :m atrices ; linear system s ; eigenvalue ; K rylov sub sp ace m ethods
第33卷第2期南 京 航 空 航 天 大 学 学 报V o l . 33N o. 2
2001年4月 A p r . 2001Jou rnal of N an jing U n iversity of A eronau tics &A stronau tics
文章编号:100522615(2001) 0220139207
求解大规模矩阵问题的K rylov 子空间方法
戴 华
(南京航空航天大学理学院 南京, 210016)
摘要 。最近几年, 其研究工作取得了许多重大进展。K rylov 子空间方法若干进展的一个概述, 。、SY MM LQ 算法、M I N R ES 算法、G M s 、QM R 算法以及这些算法的块格式; 求解大型对称特征值问题的L s L s 算法; 求解大型非对称特征值问题的L anczo s 算法、A rno ldi 算法以及。提出了一些有待进一步研究的问题。关键词:矩阵; 线性方程组; 特征值; K rylov 子空间方法中图分类号:O 241. 6 文献标识码:A
引 言
大规模矩阵问题包括线性方程组、矩阵特征值问题等数值解法的理论研究、算法设计和软件研制是当今计算数学和科学工程计算研究的重大课题, 是大规模科学与工程计算的基础和重要组成部分, 其研究具有重要的理论意义和广泛的应用价值。尤其是其中大型非对称矩阵问题的数值计算, 因其理论上的复杂性, 导致在算法的理论研究和软件设计中产生了许多本质的困难。问题的挑战性吸引了国内外许多数值分析专家从事该领域的研究, 并取得了许多重大进展。
求解线性方程组和矩阵特征值问题的K rylov 子空间方法可以追溯到50年代初L anczo s
[1, 2]
理论分析和数值试验, 充分认识到K rylov 子空间
方法是求解大型线性方程组和大型矩阵特征值问题的一类最有效的方法。本文概述K rylov 子空间方法的若干进展, 其中包括作者对这些问题的研究成果, 并提出一些有待进一步研究的问题。
1 求解大型线性方程组的K rylov 子
空间方法
对线性方程组
A x =b
(1)
其中, A ∈R n ×n , b ∈R n , 设x 0∈R n 是式(1) 解的初始估计, r 0=b -A x 0是初始残量。称K k (A , r 0) =
k -span {r 0, A r 0, …, A
1
r 0}为由A 和r 0产生的
,
K rylov 子空间。K rylov 子空间方法求式(1) 的具有
[3][4]
A rno ldi , H estenes 和Stiefel 的工作。但在其后
20年, 人们一直认为这类方法是数值不稳定的, 很
如下形式的近似解
x k ∈x 0+K k (A , r 0) 并且使残量r k =b -A x k 满足Galerk in 条件
r k ⊥L k
(2) (3)
少用于实际计算。直到70年代初Paige 在其著名的博士论文[5]中重新研究了L anczo s 方法以及80年代初Saad [6]重新研究了A rno ldi 方法, 才使人们重新认识了这两个方法。之后, 人们又做了大量的
其中L k 是R n 中的一个子空间, 或使残量r k 的范数极小化, 即
基金项目:国家自然科学基金(编号:19671043) 、江苏省自然科学基金(编号:BK 97059) 、江苏省“333工程”基金和江苏省“青蓝工程”基金资助项目。
收稿日期:2000204211; 修订日期:2000211201
作者简介:戴 华, 男, 教授, 博士导师, 1959年6月生。
140南 京 航 空 航 天 大 学 学 报第33卷
‖r k ‖=
x ∈x 0+K k (A , r 0)
m in ‖b -A x ‖(4)
T k =
Α1Β20
Β2Α2ω
ω…Β2Α2ω0…
ωωωΒk …ωωωΒk 0
T
用式(2, 3) 定义近似解x k 的K rylov 子空间方法称为正交残量(OR ) 方法。用V k 和W k 分别表示子空间K k (A , r 0) 和L k 的基为列作成的矩阵, OR 方法可行的充要条件是W T 。OR 方法产生k AV k 非奇异的近似解为
x k =x 0+V k (W k AV k ) W k r 0
-1
T
ωω
ω00
Βk Αk 0
Α1Β2
T k =
(e )
(5)
ωω0
用式(2, 4) 定义近似解x k 的K rylov 子空间方法称为极小残量(M R ) 方法。如果最小二乘问题
m in k ‖r 0-AV k y ‖
y ∈R
Βk Αk Βk +(e )
(6)
的解为y , 则M R 方法产生的近似解为
x k =x 0+V k y
k =[v 1, …, v k ],则
AV k =V k T k +Βk +1v k +1e k =V k +1T k
T
如果在式(3) A , r , 方[7]。
在式(3) L k 或在式(4) 中选取不同的范数可导出不同的算法。
1. 1 对称正定线性方程组的共轭梯度法
1952年H estenes 和Stiefel [4]提出了求解对称
(11)
SY MM LQ 算法是在OR 方法中取L k =
K k (A , r 0) , W k =V k , 则W k AV k =T k , x k =x 0+
-1
ΒV k T k e 1; M I N R ES 算法是在式(10) 意义下的M R
方法, 由式(11) , 式(6) 化为m in k ‖r 0-AV k y ‖2=
y ∈R
m in k ‖Βe 1-T k y ‖2, 则容易求得y , 并由式(7) 即
y ∈R
(e )
正定线性方程组的共轭梯度法(CG ) 。该方法可以充分利用矩阵A 的稀疏性, 每次迭代的主要计算量是向量之间的运算。共轭梯度法是一个K rylov 子空间方法。实际上, 共轭梯度法是在向量范数
2
‖x ‖A =(A x , x ) 1 下的M R 方法[7], 近似解x k 满足如下极小化性质‖b -A x k ‖A
-1
得近似解x k 。
SY MM LQ 算法和M I N R ES 算法均以三项递
推公式进行计算, 具有计算量小、存贮量少的优点, 并且M I N R ES 算法的残量范数单调下降。这两个
算法至今仍是求解大型对称不定线性方程组的最有效方法。
对具有相同系数矩阵而具有不同右端项的大型线性方程组集A x i =b i (i =1, …, m ) , 对所有方程组AX =B (其中X =[x 1, …, x m ],B =[b 1, …, b m ]) 采用块方法同时进行迭代比对各个方程组A x i =b i
[10]
使用通常的迭代法更有效。O ′L eary 提出了求解
=
x ∈x 0+K k (A , r 0)
m in ‖b -A x ‖A
-1
(8)
并且[8]
‖x -x k ‖A ≤2‖x -x 0‖A
+k
(9)
大型对称线性方程组集的块CG 方法, 并利用块L anczo s 过程
[33]
其中ϑ=K (A ) =Κm ax (A ) m in (A ) 表示A 的谱条件Κ数。
1. 2 对称不定线性方程组的SY MM LQ 算法和
M I N R ES 算法
将SY MM LQ 算法推广到求解AX
=B , 给出了块L anczo s 方法。但是对所有不同的右
端项由块L anczo s 方法产生的序列收敛速度一般是不同的, 并且残量范数未必单调下降。一个有效的块方法应该能判断、自适应地收缩收敛的方程组, 但块L anczo s 方法的收缩问题一直未解决。作者[11]利用收敛性判据识别和收缩收敛的方程组, 提出了求解大型对称线性方程组集的自适应块L anczo s 方法。为了解决块L anczo s 方法中残量范
对对称不定或非对称矩阵A , 人们通常用如下极小化性质代替式(8)
‖b -A x k ‖2=
x ∈x 0+K k (A , r 0)
m in
‖b -A x ‖2(10)
1975年Paige 和Saunders [9]提出了求解对称不定线性方程组的数值稳定的SY MM LQ 算法和M I N R ES 算法。这两个算法利用对称L anczo s 过
数的振荡, 文[11]提出了M I N R ES 方法的一个块
版本——块M I N R ES 方法, 描述了块M I N R ES 方法的收缩技术, 建立了块L anczo s 方法与块M I N R ES 方法之间的关系。
程[1]产生K k (A , r 0) 的一组标准正交基v 1, …, v k 和
矩阵
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法141
1. 3 G M R ES 及相关算法
对非对称矩阵A , A rno ldi [3]提出了产生K rylov 子空间K k (A , r 0) 的一组标准正交基的方法, 即A rno ldi 算法。该算法构造K k (A , r 0) 的一组标准正交基v 1, …, v k 和矩阵
h 11h 12h 13…h 1k
h 21
H k =
h 22h 32
h 23h 33
并将文[11]中块L anczo s 方法和块M I N R ES 方法
应用于求解对称线性方程组集(13) , 提出了求解非对称线性方程组集的块双对角化L anczo s 方法和块双对角化M I N R ES 方法。这两个方法具有计算量和存贮量小的优点。大量的数值试验[21]显示块双对角化L anczo s 方法和块双对角化M I N R ES 方法优于块G M R ES 方法。
1. 4 L anczo s M Q R 方法
[1]T
s K A , r 0) 和K k (A , r 0)
…ωω
h k , k -1
h 2k
h k -1,
k
h
11h 21
H k =
(e )
ω…
h 12h 22ω0
h 13h 23h h kk h 1k ……h k , k -1
, 即非对称L anczo s 算法。K k (A , r 0) 和K k (A T , r
0) 的一组双正交v 1, …, v k 和w 1, …, w k , 并产生矩阵
Α1
Β2Α2ω
ω…Β2Α2ωω……
0h k -1, k
…ωω
ω∆k 0
00
∆2
T k =
ωωω0
h kk 0
…00h k +1, 记V k =[v 1, …, v k ],则
(e ) T
AV k =V k H k +h k +1, k v k +1e k =V k +1H k (12) 取L k =K k (A , r 0) , W k =V k , 则W T k AV k =H k , 由OR 方法可产生完全正交化方法[12], 该方法所得近似
-1
解为x k =x 0+ΒV k H k e 1。
在式(4) 中取22范数, 则由M R 方法可产生
[13]
A xelsson 2L east 2Squares 算法、OR THOD I R 算
[15][16]
法[14]、GCR 算法和G M R ES 算法。理论上这4个算法是等价的。从数值稳定性、计算量和存贮量方面考虑, G M R ES 算法是值得推荐的。该方法利用式(12) 将最小二乘问题式(4) 化为‖r k ‖2=
(e )
m in k ‖Βe 1-H k y ‖2, 并由式(7) 确定近似解。理论
y ∈R
Β
k Αk
Α1∆2
T k =
(e )
…ωωω∆k 0
ωωω00
00
Βk Αk ∆k + 记V k =[v 1, …, v k ],W k =[w 1, …, w k ],则
(e ) δAV k =V k T k +[0, …, 0, v k +1]=V k +1T k
δ]T T
A W =W T +[0, …, 0, w
k
k
k
k +1
(14)
T
V k 和W k 满足关系W k V k =I , 并且sp an (V k ) =
上已证明对正实矩阵, G M R ES 算法收敛[16], 但其收敛率与谱分布情况密切相关。
由于A rno ldi 过程具有长递推关系, 其计算量和存贮量较大。目前已有G M R ES 算法的多个实用变形, 其中应用最广泛的是重新开始G M R ES 算法[16], 但重新开始G M R ES 算法可能收敛很慢。为了克服这些缺点, 文[17]提出了混合G M R ES 算法; 而文[18]提出了利用特征向量的重新开始
[19]
讨论了混合G M R ES 算法。最近, 钟宝江
G M R ES 算法的有效实现。
对非对称线性方程组集AX =B , V ital [20]利用块A rno ldi 过程, 提出了块G M R ES 方法。但是块
作者在文G M R ES 方法的计算量和存贮量都很大。
[21]=B 对称化为
0A Y B
(13) =T
0A X
K k (A , r 0) , span (W k ) =K k (A , r 0) 。
T
取L k =K k (A T , r 0) , 利用非对称
L anczo s 算法, 由OR 方法导出了求解式(1) 的L anczo s 双正交化算法。但仍存在两方面的问题:
Saad
[22]
一是如何有效地计算收敛性判据; 二是如何将非对称L anczo s 过程与求解非对称三对角方程组有机结合起来。文[23]给出了计算L anczo s 双正交化算法收敛性判据的一个有效方法, 并且利用三对角矩阵的LQ 分解和导出的收敛性判据将非对称L anczo s 过程与求解三对角方程组巧妙地结合起来, 提出了求解方程组(1) 的UN SY MM LQ 方法。
M R 方法求y 使得
()
‖r k ‖2=m in k ‖V k +1(Βe 1-T k e y ) ‖2(15)
y ∈R
从而得到近似解x k =x 0+V k y 。因为V k +1不是列正交规范矩阵, 求解最小二乘问题(15) 的计算量和存
142南 京 航 空 航 天 大 学 学 报第33卷
[24]
贮量都较大。F reund 和N ach tigal 建议求y 使得
()
(16) m in ‖Βe 1-T k e y ‖2
y ∈R
k
给出了求解(1) 的QM R 算法。
δj +1, w δj +1) =0, 非对称L anczo s 算法发如果(v 生中断。为了解决这个问题, Parlett 等[25]和[26]
F reund 等分别提出了两个前瞻策略。但文[25]的计算工作量较大, 并且不易推广到块格式。
[27]
F reund 和M alho tra 将QM R 算法推广到求解大型非对称线性方程组集AX =B , 给出了块QM R 算法。
似, 并用相应的R itz 向量作为A 的特征向量的近似。这就是求解对称矩阵特征值问题的L anczo s 方法。
给出了L anczo s 方法中R itz 值的一
种收敛性估计, 不过其证明有错误; Paige [5]重新研
Kan iel
[28]
究了这一问题, 不仅给出了R itz 值的收敛性估计, 而且给出了R itz 向量的收敛性估计; Saad [29]在Paige 工作的基础上对R Paige 更简单的估计。
, v v k 是相互正交的, 并且至多n 。但是在出现舍入误差的情况下, 计v 1, …, v k 会失去正交性。因此, 长期以来人们一直认为L anczo s 方法是数值不稳定的。直到1971年, Paige [5]通过仔细的舍入误差分析, 发现了失去正交性恰与近似特征值的精度提高有关。之后, 经过大量的理论分析与数值试验, 人们充分认识到L anczo s 方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一。然而, 对于A 的某一指定的特征值Κ, k 究竟取多大时T k 才有与其非常靠近的特征值? 一般说来, 这取决于A 的特征值的分布、Κ与A 的其他特征值的分离程度以及Κ在
通常对位于谱区间两端并且A 的谱区间内的位置。
分离较好的特征值Κ, 当k νn 时T k 的特征值中就含有Κ的很好近似值。
在实际计算中, 若k 较大, 因{v i }的正交性逐渐失去, 常会发生A 的同一特征值被反复计算多次的现象。为了避免这一现象的产生, 文[30]提出了带选择正交化的L anczo s 方法。当矩阵的阶数非常高时, 常取一适当的k , 迭代使用L anczo s 方法(即重新开始L anczo s 方法) 。但其收敛速度与初始向量的选取有关, 文[31]利用Chebyshev 多项式给出了加速重新开始L anczo s 方法的一个技巧, 提出了迭代Chebyshev 2L anczo s 方法; 文[32]提出了隐式重新开始L anczo s 方法。
当要求的特征值密集或其中有重特征值时, L anczo s 方法及其变形的有效性与可靠性会下降。为了克服这个缺点, U nder w ood [33]分别提出了求解对称矩阵特征值问题的块L anczo s 方法, 并且文[33, 29]研究了块L anczo s 方法中R itz 对的收敛性, 其收敛率与初始向量有关。文[34]给出了选取初始向量的一种策略, 提出了块Chebyshev 2L anczo s 方法以加速块L anczo s 方法的收敛性。目前, 该方法是求解大型稀疏对称矩阵部分极端特征对的最有效方法。
2 矩阵特征值问题的K 方法
对矩阵特征值问题
A x =Κx
(17)
其中A 是n 阶矩阵, 记K 是C n 中的一个k 维子空
间。对子空间K 的正交投影方法就是求问题(17)
～υ使得的近似特征向量
k υC 。Κ和y 是如下k 维特征值问题的特征对
H υy (19) V k AV k y =Κ
～
～
υ为R itz 值,
对子空间K 的正交投影方法就是用R itz 对υ, ～(Κ, x ) 。随着子空间K 的维
～υ～
数增大, 残量r =A
υ, ～
上, 如果残量范数充分小, 通常取R itz 对(Κ
子空间K 的最佳选择之一是由A 和v 1产生的K rylov 子空间, 即K =K k (A , v 1) =sp an (v 1,
k -1
A v 1, …, A v 1) , 由此导出的方法称为K rylov 子空间方法。从稳定性和有效性考虑, 需要构造K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 并且使k 阶矩阵V H k AV k 的特征值问题(19) 容易求解。
2. 1 对称矩阵的L anczo s 方法
对对称(或H er m ite ) 矩阵A , 利用对称L anczo s 算法可产生K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 并且V H k AV k =T k 为对称三对角矩阵。计算T k 的特征值作为A 的特征值的近
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法143
2. 2 非对称矩阵的A rno ldi 方法
对非对称矩阵A , 用A rno ldi 算法可产生K rylov 子空间K k (A , v 1) 的一组标准正交基v 1, …, v k , 且V H 。求解k AV k =H k 为上H essenberg 矩阵
特征值问题(17) 的A rno ldi 方法[3]就是求H k 的特征值作为A 的近似特征值, 并用相应的R itz 向量作为A 的近似特征向量。
用A rno ldi 算法构造K k (A , v 1) 的标准正交基
v 1, …, v k 时, 需要把所有基向量都存在主存中, 当k 较大时, A rno ldi 方法的存贮量很大甚至无法实
接近于A -1的M -1并且容易计算M -1。如果选取好的预条件矩阵, 各种K rylov 子空间算法的计算效率没有多大差别, 但如何识别并构造好的预条件矩阵是一个重要的研究课题。在这个领域已取得了一些进展[9], 但对高度非正规矩阵和对称不定矩阵, 如何选取合适的预条件矩阵仍有待进一步研究。
用K rylov 子空间方法求大型线性方程组近似。但除, 、正实矩阵, 给出理论误差, 有待进一步研究。
[40]
K rylov 子空间方法。D avidson 方法和广义
[41]
D avidson 方法都可以看作求解对称矩阵特征值
) =(M -问题的预条件方法, 这些方法用算子N (Θ-1ΘI ) (A -ΘI ) 产生投影子空间, 其中Θ是要求特征值的近似, M 是对称矩阵A 的近似。M -ΘI 可看作
) 变换为对对A -ΘA 的预条件矩阵。文[42]将N (Θ
-T T
) =L -1(A -Θ称矩阵H (ΘI ) L , 其中LL 是预条
) 应用件矩阵M -ΘI 的Cho lesky 分解, 并对H (ΘL anczo s 方法, 从而给出了预条件L anczo s 方法。
现。为了保证快速收敛, 理论上k 取越大越好; 上从存贮量考虑又不能这样做。力取最大可能的k , ldi 方法计算的R itz 。R itz 向量或几个R itz 向量的线性组合代替v 1重复使用A rno ldi 方法, 这个过程称为迭代A rno ldi 方法或重新开始A rno ldi 方法。重新开始A rno ldi 方法的收敛速度可能很慢, 文[35]提出利用Chebyshev 多项式加速重新开始A rno ldi 方法, 但至今尚无有效方法确定Chebyshev 迭代中的参数。因此, 如何利用已积累的信息逐次重新构造更好的投影子空间是提高重新开始A rno ldi 方法有效性的关键。文[36]提出了一种隐式重新开始A rno ldi 方法。
对大型非对称矩阵A , 当所要求的特征值密集时, 为了使A rno ldi 方法收敛, K k (A , v 1) 的维数k 会大到存贮量和计算量令人无法接受的程度; 另一方面, 当A 有重特征值时, 如何确定其重数及其对应的特征空间是数值计算中的一大难题, A rno ldi 方法在理论和算法上不能解决这一难题。为此, 作者提出了块A rno ldi 方法[37], 该方法不仅可以求密集特征值, 而且较A rno ldi 方法具有更高的可靠性和有效性。之后, 文[38]进一步研究了块A rno ldi 方法, 并且文[39]讨论了如何利用Chebyshev 多项式加速块A rno ldi 方法。
在此基础上, 作者提出了预条件块L anczo s 方法[43]。对非对称矩阵, 如何应用预条件技术提高K rylov 子空间方法的收敛速度有待探讨。
参
考
文
献
1 L anczo s C . A n iterati on m ethod fo r the so luti on of the
eigenvalue p roblem of linear differential and integral operato rs [J ]. J R es N at Bur Stand , 1950, 45:255～282
2 L anczo s C . So luti on of system s of linear equati ons by m ini m ized iterati ons [J ]. J R es N at Bur Stand , 1952, 49:33～53
3 A rno ldiW E . T he p rinci p le of m ini m ized iterati ons in
the so luti on of the m atrix eigenvalue p roblem [J ].
～29Q uart A pp lM ath , 1951, 9:17
4 H estenes M R , Stiefel E . M ethods of conjugate
gradients fo r so lving linear system s [J ]. J R es N at Bur
3 讨 论
K rylov 子空间方法是求解大型稀疏线性方程
～436Stand , 1952, 49:409
5 Paige C C .
T he computati on of eigenvalues and
eigenvecto rs of very large sparse m atrices [D ]:[dissertati on ]. L ondon :U niversity of L ondon , 19716 Saad Y . V ariati ons on A rno ldi ′s m ethod fo r computing
eigenelem ents of large unsymm etric m atrices [J ].
组的一类有效算法, 其收敛速度依赖于矩阵特征值
或奇异值的分布。改善特征值或奇异值分布的方法称为预条件技术, 即选取预条件矩阵M 使M -1A 尽可能接近于单位矩阵。通常的做法是令M 在某种意义下接近A 并且M -1的计算易于实现或选取
～295L inear A lgebra A pp l , 1980, 34:269
7 Saad Y ,
Schultz M
H .
Conjugate gradient 2like
144南 京 航 空 航 天 大 学 学 报第33卷
algo rithm s fo r so lving nonsymm etric linear system s [J ]. M ath Comp , 1985, 44:417～4248 Saad Y .
Iterative m ethods fo r sparse linear system s
[M ]. P W S Publish ing Company , 1996
9 Paige C C , SaundersM A . So luti on of sparse indefinite
system s of linear equati ons [J ]. S I AM J N um er A nal , 1975, 12:617～629
10 O ′L eary D P . T he block conjugate gradient algo rithm
and related m ethods [J ]. L inear A lgebra A pp l , 1980, 29:293～32211 D ai H ua .
Tw o algo rithm s fo r symm etric linear
system s w ith m ulti p le righ t 2hand sides [J ]. N um er
unsymm etric system s [J ]. S I AM J N um er A nal , 1982, 19:485～506
23 周树荃, 戴 华. 求解大型非对称线性方程组的
UN SY MM LQ 方法[J ]. 高等学校计算数学学报, 1989, 10(2) :118～130
24 Parlett B N , T aylo r D R , L iu Z A . A look 2ahead
L anczo s algo rithm fo r unsymm etric m atrices [J ].
～124M ath Comp , 1985, 44:105
25 F reund R W , Gutknech t , N ach tigal N M . A n
i m p lem on L anczo s algo rithm r 2H er m ].
S I AM J Sci Stat
158, R , N ach tigal N M . QM R :a quasi 2m ini m al
residual m ethod fo r non 2H er m itian linear system s [J ].
～M ath 2A J of Ch inese U niversities , 2000, 9:91
12 Saad Y . K rylov subspace m ethods r unsymm etric linear J M , 37:105～126
13 A xelsson O . gradient type m ethods fo r
unsymm etric and
inconsistent
system s of
linear
～339N um er M ath , 1991, 60:315
27 F reund R W , M alho tra M . A block QM R algo rithm
fo r non 2H er m itian linear system s w ith m ulti p le righ t 2
～hand sides [J ]. L inear A lgebra A pp l , 1997, 254:119
157
28 Kaniel S . E sti m ates fo r som e computati onal tech 2
niques in linear algebra [J ]. M ath Comp , 1966, 20:369
～16equati ons [J ]. L inear A lgebra A pp l , 1980, 29:1
14 Young D M , Jea K C . Generalized conjugate 2gradient
accelerati on of nonsymm etriable m ethods [J ]. L inear
～194A lgebra A pp l , 1980, 34:159
15 E isentat S C , E l m an H C , Schultz M H . V ariati onal
iterative m ethods fo r nonsymm etric system s of linear equati ons [J ]. S I AM J N um er A nal , 1983, 20:345～357
16 Saad Y , Schultz M H .
G M R ES :a generalized
m ini m al residual algo rithm fo r so lving nonsymm etric linear system s [J ]. S I AM J Sci Stat Comput , 1986, 7:856～869
17 N ach tigal N M , R eichel L , T refethen L N . A hybrid
G M R ES algo rithm fo r nonsymm etric linear system s [J ]. S I ～825AM J M atrix A nal A pp l , 1992, 13:79618 M o rgan R B . A restarted G M R ES m ethod augm ented
w ith eigenvecto rs [J ]. S I AM J M atrix A nal A pp l , 1995, 16:1154～1171
19 Zhong Bao jiang . I mp lem etati on of the hybrid G M R ES
algo rithm using Househo lder transfo r m ati ons [J ]. T ransacti ons of N anjing U niversity of A eronatuics &
～378
29 Saad Y . O n the rates of convergence of the L anczo s
and the block L anczo s m ethods [J ]. S I AM J N um er
～706A nal , 1980, 17:689
30 Parlett B N , Sco tt D S . T he L anczo s algo rithm w ith
selective o rthogonalizati on [J ]. M ath Comp , 1979, 33:217～238
31 戴 华, 周树荃. 求解大型对称特征值问题的迭代
Chebyshev 2L anczo s 方法[J ]. 南京航空航天大学学
报, 1986, 18(4) :25～34
32 Calvetti D , R eichel L , So rensen D C . A n i m p licitly
restarted L anczo s m ethod
fo r
large
symm etric
～21eigenvalue p roblem s [J ]. ETNA , 1994, 2:1
33 U nder w ood R R . A n iterative block L anczo s m ethod
fo r
the
so luti on
of
large
sparse
symm etric Califo rnia :
eigenp roblem s [D ]:
[dissertati on ].
Stanfo rd U niversity , 1975
34 周树荃, 戴 华. 求解大型对称特征值问题的块
Chebyshev 2L anczo s 方法[J ]. 南京航空航天大学学
～152A stronautics , 1997, 14:146
20 V ital B . E tude de quelques m ethodes de reso luti on de
p roblem es
lineaires
de
grande
taille
sur
m ulti p rocesseur [D ]:
[dissertati on ].
R ennes :
报, 1989, 21(4) :22～28
35 Ho D , Chatelin F , Bennani M . A rno ldi 2Chebyshev
m ethod fo r large scale nonsymm etric m atrices [J ]. RA I RO M ath M odelling and N um er A nal , 1990, 24:53～65
36 So rensen D C .
I mp licit app licati on of po lynom ial
filters in a k 2step A rno ldim ethod [J ]. S I AM J M atrix
U niversite de R ennes , 1990
21 D ai H ua . B lock bidiagonalizati on m ethods fo r m ulti p le
nonsymm etric linear system s [C ]. CER FA CS , TR 98 35, 1998PA
22 Saad Y . T he L anczo s bi o rthogonalizati on algo rithm
and o ther oblique p ro jecti on m ethods fo r so lving large
～385A nal A pp l , 1992, 13:357
37 戴 华. 求解大型非对称特征问题的块A rno ldi 方法
第2期戴 华:求解大规模矩阵问题的K rylov 子空间方法
1975, 17:87～94
145
[J ]. 南京航空航天大学学报, 1987, 19(1) :100～10838 Sadkane M . B lock 2A rno ldi and D avidson m ethods fo r
unsymm etric large eigenvalue p roblem s [J ]. N um er
41 M o rgan R B , Sco tt D S . Generalizati on of D avidson ′s
m ethod
fo r
computing
eigenvalues
of
sparse
symm etric m atrices [J ]. S I AM J Sci Stat Comput , 1986, 7:817～825
42 M o rgan R B , Sco tt D S . P reconditi oning the L anczo s
algo rithm fo r sparse symm etric eigenvalue p roblem s [J ]. S I ～593AM J Sci Computing , 1993, 14:58543 D ai H ua , L ancaster P . P oning block L anczo s
algo fo r so symm etric eigenvalue
～211M ath , 1993, 64:195
39 Sadkane M . A block A rno ldi 2Chebyshev m ethod fo r
computing the leading eigenpairs of large sparse unsymm etric m atrices [J ]. N um er M ath , 1993, 64:181
～193
40 D avidson E R . T he iterative calculati on of a few of
the low est eigenvalues and co rresponding eigenvecto rs of large real symm etric m atrices [J ]. J Comp Phys ,
～374s [J ]. , 18:365
Krylov Subspace for Large Sca le M a tr ix Problem s
D a i H ua
(Co llege of , N anjing U niversity of A eronautics &A stronautics N anjing 210016, P . R . Ch ina )
Abstract So lving large scale m atrix p rob lem s including large linear system s , eigenvalue p rob lem s , etc . , is a vital sub ject in com pu tati onal m athem atics and scien tific engineering com pu ting . T here have been i m po rtan t advances in the sub ject in recen t years . A com p rehen sive su rvey of som e developm en ts including the au tho rs ′w o rk s regarding K rylov sub sp ace m ethods fo r so lving large linear system s and eigenvalue p rob lem s is given .
Sp ecific top ics include :con jugate gradien t algo rithm ,
SY MM LQ
algo rithm , M I N R ES algo rithm , G M R ES algo rithm , L anczo s b i o rthogonalizati on algo rithm , QM R algo rithm and their b lock versi on s fo r so lving large linear system s , L anczo s algo rithm and its b lock versi on fo r so lving large symm etric eigenvalue p rob lem s , L anczo s algo rithm , A rno ldi algo rithm and their b lock versi on s fo r so lving large un symm etric eigenvalue p rob lem s . T he accelerati on techn iques and p reconditi on ing techn iques fo r large scale m atrix p rob lem s are discu ssed . Som e p rob lem s that need to be fu rther studied are p resen ted .
Key words :m atrices ; linear system s ; eigenvalue ; K rylov sub sp ace m ethods