Krylov 子空间的定义：
定义：令N R υ∈，由1m A υυυ-L ，，，A 所生成的子空间称之为由υ与A 所生成的m 维Krylov 子空间，并记(),m K A v 。
主要思想是为各迭代步递归地造残差向量，即第n 步的残差向量 ( )
n r 通过系数矩阵A 的某个多项式与第一个残差向量()0r 相乘得到。即
()()( )
0n r p A r =。
但要注意，迭代多项式的选取应该使所构造的残差向量在某种内积意义下相互正交，从而保证某种极小性（极小残差性），达到快速收敛的目的。
Krylov 子空间方法具有两个特征：1.极小残差性，以保证收敛速度快。2.每一迭代的计算量与存储量较少，以保证计算的高效性。
投影方法
线性方程组的投影方法
方程组Ax b =,A 是n n ?的矩阵。给定初始()0x ,在m 维空间K(右子空
间）中寻找x 的近似解()1x 满足残向量()1r b Ax =-与m 维空间L(左子空间）正交，即()1b Ax L -⊥，此条件称为Petrov-Galerkin 条件。 当空间K=L 时，称相应的投影法为正交投影法，否则称为斜交投影法.
投影方法的最优性：
１. (误差投影)设A 为对称正定矩阵,()0x 为初始近似解,且K=L,则()1x 为采用投影方法得到的新近似解的充要条件是()()()()0
1min z x K
x z ??∈+= 其中,()()()1 2,z A x z x z ?=--
2.(残量投影)设A 为任意方阵,()0x 为初始近似解,且L AK =,则()1x 为采用投影方法得到的新近似解的充要条件是()()()()0
1min z x K x z ψψ∈+= 其中()() 12
2,z b Az b Az b Az ψ= -=--
矩阵特征值的投影方法
对于特征值问题Ax x λ=,其中A 是n ×n 的矩阵，斜交投影法是在m 维右子空间K 中寻找i x 和复数i λ满足i i i Ax x L λ-⊥,其中L 为m 维左子空间.当L=K 时，称此投影方法为正交投影法. 误差投影型方法： 取L=K 的正交投影法
非对称矩阵的FOM 方法（完全正交法） 对称矩阵的IOM 方法和DIOM 方法 对称矩阵的Lanczos 方法 对称正定矩阵的CG 方法 残量投影型方法： 取L=AK 时的斜交投影法
GMERS 方法（广义最小残量法）
重启型GMERS 方法、QGMERS 、DGMERS Arnoldi 方法 标准正交基方法：
Arnoldi 方法是求解非对称矩阵的一种正交投影方法。Arnoldi 算法就是对非对称矩阵A ，产生Krylov 子空间()()0,m A r K 的一组标准正交
基的方法。该算法构造()()0,m A r K 的一组标准正交基12,,,m v v v K 和Hessenberg 矩阵
111213121 2223
232331,,1 00
m m m m m m m mm h h h h h h h h h h H h h h --?? = L L
M L M O O O L ，111213121 2223 232331,,1 1,00 00
m m m m m mm m m m m h h h h h h h h h h H h h h h --+=?????? L L M
L M O O O L L
基于Gram-Schmit 正交化方法
首先，选取一个Euclid 范数为1的向量()1v ,对(),m A κν,通常可取
()() 11 2 ,ν ννν= ，在()()( )
12,,,j v v v L 已知的情况下,不妨设()()()( ) 12,,,,j
j v v v Av L 线性
无关(否则构造完毕),则可求出与每个都正交的向量 ()()()()()1 2 12j j j
j j jj w Av h v h v h v =----L 而不难看出 ()() (
),,1,2,,j i ij h Av v i j
==L ，再记()( ) ( ) 1 2 1,,j
j j j h w w +=，得到与 () () () 12,,,j v v v L 都正交的向量() ()11,j
j j j w v h ++=
，重复此过程，即可得到一组标准
正交基。若期间某个j 使得1,0j j h +=，则说明v 的次数是j ，且j κ是A 的不变子空间。 Arnoldi 算法： （1）取向量()1v ，满足
()1 1,1v j ==
（2）按（2）式计算,1,2,,ij h i j =L ，再按（1）式计算()j w （3）按（3）式计算1,j j h +，若1,0j j h +=，则停止，否则按（4）式计算()1j v + （4）若j m <，则1j j =+，转（2），否则停止
()()()()()1212j j j
j j jj w Av h v h v h v =----L （1） ()( ) ( ),,1,2,,j
i ij h Av v i j ==L （2） () () ( ) 1 2
1,,j j j j h w w += （3） ( )
()11,j j j j w v h ++= （4）
定理：如果记以()()( )
12,,,m v v v L
为列构成的矩阵为m V ,由ij h 定义的(m+1)×m
阶上Hessenberg 矩阵（假设一个n n ?阶矩阵A ，在1i j > +时，
它的0ij a =，那么这个矩阵A 就叫做Hessenberg 矩阵）为m H ,删除最后一行得到的矩阵为m H ，则：()
() () 1T
m m m m m m m AV V H w e V H +=+= T
m m m V AV H =
在Arnoldi 算法中，可能有较大舍入误差，改写： ()()( ) ()( ) 1
111,,,1,2,,j i i
ij j i j h Av h v h v v i j --=---=L L
修正的Arnoldi 算法： （1）取向量()1v ，满足()1 1,1v j ==
（2）计算()()j j w Av =
（3）依次对1,2,,i j =L ，计算()()(),j i ij h w v =与 ()()()j j i
ij w w h v =-
（4）计算1,j j h +，若1,0j j h +=，则停止，否则计算()1j v + （5）若j m <，则1j j =+，转（2），否则停止
FOM （完全正交化）方法
非对称矩阵的FOM 方法： 解方程组的投影法的矩阵表示 设n m ?阶矩阵()() ()12,,,m V v v v ??=??L
与()()()12,,,m W w w w ??=?? L 的列分别构成K
与L 的一组基。记()()10,z x x z Vy =-=，有()()00T W r AVy -= 当T W AV 非奇异时，有()
( ) 1
0T T y W AV W r -=，从而得到迭代公式： ()()()( ) 1 1
0T T x x V W AV W r -=+ FOM 算法： （1）计算() () 00r b Ax
=-，() () () 1 2 00,r r θ=，()() 01r v θ =，置1j =
（2）计算()()j j w Av =
（3）依次对1,2,,i j =L ，计算()()(),j i ij h w v =与()()()j j i ij w w h v =- （4）计算1,j j h +，若1,0j j h +=，则置m j =，转（6） （5）计算()1j v +，若j m <，则置1j j =+，转（2） （6）按下式计算()m x
( ) ()()( ) () ()0 1 1,m m
m m m x x V y y H e θ-=+=
不难看出，当采用上述FOM 算法时，需要存储所有的i v ，(i=1,2,…m)，当m 增大时，存储量以()O mn 量级增大.而FOM 计算量是()2O m n .可见其代价十分高昂.因此我们考虑重启的 FOM 算法 重启型FOM 算法： （1）计算()
()
()() ( ) ()() 1 02 00001 ,,,r r b Ax r r v θθ =-==
（2）生成()()0,m A r κ的一组标准正交基，得到m H
（3）按下式计算()m x ，若()m x 满足精度要求，则停止，否则置()()0m x x =，
转（1）。()()()()()()011
,m m m m m x x V y y H e θ-=+= IOM 方法
非对称矩阵的IOM 方法
所谓不完全正交化方法(IOM),是指在正交化过程中,()1j v +仅与最近k 个()(
) ( )
11,,,j k j j v v v -+-L
正交,这样做虽然破坏的正交性,但是降低了计算量.
当然k 选得越小，对每个j 对应的计算量也越小,但可能要选更大的m 才能取得满足精度要求的近似解.
IOM 算法仅仅是把FOM 算法的第三步改为()max 1,1,,i j k j =-+L
，
计算ij h 与()j w 。
但采用IOM 后，仍然需要存储()()( )
12,,,m v v v L ,因为在第(vi)步 ( ) ()( )
m m m x x V y =+中仍然需要这些向量.
解决这个问题可以考虑采用H 的LU 分解，通过自身分解的迭代 更新以减少每一步的存储量 DIOM 算法： （1）计算() () 00r b Ax =-，() () () 1 2 00,r r θ=，()() ( ) 010r
v r =，置1m =
（2）计算()()m m w Av = （3）对()max 1,1,,i m k m =-+L ，依次计算()()
(),m i im h w v =与()()()m m i ij w w h v =- （4）计算()()() ( ) ( ) 12
11,,,m m m m m m h w w v w hm ++==
（5）按（4）式更新m H 的LU 分解，若0mm u =，则停止 （6）按（3）式计算m ξ，按（2）式计算()m p ，其中0i ≤时，()0i im u p = （7）按（1）式计算()m x ，若符合精度要求，则停止，否则1m m =+，转（2）
( ) ( ) ( )
1m m m m x x p ξ-=+ （1） () ()()1 1
m m m i mm im p u v u p
--? =- ∑ （2） 1,110 ,1{
m m m m l m ξξξ--==-> （3） 1,1,m k m m k m u h -+-+=
,11,,1,,1im im i i i m u h l u i m k m --=-=-+-L （4） ,1,1,11,1,1
,m m m m mm mm m m m m m m h l u h l u u ------= =-
Lanczos 方法
Lanczos 方法是求解大规模稀疏对称矩阵端部特征问题的一种常用的正交投影方法。它在Krylov 子空间上对对称矩阵采用Rayleign-Ritz 方法，得到某些特征值和特征向量的近似。基本思想是给定一初始向量，逐步地构造出由该初始向量和对称阵生成的Krylov 子空间的标准正交基。通过简单的三项递推公式将大规模对称阵投影为小规模对称三对角阵，然后用此三对角阵的特征对来得出原始矩阵的近似特征对。由于三对角阵好的结构和小的维数，在准确运算下，每一步所需存储量和计算量是常量，不会随着子空间维数的增加而改变。
Lanczos 算法： （1）计算() () 00r b Ax =-，() () ()
1 2 00,r r θ=，()() 01r v θ =，置1j =
（2）计算()()()1j j j j w Av v θ-=-，其中当1j =时()10j j v θ-= （3）计算()()(),j j j w v δ=与()()()j j j j w w v δ=- （4）计算()
() () 1 2 1,j j j w w
θ+=，若10j θ+=，则置m j =转（5），否则置 ( ) () 11 j j j w
v θ++=，若j m <，则1j j =+，转（2）
（5）置()()()111,,,,,m m i i i m T tridiag V v v θδθ++??==??L ，计算
( ) () ()1 1
m m y T e θ-=
（6）计算()()()0m m m x x V y =+ Lanczos 双正交化方法 在双正交化过程中，取 ( ) ( )()()()00 01
,{,,,}m m K A r span r Ar A r -=L ( ) ( ) ()()() ()1 00
,{,,,}m T T T m L A r span r A r A r -=L
容易看出A 和T A 在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是A 和T A 的方程组.
基于Lanczos 双正交化过程的双共轭梯度法BiCG 算法： （1）计算()()00r b Ax =-，选取()0t r 使得()()()00,0t r r ≠，置方向向量()()00s r =，
()( )
0t t s r =，并置m=0 （2）计算 ()() ()()()) ,,m m m
m m t t
r r As s α=与向量( ) ( ) ( )
1m m m m x x s α+=+， ( ) ( ) ()
1m m m m r r As α+=-，()()()m m m T t t m r r A t α=-，如果()1m x +满足精度要求，则停止
（3）计算参数 ()()
()()()()11,,m m m m m t t
r r r r β++=，与向量( ) ( ) ( )
11m m m m s r s β++=+， ( ) ( )
( )
11m m m t t m t s r s β++=-，置m=m+1,转（2） CG 方法
CG 法用来解对称且正定方程组十分有效，但若是拿来解非对称或是非正定的线性方程组则会发生中断。它是借由迭代的方式产生一序列的方向向量用来更新迭代解以及残向量，虽然产生的序列会越来越大，但是却只需要存储少数的向量。当系数矩阵A 相当大而且稀疏，此时CG 法几乎就是高斯消去法。CG 法理论上虽然保证最多n 步能解出线性方程组Ax b =的解，但是若系数矩阵是病态的，此时累进误差会让CG 法在n 步后无法求得充分精确的近似解。 CG 算法：
（1）计算()()00r b Ax =-，选取()()00s r =，置0j = （2）计算参数()()()
()()),,j j j j
j r r As s α=，更新向量( ) ()( ) 1j j
j j x x s α+=+与残向量 ( ) ()( ) 1j j
j j r r As α+=-，若()j x 满足精度要求，则停止 （3）计算()()()()())( ) (
) ( )
111111,,,j j j j
j j j j j r r r r s r s ββ++++++==+，置1j j =+，转（2） CG 法是解正定对称线性方程组最有效的方法之一，该方法充分利用了矩阵A 的稀疏性，每次迭代的主要计算量是向量乘法。
GMRES 方法
GMRES 算法要求系数矩阵A 是非奇异，非对称的n 维方阵。GMRES 算法利用Arnoldi 变换构造一正交基[]12,,,m m V v v v =L 来表示Krylov 子空
间m κ。
GMRES 方法把方程组的求解化为一个极小问题的求解，不去直接求
()1
x ,转而去解此极小问题，这是GMRES 方法的独到之处。 ()()() ()( )()( ) ( ) ()( )01 1 1 122 2 2
2
m m m m m m m m m z b Az r AV y v AV y V e H y e H y ψθθθ+=-=-=-=-=- GMRES 算法的收敛性完全取决于系数矩阵A 的特征值的性质。 GMRES 算法：
（1）计算() () 00r b Ax =-，() () () 1 2 00,r r θ=，()() 01r v θ =，置1j =
（2）计算()()j j w Av =
（3）依次对1,2,,i j =L ，计算()()(),j i ij h w v =与()()()j j i ij w w h v =- （4）计算1,j j h +，若1,0j j h +=，则置m j =，转（6） （5）计算()1j v +，若j m <，则置1j j =+，转（2） （6）求()()12arg min m m z
y e H z θ=-，计算()()()0m m m x x V y =+ 求解最小二乘问题()()12arg min m m z y e H z θ=-的算法： （1）令1g θ=，1i = （2）计算()1
222 1.1,,,ii i i ii i i h h h h c s ρρρ++=+==，置ii h ρ=
（3）置向量.1:i i m t h +=，计算,1:1,1:1,1:1,1:,i i m i i m i i m i i m h ct sh h ch st +++++++=+=-（.1:i i m h +表示矩阵第i 行从i+1列到m 列的元素构成的列向量） （4）置0i t g =，计算010,i i g ct g st +==- （5）若i m <，转（2） （6）依次对
,1,,1j m m =-L ，计算() ,11j jm m j j j j jj g
h g h g g h ++---= L 所得到的
[]12,,,m g g g L 即为所求的()m y
GMRES 算法允许Krylov 子空间维数增加到n ，而且总是在最大迭代次数n 内终止运算；另一种重启型GMRES 算法严格要求子空间维数为一个定值m ，在进行了m 次迭代后，以得到的最后迭代结果m x 作为初始点重新进行Arnoldi 变换，当残余向量r A bx =-满足0T r Ar =时，终止计算。综合考虑时间和空间复杂度，重启型的GMRES 算法更适合。 重启型GMRES 算法： （1）计算()()()(
) ( ) ()() 1 02 01 ,,,r
r b Ax r r v θθ =-==
（2）生成()()0,m A r κ的一组标准正交基，得到m H （3）求()()12arg min m m z
y e H z θ=-，计算()()()0m m m x x V y =+，若()m x 满足精度
要
求，则停止，否则置()()0m x x =，转（1）
同样可以采取不完全正交的方法，在正交化过程中,()1j v +仅与最近k 个()(
) ( )
11,,,j k j j v v v -+-L
正交就能得到QGMRES 方法。为了避免存储 ()()( ) 1 2
,,,m k v v v -L ，可以对m H 进行QR 分解.这种算法被称为DQGMRES 。
QGMRES 算法： （1）计算() () ()() ( ) ()() 1 02 00001 ,,,r r b Ax r r v θθ =-==
（2）计算()()j j w Av = （3）对()max 1,1,,i j k j =-+L ，计算()()
(),j i ij h w v =与()()()j j i ij w w h v =-
（4）计算1,j j h +，若1,0j j h +=，则置m j =，转（6） （5）计算()1j v +，若j m <，则置1j j =+，转（2）
（6）求()()12 arg min m m z y e H z
θ=-，计算()()()0m m m x x V y =+ Krylov 子空间方法解矩阵特征值问题
Arnoldi 方法求解非对称矩阵特征值 由定理：()()()1T m m m m m m m AV V H w e V H +=+= T
m m m V AV H =
而()()11,m m m m w h v ++=，则有如下结论：
（1）若1,0m m h +=，则(),m A v κ是A 的不变子空间，从而m H 的所有特征值是A 的特征值子集。
（2）若1,0m m h +≠，则m H 的特征值对是(),i i y λ，即m i i i H y y λ=，由上述定理可得：()
() () () () 11,1,T T
m m m m i i m i m m i m m i AV y V y h v e y h e
y λ+++-==
以上讨论说明,可以用m H 的特征值 (又称Ritz 值)来近似A 的特征值,用向量i m i x V y =(又称Ritz 向量)来近似相应的特征向量.事实上, m H 为A
在Krylov 子空间上的正交投影.
由于m H 是上Hessenberg 阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解. 求解矩阵特征值的Arnoldi 方法： （1）给定m ，取向量()1v ，满足
()1 1,1v j ==
（2）计算()()j j w Av =
（3）依次对1,2,,i j =L ，计算()()(),j i ij h w v =与 ()()()j j i
ij w w h v =-
（4）计算1,j j h +，若1,0j j h +=，则停止，否则计算()1j v + （5）若j m <，则1j j =+，转（2）
（6）计算()m ij H h =的特征值对(),i i y λ及A 关于i λ的Ritz 向量i m i x V y = Arnoldi 算法构造标准正交基存在的问题： 1需要存储所有的基向量,当m 很大时,存储量大 2理论上为了保证收敛速度,m 越大越好 Lanczos 方法求解对称矩阵特征值
A 是对称阵时,m H 是三对角阵，仍然采用m H 的特征值来近似A 的特征值，相应的Ritz 向量来近似A 的特征向量.
问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，QR 法等，可以在并行机上达到tO(logN)的量级.
求解对称矩阵特征值的Lanczos 方法： （1）计算() () 00r
b Ax =-，() () () 1 2 00,r r θ=，()() 01r v θ =，置1j =
（2）计算()()()1j j j j w Av v θ-=-，其中当1j =时()10j j v θ-= （3）计算()()(),j j j w v δ=与()()()j j j j w w v δ=-
（4）计算()()()1
21,j j j w w θ+=，若10j θ+=，则置m j =转（5），否则置 () ( ) 11 j j j w v
θ++=，若j m <，则1j j =+，转（2）
（5）计算()m ij H h =的特征对(),i i y λ及A 的关于i λ的Ritz 向量i m i x V y = 预条件法
Krylov 子空间方法的收敛性一般都和矩阵的特征值分布有关，特征值分布越集中，收敛速度越快，若特征值分布较分散，则收敛的很慢，此时可以采用预条件子来改善矩阵的特征值分布.
选取矩阵T M LL =使得A 的特征值比A 集中，并且1M -容易求得，于是将原问题化为问题11M Ax M b --=,这是左预条件法，还可
采用右预条件法和分裂预条件法
在求解矩阵特征值问题时，可以利用chebyshev 多项式来加快收敛
所谓多项式加速,就是采用()
()()11v q A v =%作为迭代初始向量,其中()q ? 为多项式,如果()11 n i i i v p ξ ==
∑, 其中i p 为A 的特征值对应的特征向量 则() (1) 1 1 1 ()()()()n r n
n i i i i i i i i i i i i r v q A v q p q p q p ξλξλξλ===+===+∑∑∑。 将特征值按实部递减顺序排列,其中为r 个“想要”的特征值,将“不想要”的特征值点集用S 表示.如果多项式()q ?,能使得S 尽可能小,则相当于让求解过程产生加速.而Chebyshev 多项式由于它的特殊性质,恰好能满足这种要求.不过，至今,仍无有效方法确定Chebyshev 多项式的某些参数.
块方法、调和投影法
当矩阵的特征值分布较密集或有重特征值时，常规的Arnoldi 方法或Lanczos 方法就必须取m 很大时才能收敛，此时可以采用块方法 取1(,){,,....}m
m A V span V AV A V ξ-=可以使得无须很大的m 就可让迭代收 敛，并且此方法适合改为并行算法。
当特征值密集时，还可采用调和投影法，这也是最近研究的一个 热点。当左右子空间的基满足m m W AV 时，称为调和投影，它除了可以通过选取位移点改善原矩阵的特征值分布外，还可以将内部特征值问题化为外部特征值问题.
因篇幅问题不能全部显示，请点此查看更多更全内容