krylov子空间算法
Krylov子空间方法II
max {q (λi )2 }
2 yi λi
min
max {q (λi ) } y Λy
5/76
= = 即
q ∈Pk , q (0)=1 1≤i≤n q ∈Pk , q (0)=1 1≤i≤n
min
max {q (λi )2 } ϵ0 Aϵ0 max {q (λi )2 } ∥ϵ0 ∥2 A,
⊺
min
11/76
(4.6)
又 Λ 是对角矩阵, 所以 ∥q (Λ)∥2 = max |q (λi )|.
1≤i≤n
设 x(k) 是由 GMRES 方法得到的近似解. 由 GMRES 方法的最优性可 知, x(k) 极小化残量的 2 范数. 因此, ∥b − Ax(k) ∥2 = = ≤
x∈x(0) +Kk (A,r0 ) q ∈Pk , q (0)=1 q ∈Pk , q (0)=1
1+ε δ
10/76
5.2 GMRES 方法的收敛性
正规矩阵情形 设 A 是正规矩阵, 即 A = U ΛU ∗ , 其中 Λ = diag(λ1 , λ2 , . . . , λn ) 的对角线元素 λi ∈ C 为 A 的特征值. 设 x ∈ x(0) + Kk (A, r0 ), 则存在多项式 p(t) ∈ Pk−1 使得 x = x(0) + p(A)r0 . 于是 b − Ax = b − Ax(0) − Ap(A)r0 = (I − Ap(A))r0 ≜ q (A)r0 , 其中 q (t) = 1 − t p(t) ∈ Pk 满足 q (0) = 1. 直接计算可知 ∥b − Ax∥2 = ∥q (A)r0 ∥2 = ∥U q (Λ)U ∗ r0 ∥2 ≤ ∥U ∥2 ∥U ∗ ∥2 ∥q (Λ)∥2 ∥r0 ∥2 = ∥q (Λ)∥2 ∥r0 ∥2 .
keryolv子空间迭代法
keryolv子空间迭代法Krylov子空间迭代法是一种求解大规模线性方程组的有效方法。
它的基本思想是利用一个初始向量和一个矩阵来构造一个Krylov子空间，然后在这个子空间中寻找一个近似解。
这种方法通常比直接求解线性方程组的方法更快，尤其是当矩阵非常大时。
下面将从以下几个方面详细介绍Krylov子空间迭代法：1. Krylov子空间的定义和构造Krylov子空间是由一个向量v和一个矩阵A产生的一组向量集合，表示为：K(A,v) = span{v, Av, A^2v, ..., A^(k-1)v}其中k是任意正整数。
这个集合包含了所有由v和A作用k次得到的向量的线性组合。
2. Arnoldi过程Arnoldi过程是一种构造Krylov子空间的算法。
它通过对向量集合进行正交化来构造一个Hessenberg矩阵，该矩阵描述了向量在Krylov 子空间中的投影。
Arnoldi过程可以表示为以下步骤：(1) 选择初始向量v，并令q1 = v/||v||。
(2) 对于k = 1, 2, ..., n，执行以下步骤：(a) 计算w = Aqk。
(b) 对于j = 1, 2, ..., k，计算hj,k = qj^Tw，并令w = w - hj,kqj。
(c) 计算hk+1,k = ||w||，如果hk+1,k=0，则停止迭代。
(d) 如果hk+1,k≠0，则令qk+1 = w/hk+1,k，并将(h1,1, h2,1, ..., hk+1,k)作为Hessenberg矩阵的第k列。
3. GMRES方法GMRES是一种基于Krylov子空间的迭代方法，用于求解线性方程组Ax=b。
它通过在Krylov子空间中寻找一个最小化残差的向量来逼近解向量。
GMRES可以表示为以下步骤：(1) 选择初始向量x0和r0=b-Ax0。
(2) 构造Krylov子空间K(A,r0)，并使用Arnoldi过程构造Hessenberg矩阵H和正交矩阵Q。
Krylov子空间方法
由于 x ˜ ∈ x(0) + K, 因此存在向量 y ∈ Rm 使得 x ˜ = x(0) + V y 由正交性条件 (4.4) 可知 r0 − AV y ⊥ wi , i = 1, 2, . . . , m , 即 W ⊺ AV y = W ⊺ r0 .
x ˆ≜x ˜ − x(0) = V y
9/115
Arnoldi 过程: 计算 Km 的一组正交基
算法 2.1 基于 Gram-Schmidt 正交化的 Arnoldi 过程
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13:
给定非零向量 r, 计算 v1 = r/∥r∥2 for j = 1, 2, . . . , m − 1 do wj = Avj for i = 1, 2, . . . , j do hij = (wj , vi ) end for j ∑ wj = wj − hij vi hj +1,j = ∥wj ∥2 if hj +1,j = 0 then break end if vj +1 = wj /hj +1,j end for
若给定初值 x(0) ∈ Rn , 则改用仿射空间 x(0) + K, 即 find x ˜ ∈ x(0) + K such that b − Ax ˜ ⊥ L. (4.3)
好的初值一般都包含有价值 的信息
事实上, 如果将 x ˜ 写成: x ˜ = x(0) + x ˆ, 其中 x ˆ ∈ K, 则 (4.3) 就等价于 find x ˆ∈K such that r0 − Ax ˆ ⊥ L, (4.4)
定解条件
r = b − Ax ˜⊥L 其中 x ˜ 是近似解, L 是另一个 m 维子空间. 不同的 L 对应不同的投影方法 当 L = K 时, 我们称为 正交投影法 , 否则称为 斜投影法
Krylov子空间迭代法
采用IOM后，仍然需要存储v(1), v(2), …v(m),因为在第(vi)步 中仍然需要这些向量. 解决这个问题可以考虑采用H的LU分解，通过自身分解的迭代更新以减少每 一步的存储量 使xm的更新依赖于xm-1,
14
Arnoldi方法-DIOM
lower bidiagonal
banded upper triangular
15
Arnoldi方法-DIOM
16
Arnoldi方法-DIOM
17
Thanks for your time !
18
得到基于Galerkin原 理构成的算法
5
Arnoldi方法-基本算法
6
Arnoldi方法-基本算法
7
Arnoldi方法-MGS
8
Arnoldi方法-HO
9
Arnoldi方法-FOM
10
Arnoldi方法-FOM
11
Arnoldi方法-FOM（m）
12
Arnoldi方法-IOM
13
Arnoldi方法-DIOM
Krylov子空间方法
March 23, 2016
内
• Arnoldi算法
– Arnoldi过程 – Gram-Schmidt Arnoldi – HouseHolder Arnoldi
容
• 子空间和Krylov子空间
• FOM
– IOM – DIOM
2
子空间
• 空间
– 集合，元素都是向量 – 线性空间（向量空间）
• 线性空间（交换律，结合律，幺元性，零元性，可 逆性，数乘分配律等）
• 子空间
– 线性空间的非空子集
Krylov子空间、优化问题与共轭梯度法
Krylov 子空间、优化问题与共轭梯度法自动化 富晓鹏工程实践中经常需要求解大型线性系统KU=F 。
在很多情况下矩阵K 是非常稀疏的，比如来自偏微分方程的离散化等，此时矩阵中每行仅有较少的非零元素。
面临这样的问题，我们首先面对的问题是，应该采用直接消元法还是迭代方法。
对前者来说，为充分利用系数特性，节点重编号是重要的；而对后者来说，适当的预处理是关键。
本文将重点放在后一类方法中的一种进行介绍与分析，即共轭梯度法。
共轭梯度法适用于矩阵K 为对称阵的情况，算法本身简洁高效，且与一些其他的数学理论、概念相紧密联系，本文分析了共轭梯度法与Krylov 子空间，以及优化问题之间隐含的联系，并简要给出算法框架。
1. 线性方程组迭代解法与Krylov 子空间我们考虑迭代法求解线性方程组Ax=b 。
假定未采用预处理矩阵P ，或P 矩阵已经隐含在A 与b 中。
迭代法求解格式如下：1()k k P x P A x b +⋅=-⋅+ (1)为说明问题，我们考虑简单的迭代格式P=I ，并且x 1=b 。
则迭代的最初几步为：2()2x I A b b b Ab =-+=- (2)232()33x I A x b b Ab A b =-+=-+ (3) …由上面几个式子可得，以上迭代格式第j 步的解x j 是b ，Ab ，…，A j -1b 的线性组合。
当A 矩阵稀疏时，这些向量可以采用矩阵向量乘法的稀疏技巧很快得到。
以上发现自然与Krylov 子空间的概念相联系起来。
Krylov 矩阵： K j = [b Ab A 2b … A j -1b]Krylov 子空间：K j = b ，Ab ，…，A j -1b 的所有线性组合Krylov 命名了向量b ，Ab ，…，A j -1b 的全部线性组合构成的子空间，并认为在这一子空间中，有比上例中特定元素更与线性方程组的解相接近的元素。
共轭梯度法就是在这一子空间中，每一步迭代都依照某种标准寻求最优元素的线性方程组解法。
krylov子空间迭代法
krylov子空间迭代法Krylov子空间迭代法是一种有效的求解线性方程组的迭代方法，因Krylov于1908年提出而得名。
它是一种基于子空间的迭代方法，可以在较少的计算量下，解决高维线性方程组的较大特征值的问题。
Krylov子空间迭代法的基本思想是：将线性方程组中的高维系数矩阵P划分为n个受限的Krylov子空间，用这些子空间来模拟矩阵P的特征值的变化趋势。
这样，可以使线性方程组的解从低维子空间转移到高维子空间，从而求出线性方程组的解。
Krylov子空间迭代法具有以下优点：（1）采用Krylov子空间技术可以降低计算维度，减少计算量，提高计算效率；（2）将子空间技术与迭代法相结合，实现了近似求解线性方程组的解；（3）Krylov子空间迭代法能有效收敛，解的可靠性高；（4）运行简便，无需调整参数；（5）可用于求解各种类型的线性方程组。
由于Krylov子空间迭代法的优越性，它已经广泛应用于工程、数学、物理、生物等多学科的计算和仿真中。
从根本上讲，Krylov子空间迭代法是一种非常有效的迭代方法，它可以有效地解决线性方程组的特征值问题。
下面我们将介绍Krylov 子空间迭代法的算法步骤：（1）输入高维系数矩阵P、初始向量v、迭代次数m及收敛准则ε；（2）构造Krylov子空间：V＝[v，Pv, Pv,……，P^m-1v]；（3）用V中的向量代替P，将Pv-λv转化为V的线性方程；（4）求解V线性方程组；（5）求出V的特征值λ；（6）利用第4步求出的解v，求出线性方程组的解x；（7）若特征值收敛，则停止迭代；（8）重复第2至第7步，直至特征值收敛；（9）输出计算结果。
以上就是Krylov子空间迭代法的算法步骤。
Krylov子空间迭代法的算法实现起来相对简单，只需要实现以上的几个步骤即可。
由于Krylov子空间迭代法的有效性，它已经被广泛应用于工程、数学、医学、物理、生物等多学科的计算和仿真中。
总之，Krylov子空间迭代法是一种高效的求解线性方程组的迭代方法，它可以有效收敛，具有较高的求解精确度和计算效率。
krylov子空间算法
Krylov 子空间的定义：定义：令N R υ∈，由1m A υυυ-，，，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 为采用投影方法得到的新近似解的充要条件是()()()()01min z x Kx z ϕϕ∈+=其中,()()()12,z A x z x z ϕ=--2.(残量投影)设A 为任意方阵,()0x 为初始近似解,且L AK =,则()1x 为采用投影方法得到的新近似解的充要条件是()()()()01min z x Kx z ψψ∈+=其中()()122,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 、DGMERSArnoldi 方法标准正交基方法：Arnoldi 方法是求解非对称矩阵的一种正交投影方法。
(完整版)Krylov子空间迭代法
February 10, 2020
内容
• 子空间和Krylov子空间
• Arnoldi算法
– Arnoldi过程 – Gram-Schmidt Arnoldi – HouseHolder Arnoldi
• FOM
– IOM – DIOM
2
子空间
• 空间
– 集合，元素都是向量 – 线性空间（向量空间）
根据Cayley-Hamilton定理有
������������ + ������������−1������������−1+. . . +������1������1 + ������0������0 = 0
即
VP= -������������ 其中������ = [������0, ������1, . . . , ������������−1ሿ，������ = ������0, ������1, . . . , ������������−1 ������ Krylov子空间： ������������(������, ������) = ������������������������{������, ������������, . . . , ������������−1������ሽ Krylov矩阵： ������������(������, ������) = [������, ������������, . . . , ������������−1������ሿ
• 线性空间（交换律，结合律，幺元性，零元性，可 逆性，数乘分配律等）
• 子空间
– 线性空间的非空子集
• 包含零元素，并且满足加法和乘法的封闭性
– 扩张（符合记作span）
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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 为采用投影方法得到的新近似解的充要条件就是()()()()01min z x Kx z ϕϕ∈+=其中,()()()12,z A x z x z ϕ=--2.(残量投影)设A 为任意方阵,()0x 为初始近似解,且L AK =,则()1x 为采用投影方法得到的新近似解的充要条件就是()()()()01min z x Kx z ψψ∈+=其中()()122,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 、DGMERSArnoldi 方法标准正交基方法:Arnoldi 方法就是求解非对称矩阵的一种正交投影方法。
Arnoldi 算法就就是对非对称矩阵A,产生Krylov 子空间()()0,m A r K 的一组标准正交基的方法。
该算法构造()()0,m A r K 的一组标准正交基12,,,m v v v K 与Hessenberg 矩阵1112131212223232331,,100m m m m m m m mm h h h h h h h h h h H h h h --⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦L LM L M O O O L,1112131212223232331,,11,0000m m m m m mm m m m m h h h h hh h h h h H h h h h --+⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦L LML M O O O L L基于Gram-Schmit 正交化方法首先,选取一个Euclid 范数为1的向量()1v ,对(),m A κν,通常可取()()112,νννν=,在()()()12,,,j v v v L 已知的情况下,不妨设()()()()12,,,,jj v v v Av L 线性无关(否则构造完毕),则可求出与每个都正交的向量()()()()()1212jjjj j jj w Av h v h v h v =----L而不难瞧出()()(),,1,2,,j i ij h Av vi j==L ,再记()()()121,,jj j j h w w +=,得到与()()()12,,,j v v vL 都正交的向量()()11,jj j jw vh ++=,重复此过程,即可得到一组标准正交基。
若期间某个j 使得1,0j j h +=,则说明v 的次数就是j,且j κ就是A 的不变子空间。
Arnoldi 算法: （1）取向量()1v ,满足()11,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 jj j jj w Av h v h v h v =----L (1)()()(),,1,2,,ji ij h Av v i j ==L(2)()()()121,,j j j j h ww+=(3) ()()11,jj j jw vh ++=(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 ,则:()()()1Tm m m m m m m AV V H w e V H +=+=Tm m m V AV H =在Arnoldi 算法中,可能有较大舍入误差,改写:()()()()()1111,,,1,2,,ji iij j i j h Av h v h v v i j --=---=L L修正的Arnoldi 算法: （1）取向量()1v ,满足()11,1v j ==（2）计算()()j j w Av =（3）依次对1,2,,i j =L ,计算()()(),j i ij h w v =与()()()jjiij 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,,,mV 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 非奇异时,有()()10T T y W AV W r -=,从而得到迭代公式:()()()()110T T x x V W AV W r -=+FOM 算法: (1)计算()()00rb Ax=-,()()()1200,rrθ=,()()01rv θ=,置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()()()()()()011,m mm 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）计算()()()()()()()10200001,,,rrb Ax r rv θθ=-==（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）计算()()00rb Ax=-,()()()1200,rrθ=,()()()010rv r =,置1m =（2）计算()()m m w Av = （3）对()max 1,1,,i m k m =-+L ,依次计算()()(),m i im h w v =与()()()mmi ijw w h v =-（4）计算()()()()()1211,,,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)()()()11m m m i mmimpuv up--⎛⎫ ⎪=- ⎪⎝⎭∑ (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 方法就是求解大规模稀疏对称矩阵端部特征问题的一种常用的正交投影方法。