6.4 加速方法
当迭代解 x(0),x(1),x(2),…,x(k) 已经计算出来后,我们可以对其进行组合,得到一个新的近似解,这样就可以对原方法进行加速。
6.4.1 外推技术
设原迭代格式为
x(k+1)=Gx(k)+g.(6.23) 由 x(k) 和 x(k+1) 加权组合后可得新的近似解
x(k+1)=(1−ω)x(k)+ω(Gx(k)+g),(6.24) 其中 ω 是参数. 这种加速方法就称为外推方法
如果原迭代格式是 Gauss-Seidel 迭代, 虽然 SOR 与外推方法 (6.17) 都是通过加权进行加速,但格式是不一样的: 前者是在单个分量计算出来后就进行加权, 而后者则是当所有分量都计算出来后才加权.
为了使得迭代格式(6.17)尽可能快地收敛, 需要选择 ω 使得其迭代矩阵 Gω≜(1−ω)I+ωG 的谱半径尽可能地小. 假设 G 的特征值都是实数, 且最大特征值和最小特征值分别为 λ1 和 λn . 于是
ρ(Gω)=λ∈σ(G)max∣(1−ω)+ωλ∣=max{∣1−ω+ωλ1∣,∣1−ω+ωλn∣}. 定理6.28 设 G 的特征值都是实数, 其最大和最小特征值分别为 λ1 和 λn , 且 1∈/[λn,λ1] , 则
ω∗=argωminρ(Gω)=2−(λ1+λn)2, 此时
ρ(Gω∗)=1−∣ω∗∣d, 其中 d 是1到 [λn,λ1] 的距离,即当 λn≤λ1<1 时, d=1−λ1 ,当 λ1≥λn>1 时, d=λn−1
(板书)
证明. 易知 ∣1−ω+ωλ1∣=∣1−ω+ωλn∣ 当且仅当 λ1=λn 或 ω=0 或 ω=ω∗=2−(λ1+λn)2 .
若 λ1=λn ,则当 ω=1−λ11=ω∗ 时 ρ(Gω∗) 达到最小值0.
下面假定 λ1=λn . 先考虑 λ1<1 的情形. 此时有
λn≤λ≤λ1max∣(1−ω)+ωλ∣=⎩⎨⎧1−ω+λnω,1−ω+λ1ω,−(1−ω+λnω),ω≤0;0<ω≤ω∗;ω>ω∗. 通过函数图像可知, 当 ω=ω∗ 时取最小值. 此时
ρ(Gω∗)=1−(1−λ1)ω∗=1−2−(λ1+λn)2(1−λ1)=2−(λ1+λn)λ1−λn. 对于 λn>1 的情形,可以进行类似的讨论
由定理6.28可知, ρ(Gω∗)=1−∣ω∗∣d, 且当 ω∗=1 时,外推送代(6.17)比原迭代法收敛要更快一些.
最优参数依赖于原迭代矩阵 G 的特征值, 因此实用性不强. 在实际应用时可以估计特征值所在的区间 [α,β] , 然后用 α,β 来代替 λn 和 λ1 .
JOR方法
对Jacobi迭代进行外推加速, 则可得JOR(Jacobi over-relaxation)方法:
x(k+1)=(1−ω)x(k)+ω(D−1(L+U)x(k)+D−1b)=x(k)+ωD−1(b−Ax(k)),k=0,1,2,…. 定理6.29 设 A 对称正定. 若
0<ω<ρ(D−1A)2, 则JOR方法收敛
6.4.2 Chebyshev多项式加速
本节对外推技巧进行推广.
假定通过迭代格式 (6.16) 已经计算出 x(0),x(1),…,x(k) , 下面考虑如何将这些近似解进行组合, 以便得到更精确的近似解.
记 εk=x(k)−x∗ 为第 k 步迭代解的误差, 则有
εk=Gεk−1=G2εk−2=⋯=Gkε0. 设 x~(k) 为 x(0),x(1),…,x(k) 的一个线性组合,即
x~(k)=α0x(0)+α1x(1)+⋯+αkx(k),(6.25) 其中 αi 为待定系数,且满足 ∑i=0kαi=1. 于是
x~(k)−x∗=α0ε0+α1Gε0+⋯+αkGkε0≜pk(G)ε0,(6.26) 其中 pk(t)=∑i=0kαiti 为 k 次多项式, 且满足 pk(1)=1 .
我们希望通过适当选取参数 αi , 使得 x~(k)−x∗ 尽可能地小, 即使得 x~(k) 收敛到 x∗ 速度远远快于 x(k) 收敛到 x∗ 速度. 这种加速方法就称为多项式加速或半迭代法 (semi-iterative method).
例6.4 设 pn(t) 为 G 的特征多项式, 令 p~n(t)≜pn(t)/pn(1) , 则 p~n(1)=1 且 p~n(G)=0 , 所以选取 αi 为 p~n 的系数, 则 x~(n)−x∗=0 . 但这种选取方法不实用, 原因是:
(1) p~n(t) 的系数并不知道;
(2) 我们通常希望收敛所需的迭代步数 ≪n .
下面讨论参数 αi 的较实用的选取方法.由(6.19)可知
∥x~(k)−x∗∥2=∥pk(G)ε0∥2≤∥pk(G)∥2⋅∥ε0∥2. 因此我们需要求解下面的极小化问题
p∈Pk,p(1)=1min∥p(G)∥2,(6.27) 其中 Pk 表示所有次数不超过 k 的多项式组成的集合. 一般来说, 这个问题是非常困难的. 但在一些特殊情况下, 我们可以给出其 (近似) 最优解.
假设迭代矩阵 G 是对称矩阵,即 G 存在特征值分解
G=QΛQT, 其中 Λ 是对角矩阵, 且对角线元素都是实的, Q 是正交矩阵. 于是有
minp∈Pk,p(1)=1∥p(G)∥2=minp∈Pk,p(1)=1∥p(Λ)∥2=minp∈Pk,p(1)=1max1≤i≤n{∣p(λi)∣}≤minp∈Pk,p(1)=1maxλ∈[λn,λ1]{∣p(λ)∣},(6.28) 其中 λ1,λn 分别表示 G 的最大和最小特征值. 这是带归一化条件的多项式最佳一致逼近问题 (与零的偏差最小). 该问题的解与著名的 Chebyshev 多项式有关.
由于所有算子范数 ∥pk(G)∥ 的下确界是 ρ(pk(G)) ,因此,一种较实用的选取方法是使得 pk(G) 的谱半径尽可能地小.
Chebyshev多项式
Chebyshev多项式是一类很重要的正交多项式, 在函数逼近, 函数插值, 数值积分等方面都有着重要的应用.
Chebyshev多项式 Tk(t) 可以通过下面的递归方式来定义:
T0(t)=1,T1(t)=t, Tk(t)=2tTk−1(t)−Tk−2(t),k=2,3,…,(6.29) 也可以直接由下面的式子定义
Tk(t)={cos(karccost),cosh(karccosht),∣t∣≤1∣t∣>1, 其中 cosh 为双曲余弦, 即 cosh(t)=2et+e−t .
引理6.30 Chebyshev多项式具有下面的性质:
(1) Tk(1)=1
(2) Tk(t) 的首项系数为 2k−1 ;
(3) Tk(−t)=(−1)kTk(t) , 即 T2k(t) 只含偶次项, T2k+1(t) 只含奇次项;
(4) 当 ∣t∣≤1 时 ∣Tk(t)∣≤1 ; 当 ∣t∣>1 时 ∣Tk(t)∣>1 ;
(5) Tk(t)=0 的解为 ti=cos2k(2i−1)π,i=1,…,k;
(6) Tk(t) 有 n−1 个极值点: ti=cosnkπ,i=1,2,…,k−1 ;
下面的结论表明, 在所有首项系数为 1 的 k 次多项式中, pk(t)=2k−11Tk(t) 在 [−1,1] 上与零的偏差是最小的 (在无穷范数意义下).
定理6.31 设 pk(t)=2k−11Tk(t) , 则
−1≤t≤1max∣pk(t)∣≤−1≤t≤1max∣p(t)∣,∀p(t)∈P~k, 其中 P~k 表示所有首项系数为1的 k 次多项式组成的集合,即
∥pk(t)∥∞=p(t)∈Pkmin∥p(t)∥∞. 这里的范数 ∥⋅∥∞ 是 C[−1,1] 上的无穷范数,即 ∥f∥∞=maxx∈[−1,1]∣f(x)∣ , f(x)∈C[−1,1]
该性质可用于计算首项系数非零的 n 次多项式在 [−1,1] 上的 n−1 次最佳一致逼近多项式. 利用这个性质, 我们可以采用 Chebyshev 多项式的零点作为节点进行多项式插值, 以使得插值的总体误差达到最小化.
Chebyshev 另一个重要性质是下面的最小最大性质.
定理6.32 设 η∈R 满足 ∣η∣>1 ,则下面的最小最大问题
的唯一解为
p(t)∈Pk,p(η)=1min−1≤t≤1max∣p(t)∣ T~k(t)≜Tk(η)Tk(t). 通过简单的仿射变换, 该定理的结论可以推广到一般区间.
定理6.33 设 α,β,η∈R 满足 α<β 且 ∣η∣∈/[α,β] . 则下面的最小最大问题
p(t)∈Pk,p(η)=1minα≤x≤βmax∣p(t)∣ 的唯一解为
T^k(t)≜Tk(β−α2η−(β+α))Tk(β−α2t−(β+α)), 其中 Pk 表示所有次数不超过 k 的实系数多项式组成的集合.
Chebyshev加速方法
考虑迭代格式(6.16),我们假定:
(1) 迭代矩阵 G 的特征值都是实数;
(2) 迭代矩阵谱半径 ρ=ρ(G)<1 , 故 λ(G)∈[−ρ,ρ]⊂(−1,1) .
于是最小最大问题(6.21)就转化为
p∈Pk,p(1)=1minλ∈[−ρ,ρ]max{∣p(λ)∣}. 由于 1=[−ρ,ρ] , 根据定理 6.33, 上述问题的解为
pk(t)=Tk(1/ρ)Tk(t/ρ). 下面考虑 x~(k) 的计算. 我们无需先计算出 x(0),x(1),…,x(k) , 然后再通过线性组合 (6.18) 来计算 x~(k) . 事实上, 我们可以通过 Chebyshev 多项式的三项递推公式 (6.22), 由 x~(k−1) 和 x~(k−2) 直接计算出 x~(k) . 这样做的另一个好处是无需存储所有的 x~(i) . 下面给出具体的推导公式.
令 μk=Tk(1/ρ)1 即 Tk(1/ρ)=μk1. 由三项递推公式(6.22)可得
μk1=ρ2⋅μk−11−μk−21. 所以
x~(k)−x∗=pk(G)ε0=μkTk(G/ρ)ε0=μk[ρ2G⋅Tk−1(G/ρ)−Tk−2(G/ρ)]ε0=μk[ρ2G⋅μk−11pk−1(G)ε0−μk−21pk−2(G)ε0]=μk[ρ2G⋅μk−11(x~(k−1)−x∗)−μk−21(x~(k−2)−x∗)]. 整理后可得
x~(k)=μk−12μk⋅ρGx~(k−1)−μk−2μkx~(k−2)+dk, 其中
dk=x∗−μk−12μk⋅ρGx∗+μk−2μkx∗=x∗−μk−12μk⋅ρx∗−g+μk−2μkx∗=μk(μk1−ρμk−12+μk−21)x∗+μk−1ρ2μkg=μk−1ρ2μkg. 由此, 我们可以得到迭代格式 (6.16) 的 Chebyshev 加速方法
算法6.9.Chebyshev加速方法
1: Set μ0=1 , μ1=ρ=ρ(G) , x~(0)=x(0) , k=1
2: compute x~(1)=Gx(0)+g
3: while not converge do
4: k=k+1
5: μk=(ρ2⋅μk−11−μk−21)−1
6: x~(k)=μk−12μk⋅ρGx~(k−1)−μk−2μkx~(k−2)+μk−1ρ2μk⋅g
7: end while
该方法的每步迭代中只有一次矩阵向量乘积, 故方法每个迭代步的整体运算量与原迭代格式的每个迭代步的运算量基本相当.
设 λ(G)∈[α,β] ,且 −1<α≤β<1, 则我们也可以构造出相应的Chebyshev加速方法