标签：Chebyshev frac bm 漫谈 Krylov Ti y0 left lambda
漫谈：Chebyshev多项式，Krylov子空间，Chebyshev迭代，共轭梯度方法
Chebyshev迭代和共轭梯度方法的收敛速度（后者或称误差分析）都与Chebyshev多项式有着紧密联系，因此做一些整理，以期把其中的逻辑理清，推导理顺。只覆盖关键要点，不求面面俱到。
Chebyshev 多项式
标准Chebyshev多项式
与权函数
ρ
(
y
)
=
1
1
−
y
2
\rho(y)=\frac{1}{\sqrt{1-y^2}}
ρ(y)=1−y2
1对应的正交多项式为Chebyshev多项式。标准的Chebyshev多项式分两段定义：
T
n
(
y
)
=
{
c
o
s
(
n
arccos
(
y
)
)
,
y
∈
[
−
1
,
1
]
1
2
[
(
y
+
y
2
−
1
)
n
+
(
y
−
y
2
−
1
)
n
]
,
y
∈
R
∖
[
−
1
,
1
]
T_n(y)=\left\{ \begin{aligned} &\mathrm{cos}(n\space\arccos(y)) &\space, y\in [-1,1]\\ &\frac{1}{2}\left[\left(y+\sqrt{y^2-1}\right)^n+\left(y-\sqrt{y^2-1}\right)^n\right] &\space, y \in \mathbf{R} \setminus [-1,1] \end{aligned} \right.
Tn(y)=⎩⎨⎧cos(n arccos(y))21[(y+y2−1
)n+(y−y2−1
)n] ,y∈[−1,1] ,y∈R∖[−1,1]
概述一些基本性质。 n n n次的Chebyshev多项式在开区间 ( − 1 , 1 ) (-1,1) (−1,1)内有 n n n个零点，在 [ − 1 , 1 ] [-1,1] [−1,1]内有 n + 1 n+1 n+1个符号交错的极值点（ − 1 -1 −1和 + 1 +1 +1交替）。Chebyshev多项式不是奇函数就是偶函数，且随 n n n交替变化。
对于任意位置区间
[
a
,
b
]
[a,b]
[a,b]的Chebyshev多项式，都可以通过平移和缩放来化为标准区间
[
−
1
,
1
]
[-1,1]
[−1,1]上的标准Chebyshev多项式。对于
x
∈
[
a
,
b
]
x\in[a,b]
x∈[a,b]，做变换
y
=
x
−
1
2
(
b
+
a
)
1
2
(
b
−
a
)
y=\frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)}
y=21(b−a)x−21(b+a)，即有
y
∈
[
−
1
,
1
]
y\in[-1,1]
y∈[−1,1]：
T
n
~
(
x
)
=
T
n
(
y
=
x
−
1
2
(
b
+
a
)
1
2
(
b
−
a
)
)
\widetilde{T_n}(x)=T_n\left(y=\frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)}\right)
Tn
(x)=Tn(y=21(b−a)x−21(b+a))
相应地，权函数也做变换：
ρ
~
(
x
)
=
ρ
(
y
=
x
−
1
2
(
b
+
a
)
1
2
(
b
−
a
)
)
=
1
1
−
(
x
−
1
2
(
b
+
a
)
1
2
(
b
−
a
)
)
2
\widetilde{\rho}(x)=\rho\left(y=\frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)}\right)= \frac{1}{\sqrt{1-\left( \frac{x-\frac{1}{2}(b+a)}{\frac{1}{2}(b-a)} \right)^2}}
ρ
(x)=ρ(y=21(b−a)x−21(b+a))=1−(21(b−a)x−21(b+a))2
1
三项递推关系
记
θ
=
a
r
c
c
o
s
y
\theta=\mathrm{arccos}y
θ=arccosy，利用和差化积公式有
cos
(
(
n
+
1
)
θ
)
+
cos
(
(
n
−
1
)
θ
)
=
2
cos
(
n
θ
)
cos
(
θ
)
\cos((n+1)\theta)+\cos((n-1)\theta)=2\cos(n\theta)\cos(\theta)
cos((n+1)θ)+cos((n−1)θ)=2cos(nθ)cos(θ)
故有
T
n
+
1
(
y
)
=
2
y
T
n
(
y
)
−
T
n
−
1
(
y
)
T_{n+1}(y)=2yT_n(y)-T_{n-1}(y)
Tn+1(y)=2yTn(y)−Tn−1(y)
逼近性质
对之前提到的定义在
x
∈
[
a
,
b
]
x\in[a,b]
x∈[a,b]上的Chebyshev多项式，按照区间外一点
c
∉
[
a
,
b
]
c\notin[a,b]
c∈/[a,b]来做归一化：
T
n
^
(
x
)
=
T
n
~
(
x
)
T
n
~
(
c
)
\widehat{T_n}(x)=\frac{\widetilde{T_n}(x)}{\widetilde{T_n}(c)}
Tn
(x)=Tn
(c)Tn
(x)
显然有
T
n
^
(
x
)
∈
Φ
n
,
c
=
{
f
∈
P
∣
f
(
c
)
=
1
}
\widehat{T_n}(x)\in\Phi_{n,c}=\{f\in\mathbf{P}\space|\space f(c)=1 \}
Tn
(x)∈Φn,c={f∈P ∣ f(c)=1}。后者表示所有在
c
c
c这一点值为1的多项式的集合。Chebyshev多项式具有如下的逼近性质。
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
=
inf
f
∈
Φ
n
,
c
∣
∣
f
∣
∣
C
[
a
,
b
]
||\widehat{T_n}(x)||_{C[a,b]}={\underset {f\in\Phi_{n,c}}{\operatorname {inf} }} ||f||_{C[a,b]}
∣∣Tn
(x)∣∣C[a,b]=f∈Φn,cinf∣∣f∣∣C[a,b]
即所有在
c
c
c这一点值为1的多项式里，Chebyshev多项式是取到最小范数的那一个。由有限维空间内范数的等价性，从无穷范数入手证明。
首先说明 T n ^ ( x ) \widehat{T_n}(x) Tn (x)是well-defined的。因为 T n ( y ) T_n(y) Tn(y)的零点都在开区间 ( − 1 , 1 ) (-1,1) (−1,1)内，因此平移和缩放后 T n ~ ( x ) \widetilde{T_n}(x) Tn (x)的零点都在 [ a , b ] [a,b] [a,b]内， T n ~ ( c ) ≠ 0 \widetilde{T_n}(c)\neq0 Tn (c)=0。
反证法。假设 Y n ∈ Φ n , c Y_n\in\Phi_{n,c} Yn∈Φn,c，且满足 ∣ ∣ Y n ∣ ∣ C [ a , b ] < ∣ ∣ T n ^ ∣ ∣ C [ a , b ] ||Y_n||_{C[a,b]}<||\widehat{T_n}||_{C[a,b]} ∣∣Yn∣∣C[a,b]<∣∣Tn ∣∣C[a,b]。则对于 z n ( x ) = T n ^ ( x ) − Y n ( x ) z_n(x)=\widehat{T_n}(x)-Y_n(x) zn(x)=Tn (x)−Yn(x)，它有一个零点是 c c c。
对于标准Chebyshev多项式
T
n
y
T_n{y}
Tny，在
[
−
1
,
1
]
[-1,1]
[−1,1]上的（交错）极值点为
cos
(
k
π
n
)
,
k
=
0
,
1
,
.
.
.
,
n
\cos(\frac{k\pi}{n}), k=0,1,...,n
cos(nkπ),k=0,1,...,n。所以
T
n
~
(
x
)
\widetilde{T_n}(x)
Tn
(x)在
[
a
,
b
]
[a,b]
[a,b]上的（交错）极值点为
x
k
=
b
+
a
2
+
b
−
a
2
cos
(
k
π
n
)
,
k
=
0
,
1
,
.
.
.
,
n
x_k=\frac{b+a}{2}+\frac{b-a}{2}\cos(\frac{k\pi}{n}), k=0,1,...,n
xk=2b+a+2b−acos(nkπ),k=0,1,...,n。因为
∣
∣
Y
n
(
x
)
∣
∣
C
[
a
,
b
]
<
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
||Y_n(x)||_{C[a,b]}<||\widehat{T_n}(x)||_{C[a,b]}
∣∣Yn(x)∣∣C[a,b]<∣∣Tn
(x)∣∣C[a,b]，那么自然有
∣
Y
n
(
x
k
)
∣
≤
∣
∣
Y
n
(
x
)
∣
∣
C
[
a
,
b
]
<
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
=
∣
T
n
^
(
x
k
)
∣
,
k
=
0
,
1
,
.
.
.
,
n
|Y_n(x_k)|\leq||Y_n(x)||_{C[a,b]}<||\widehat{T_n}(x)||_{C[a,b]}=|\widehat{T_n}(x_k)|\space, k=0,1,...,n
∣Yn(xk)∣≤∣∣Yn(x)∣∣C[a,b]<∣∣Tn
(x)∣∣C[a,b]=∣Tn
(xk)∣ ,k=0,1,...,n
所以
T
n
^
(
x
k
)
\widehat{T_n}(x_k)
Tn
(xk)与
T
n
^
(
x
k
)
−
Y
n
(
x
k
)
\widehat{T_n}(x_k)-Y_n(x_k)
Tn
(xk)−Yn(xk)同号！而
{
x
k
}
\{x_k\}
{xk}是
T
n
^
(
x
)
\widehat{T_n}(x)
Tn
(x)的一组符号交错点，那么也是
T
n
^
(
x
)
−
Y
n
(
x
)
\widehat{T_n}(x)-Y_n(x)
Tn
(x)−Yn(x)的一组符号交错点。所以在开区间
(
x
k
,
x
k
+
1
)
(x_k,x_{k+1})
(xk,xk+1)中必至少有
z
n
(
x
)
z_n(x)
zn(x)的一个零点，则
(
a
,
b
)
(a,b)
(a,b)内至少有
z
n
(
x
)
z_n(x)
zn(x)的
n
n
n个不同零点。
加上前述的 c ∉ [ a , b ] c\notin[a,b] c∈/[a,b]也是 z n ( x ) z_n(x) zn(x)的一个零点，故它至少有 n + 1 n+1 n+1个零点，这与它是 n n n次多项式矛盾。证毕。
再来估计
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
=
inf
f
∈
Φ
n
,
c
∣
∣
f
∣
∣
C
[
a
,
b
]
||\widehat{T_n}(x)||_{C[a,b]}={\underset {f\in\Phi_{n,c}}{\operatorname {inf} }} ||f||_{C[a,b]}
∣∣Tn
(x)∣∣C[a,b]=f∈Φn,cinf∣∣f∣∣C[a,b]的大小。
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
=
∣
∣
T
n
~
(
x
)
T
n
~
(
c
)
∣
∣
C
[
a
,
b
]
=
∣
∣
T
n
~
(
x
)
∣
∣
C
[
a
,
b
]
∣
T
n
~
(
c
)
∣
=
1
∣
T
n
~
(
c
)
∣
||\widehat{T_n}(x)||_{C[a,b]}=\left|\left|\frac{\widetilde{T_n}(x)}{\widetilde{T_n}(c)}\right|\right|_{C[a,b]}= \frac{||\widetilde{T_n}(x)||_{C[a,b]}}{|\widetilde{T_n}(c)|}= \frac{1}{|\widetilde{T_n}(c)|}
∣∣Tn
(x)∣∣C[a,b]=∣∣∣∣∣∣∣∣∣∣Tn
(c)Tn
(x)∣∣∣∣∣∣∣∣∣∣C[a,b]=∣Tn
(c)∣∣∣Tn
(x)∣∣C[a,b]=∣Tn
(c)∣1
而
T
n
~
(
c
)
=
T
n
(
c
−
b
+
a
2
b
−
a
2
)
=
T
n
(
−
(
b
−
c
)
+
(
a
−
c
)
2
(
b
−
c
)
−
(
a
−
c
)
2
)
=
(
−
1
)
n
T
n
(
(
b
−
c
)
+
(
a
−
c
)
(
b
−
c
)
−
(
a
−
c
)
)
\widetilde{T_n}(c)= T_n\left(\frac{c-\frac{b+a}{2}}{\frac{b-a}{2}} \right)= T_n\left(\frac{-\frac{(b-c)+(a-c)}{2}}{\frac{(b-c)-(a-c)}{2}} \right)= (-1)^nT_n\left( \frac{(b-c)+(a-c)}{(b-c)-(a-c)} \right)
Tn
(c)=Tn(2b−ac−2b+a)=Tn(2(b−c)−(a−c)−2(b−c)+(a−c))=(−1)nTn((b−c)−(a−c)(b−c)+(a−c))
所以
∣
T
n
~
(
c
)
∣
=
T
n
(
∣
(
b
−
c
)
+
(
a
−
c
)
(
b
−
c
)
−
(
a
−
c
)
∣
)
=
T
n
(
λ
+
1
λ
−
1
)
|\widetilde{T_n}(c)|= T_n\left(\left| \frac{(b-c)+(a-c)}{(b-c)-(a-c)} \right|\right)= T_n\left(\frac{\lambda+1}{\lambda-1} \right)
∣Tn
(c)∣=Tn(∣∣∣∣(b−c)−(a−c)(b−c)+(a−c)∣∣∣∣)=Tn(λ−1λ+1)
其中
λ
=
max
(
∣
a
−
c
∣
∣
b
−
c
∣
,
∣
b
−
c
∣
∣
a
−
c
∣
)
>
1
\lambda=\max(\frac{|a-c|}{|b-c|},\frac{|b-c|}{|a-c|})>1
λ=max(∣b−c∣∣a−c∣,∣a−c∣∣b−c∣)>1。所以根据标准Chebyshev多项式在
[
−
1
,
1
]
[-1,1]
[−1,1]区间外的表达式，有
T
n
(
λ
+
1
λ
−
1
)
≥
1
2
[
λ
+
1
λ
−
1
+
(
λ
+
1
λ
−
1
)
2
−
1
]
n
=
1
2
(
λ
+
1
λ
−
1
)
n
T_n\left(\frac{\lambda+1}{\lambda-1} \right)\geq \frac{1}{2}\left[ \frac{\lambda+1}{\lambda-1}+\sqrt{ \left(\frac{\lambda+1}{\lambda-1}\right)^2-1 } \right]^n= \frac{1}{2}\left(\frac{\sqrt{\lambda}+1}{\sqrt{\lambda}-1}\right)^n
Tn(λ−1λ+1)≥21⎣⎡λ−1λ+1+(λ−1λ+1)2−1
⎦⎤n=21(λ
−1λ
+1)n
所以最后得到的估计是
∣
∣
T
n
^
(
x
)
∣
∣
C
[
a
,
b
]
=
1
∣
T
n
~
(
c
)
∣
=
1
T
n
(
λ
+
1
λ
−
1
)
≤
2
(
λ
−
1
λ
+
1
)
n
||\widehat{T_n}(x)||_{C[a,b]}= \frac{1}{|\widetilde{T_n}(c)|}= \frac{1}{T_n\left(\frac{\lambda+1}{\lambda-1} \right)}\leq 2\left(\frac{\sqrt{\lambda}-1}{\sqrt{\lambda}+1}\right)^n
∣∣Tn
(x)∣∣C[a,b]=∣Tn
(c)∣1=Tn(λ−1λ+1)1≤2(λ
+1λ
−1)n
并且这个估计是较紧凑的。将会在Chebyshev迭代和共轭梯度的收敛速度分析中用到这个结论。
Chebyshev迭代
Chebyshev迭代是线性非定常迭代，其迭代使用的矩阵/系数在每个迭代步都会发生变化。而任意的单步线性定常迭代都可以表述为最简单的Richardson迭代+预处理的组合。因此从Richardson迭代出发，引入Krylov子空间的概念，然后加以适当改造得到Chebyshev迭代。
Richardson迭代和Krylov子空间
Richardson迭代是最简单的迭代格式，它在每一步都做一个修正量，此修正量就简单地取为当前迭代结果的残差：
x
(
i
+
1
)
=
x
(
i
)
+
r
(
i
)
=
x
(
i
)
+
b
−
A
x
(
i
)
=
(
I
−
A
)
x
(
i
)
+
b
\bm{x}^{(i+1)}=\bm{x}^{(i)}+\bm{r}^{(i)}=\bm{x}^{(i)}+\bm{b}-\bm{A}\bm{x}^{(i)}=(\bm{I}-\bm{A})\bm{x}^{(i)}+\bm{b}
x(i+1)=x(i)+r(i)=x(i)+b−Ax(i)=(I−A)x(i)+b
r
(
i
+
1
)
=
b
−
A
x
(
i
+
1
)
=
b
−
A
x
(
i
)
−
A
r
(
i
)
=
(
I
−
A
)
r
(
i
)
=
.
.
.
=
(
I
−
A
)
i
+
1
r
(
0
)
\bm{r}^{(i+1)}=\bm{b}-\bm{A}\bm{x}^{(i+1)}=\bm{b}-\bm{A}\bm{x}^{(i)}-\bm{A}\bm{r}^{(i)} =(\bm{I}-\bm{A})\bm{r}^{(i)}=...=(\bm{I}-\bm{A})^{i+1}\bm{r}^{(0)}
r(i+1)=b−Ax(i+1)=b−Ax(i)−Ar(i)=(I−A)r(i)=...=(I−A)i+1r(0)
由残差与误差的关系有
e
(
i
)
=
A
−
1
r
(
i
)
\bm{e}^{(i)}=\bm{A}^{-1}\bm{r}^{(i)}
e(i)=A−1r(i)，可知相应有
e
(
i
+
1
)
=
(
I
−
A
)
e
(
i
)
\bm{e}^{(i+1)}=(\bm{I}-\bm{A})\bm{e}^{(i)}
e(i+1)=(I−A)e(i)。于是在第
m
m
m步得到的迭代结果为
x
(
m
)
=
x
(
0
)
+
∑
i
=
0
m
−
1
r
(
i
)
=
x
(
0
)
+
∑
i
=
0
m
−
1
(
I
−
A
)
i
r
(
0
)
∈
x
(
0
)
+
K
m
(
r
(
0
)
)
\bm{x}^{(m)}=\bm{x}^{(0)}+\sum_{i=0}^{m-1}\bm{r}^{(i)}=\bm{x}^{(0)}+\sum_{i=0}^{m-1}(\bm{I}-\bm{A})^{i}\bm{r}^{(0)}\in\bm{x}^{(0)}+\bm{K_m}(\bm{r}^{(0)})
x(m)=x(0)+i=0∑m−1r(i)=x(0)+i=0∑m−1(I−A)ir(0)∈x(0)+Km(r(0))
其中
K
m
(
r
(
0
)
)
\bm{K_m}(\bm{r}^{(0)})
Km(r(0))是由
r
(
0
)
\bm{r}^{(0)}
r(0)生成的Krylov子空间，其中的元素是用
A
\bm{A}
A反复作用于
r
(
0
)
\bm{r}^{(0)}
r(0)得到的
r
(
0
)
,
A
r
(
0
)
,
A
2
r
(
0
)
,
.
.
.
\bm{r}^{(0)},\bm{Ar}^{(0)},\bm{A^2r}^{(0)},...
r(0),Ar(0),A2r(0),...等的线性组合。
任何的单步定常迭代格式，比如给定
A
=
M
−
N
\bm{A}=\bm{M}-\bm{N}
A=M−N，其中
M
\bm{M}
M为预处理矩阵，则迭代可以写成
x
(
i
+
1
)
=
M
−
1
N
x
(
i
)
+
M
−
1
b
=
M
−
1
(
M
−
A
)
x
(
i
)
+
M
−
1
b
\bm{x}^{(i+1)}=\bm{M}^{-1}\bm{N}\bm{x}^{(i)}+\bm{M}^{-1}\bm{b}=\bm{M}^{-1}(\bm{M}-\bm{A})\bm{x}^{(i)}+\bm{M}^{-1}\bm{b}
x(i+1)=M−1Nx(i)+M−1b=M−1(M−A)x(i)+M−1b
而对
A
x
=
b
\bm{Ax}=\bm{b}
Ax=b做预处理后得到的
M
−
1
A
x
=
M
−
1
b
\bm{M^{-1}Ax}=\bm{M^{-1}b}
M−1Ax=M−1b，再进行Richardson迭代，可以得到
x
(
i
+
1
)
=
x
(
i
)
+
r
(
i
)
=
x
(
i
)
+
M
−
1
b
−
M
−
1
A
x
(
i
)
\bm{x}^{(i+1)}=\bm{x}^{(i)}+\bm{r}^{(i)}=\bm{x}^{(i)}+\bm{M^{-1}b}-\bm{M^{-1}}\bm{A}\bm{x}^{(i)}
x(i+1)=x(i)+r(i)=x(i)+M−1b−M−1Ax(i)
可见是等价的，所以任何的单步定常迭代形式都可以写成预处理后的Richardson迭代（即当前迭代值+修正量）的形式。
单步迭代格式
原始的Richardson迭代是最简单地取了
M
=
I
\bm{M}=\bm{I}
M=I（最好计算，但最不好近似）来做预处理。现在采用非定常，即考虑每步迭代不那么naive，而是选用一个与当前步
i
i
i有关地矩阵
M
i
\bm{M_i}
Mi做预处理，假设
M
i
=
τ
i
M
\bm{M_i}=\tau_i\bm{M}
Mi=τiM由同一个矩阵
M
\bm{M}
M生成。同样类似Richardson迭代的好计算的想法，取
M
=
I
\bm{M}=\bm{I}
M=I，故有
M
i
=
τ
i
I
\bm{M_i}=\tau_i\bm{I}
Mi=τiI。所以残量的递推关系成为
r
(
i
+
1
)
=
(
I
−
τ
i
A
)
r
(
i
)
=
(
I
−
τ
i
A
)
(
I
−
τ
i
−
1
A
)
.
.
.
(
I
−
τ
0
A
)
r
(
0
)
\bm{r}^{(i+1)}=(\bm{I}-\tau_i\bm{A})\bm{r}^{(i)}=(\bm{I}-\tau_i\bm{A})(\bm{I}-\tau_{i-1}\bm{A})...(\bm{I}-\tau_0\bm{A})\bm{r}^{(0)}
r(i+1)=(I−τiA)r(i)=(I−τiA)(I−τi−1A)...(I−τ0A)r(0)
对于误差
e
(
i
)
\bm{e}^{(i)}
e(i)也同样有此关系，所以有
∣
∣
e
(
i
)
∣
∣
∣
∣
e
(
0
)
∣
∣
≤
∣
∣
(
I
−
τ
i
−
1
A
)
.
.
.
(
I
−
τ
0
A
)
∣
∣
\frac{||\bm{e}^{(i)}||}{||\bm{e}^{(0)}||}\leq|| (\bm{I}-\tau_{i-1}\bm{A})...(\bm{I}-\tau_0\bm{A}) ||
∣∣e(0)∣∣∣∣e(i)∣∣≤∣∣(I−τi−1A)...(I−τ0A)∣∣
我们希望构造出来的结果（即
τ
0
,
τ
1
,
.
.
.
\tau_0,\tau_1,...
τ0,τ1,...等的取值）能使残量（误差）收缩得最小，即求解如下的优化问题：
inf
τ
0
,
.
.
.
,
τ
i
−
1
∣
∣
(
I
−
τ
i
−
1
A
)
.
.
.
(
I
−
τ
0
A
)
∣
∣
{\underset {\tau_0,...,\tau_{i-1}} {\operatorname {inf} }} || (\bm{I}-\tau_{i-1}\bm{A})...(\bm{I}-\tau_0\bm{A}) ||
τ0,...,τi−1inf∣∣(I−τi−1A)...(I−τ0A)∣∣
假设
A
\bm{A}
A是对称阵（这些强加的适用条件后面汇总来看），那么存在正交阵
P
\bm{P}
P（如果
A
\bm{A}
A是Hermite矩阵，则
P
\bm{P}
P为酉矩阵）使得
A
\bm{A}
A能相似对角化，且取2-范数时有
∣
∣
P
∣
∣
2
=
1
||\bm{P}||_2=1
∣∣P∣∣2=1。
A
=
P
d
i
a
g
(
λ
1
,
λ
2
,
.
.
.
,
λ
n
)
P
T
=
P
Λ
P
T
\bm{A}=\bm{P} \mathrm{diag}(\lambda_1,\lambda_2,...,\lambda_n)\bm{P}^T=\bm{P}\bm{\Lambda} \bm{P}^T
A=Pdiag(λ1,λ2,...,λn)PT=PΛPT
那么有
∣
∣
(
I
−
τ
i
−
1
A
)
.
.
.
(
I
−
τ
0
A
)
∣
∣
2
=
∣
∣
P
(
I
−
τ
i
−
1
Λ
)
.
.
.
(
I
−
τ
0
Λ
)
P
T
∣
∣
2
=
∣
∣
(
I
−
τ
i
−
1
Λ
)
.
.
.
(
I
−
τ
0
Λ
)
∣
∣
2
=
max
λ
∈
σ
(
A
)
∣
(
1
−
τ
i
−
1
λ
)
.
.
.
(
1
−
τ
0
λ
)
∣
\begin{aligned} || (\bm{I}-\tau_{i-1}\bm{A})...(\bm{I}-\tau_0\bm{A}) ||_2=& ||\bm{P}(\bm{I}-\tau_{i-1}\bm{\Lambda})...(\bm{I}-\tau_0\bm{\Lambda})\bm{P}^T||_2\\ =&||(\bm{I}-\tau_{i-1}\bm{\Lambda})...(\bm{I}-\tau_0\bm{\Lambda})||_2\\ =&{\underset {\lambda\in\sigma(\bm{A})} {\operatorname {max} }} | (1-\tau_{i-1}\lambda)...(1-\tau_0\lambda)| \\ \end{aligned}
∣∣(I−τi−1A)...(I−τ0A)∣∣2===∣∣P(I−τi−1Λ)...(I−τ0Λ)PT∣∣2∣∣(I−τi−1Λ)...(I−τ0Λ)∣∣2λ∈σ(A)max∣(1−τi−1λ)...(1−τ0λ)∣
而
A
\bm{A}
A的谱是可以估计的，假设已算出
σ
(
A
)
∈
[
λ
‾
,
λ
ˉ
]
\sigma(\bm{A})\in[\underline{\lambda},\bar{\lambda}]
σ(A)∈[λ,λˉ]。那么上式可以写成：
max
λ
∈
σ
(
A
)
∣
(
1
−
τ
i
−
1
λ
)
.
.
.
(
1
−
τ
0
λ
)
∣
≤
∣
∣
(
1
−
τ
i
−
1
λ
)
.
.
.
(
1
−
τ
0
λ
)
∣
∣
∞
,
[
λ
‾
,
λ
ˉ
]
=
∣
∣
f
(
λ
)
∣
∣
∞
,
[
λ
‾
,
λ
ˉ
]
\begin{aligned} {\underset {\lambda\in\sigma(\bm{A})} {\operatorname {max} }} | (1-\tau_{i-1}\lambda)...(1-\tau_0\lambda)|&\leq& ||(1-\tau_{i-1}\lambda)...(1-\tau_0\lambda)||_{\infty,[\underline{\lambda},\bar{\lambda}]}\\ &=&||f(\lambda)||_{\infty,[\underline{\lambda},\bar{\lambda}]}\\ \end{aligned}
λ∈σ(A)max∣(1−τi−1λ)...(1−τ0λ)∣≤=∣∣(1−τi−1λ)...(1−τ0λ)∣∣∞,[λ,λˉ]∣∣f(λ)∣∣∞,[λ,λˉ]
左边的原式是一个绝对值，现在右边看成是一个关于
λ
\lambda
λ的多项式函数
f
(
λ
)
=
(
1
−
τ
i
−
1
λ
)
.
.
.
(
1
−
τ
0
λ
)
f(\lambda)=(1-\tau_{i-1}\lambda)...(1-\tau_0\lambda)
f(λ)=(1−τi−1λ)...(1−τ0λ)的无穷范数。该多项式函数有
i
i
i个实单根：
1
τ
0
,
.
.
.
,
1
τ
i
−
1
\frac{1}{\tau_0},...,\frac{1}{\tau_{i-1}}
τ01,...,τi−11。
而由之前的结论，定义在 [ λ ‾ , λ ˉ ] [\underline{\lambda},\bar{\lambda}] [λ,λˉ]上的Chebyshev多项式 T i ^ ( λ ) \widehat{T_i}(\lambda) Ti (λ)在开区间 ( λ ‾ , λ ˉ ) (\underline{\lambda},\bar{\lambda}) (λ,λˉ)内有 i i i个实单根，且 T i ^ ( λ ) = inf P i ∈ P i , P i ( 0 ) = 1 ∣ ∣ P i ∣ ∣ ∞ , [ λ ‾ , λ ˉ ] \widehat{T_i}(\lambda)={\underset {P_i\in\mathbf{P_i},P_i(0)=1}{\operatorname {inf} }} ||P_i||_{\infty,[\underline{\lambda},\bar{\lambda}]} Ti (λ)=Pi∈Pi,Pi(0)=1inf∣∣Pi∣∣∞,[λ,λˉ].
注意当
A
\bm{A}
A为对称正定阵时，
0
∉
[
λ
‾
,
λ
ˉ
]
0\notin[\underline{\lambda},\bar{\lambda}]
0∈/[λ,λˉ]，所以
inf
P
i
∈
P
i
,
P
i
(
0
)
=
1
∣
∣
P
i
∣
∣
∞
,
[
λ
‾
,
λ
ˉ
]
{\underset {P_i\in\mathbf{P_i},P_i(0)=1}{\operatorname {inf} }} ||P_i||_{\infty,[\underline{\lambda},\bar{\lambda}]}
Pi∈Pi,Pi(0)=1inf∣∣Pi∣∣∞,[λ,λˉ]的下界恰好能被原问题的
f
(
λ
)
f(\lambda)
f(λ)取到。此时
P
i
∗
(
λ
)
=
T
i
^
(
λ
)
=
T
i
~
(
λ
)
T
i
~
(
0
)
=
T
i
(
λ
−
1
2
(
λ
ˉ
+
λ
‾
)
1
2
(
λ
ˉ
−
λ
‾
)
)
T
i
(
0
−
1
2
(
λ
ˉ
+
λ
‾
)
1
2
(
λ
ˉ
−
λ
‾
)
)
=
(
1
−
τ
i
−
1
(
i
)
λ
)
.
.
.
(
1
−
τ
0
(
i
)
λ
)
P_i^*(\lambda)=\widehat{T_i}(\lambda)=\frac{\widetilde{T_i}(\lambda)}{\widetilde{T_i}(0)} =\frac{T_i\left( \frac{\lambda-\frac{1}{2}(\bar{\lambda}+\underline{\lambda})} {\frac{1}{2} (\bar{\lambda}-\underline{\lambda})} \right)}{T_i\left( \frac{0-\frac{1}{2}(\bar{\lambda}+\underline{\lambda})} {\frac{1}{2} (\bar{\lambda}-\underline{\lambda})} \right)}= (1-\tau_{i-1}^{(i)}\lambda)...(1-\tau_{0}^{(i)}\lambda)
Pi∗(λ)=Ti
(λ)=Ti
(0)Ti
(λ)=Ti(21(λˉ−λ)0−21(λˉ+λ))Ti(21(λˉ−λ)λ−21(λˉ+λ))=(1−τi−1(i)λ)...(1−τ0(i)λ)
其中
τ
j
(
i
)
\tau_j^{(i)}
τj(i)的上标表示与迭代次数
i
i
i有关。所以原问题的
i
i
i个实单根就是定义在
[
λ
‾
,
λ
ˉ
]
[\underline{\lambda},\bar{\lambda}]
[λ,λˉ]上的Chebyshev多项式的根。单步的迭代格式写成
x
(
1
)
=
x
(
0
)
+
τ
0
(
i
)
r
(
0
)
.
.
.
x
(
j
)
=
x
(
j
−
1
)
+
τ
j
−
1
(
i
)
r
(
j
−
1
)
.
.
.
x
(
i
)
=
x
(
i
−
1
)
+
τ
i
−
1
(
i
)
r
(
i
−
1
)
\begin{aligned} \bm{x}^{(1)}&=&\bm{x}^{(0)}+\tau_0^{(i)}\bm{r}^{(0)}\\ &...&\\ \bm{x}^{(j)}&=&\bm{x}^{(j-1)}+\tau_{j-1}^{(i)}\bm{r}^{(j-1)}\\ &...&\\ \bm{x}^{(i)}&=&\bm{x}^{(i-1)}+\tau_{i-1}^{(i)}\bm{r}^{(i-1)}\\ \end{aligned}
x(1)x(j)x(i)=...=...=x(0)+τ0(i)r(0)x(j−1)+τj−1(i)r(j−1)x(i−1)+τi−1(i)r(i−1)
收敛速度
Chebyshev迭代的收敛速度优势，应当通过进行一次 i i i步的Chebyshev迭代，和进行 i i i次的单步定常迭代来比较。
由Chebyshev的逼近性质，可以得到Chebyshev的迭代收敛速度为
∣
∣
e
(
i
)
∣
∣
∣
∣
e
(
0
)
∣
∣
≤
∣
∣
(
1
−
τ
i
−
1
λ
)
.
.
.
(
1
−
τ
0
λ
)
∣
∣
∞
,
[
λ
‾
,
λ
ˉ
]
=
1
∣
T
i
~
(
0
)
∣
=
1
∣
T
i
(
0
−
λ
ˉ
−
λ
‾
2
λ
ˉ
−
λ
‾
2
)
∣
=
1
∣
T
i
(
λ
ˉ
λ
‾
+
1
λ
ˉ
λ
‾
−
1
)
∣
=
1
∣
T
i
(
α
+
1
α
−
1
)
∣
≤
≈
2
(
α
−
1
α
+
1
)
i
\begin{aligned} \frac{||\bm{e}^{(i)}||}{||\bm{e}^{(0)||}}&\leq&||(1-\tau_{i-1}\lambda)...(1-\tau_0\lambda)||_{\infty,[\underline{\lambda},\bar{\lambda}]}\\ &=&\frac{1}{|\widetilde{T_i}(0)|}=\frac{1}{\left|T_i\left(\frac{0-\frac{\bar{\lambda}-\underline{\lambda}}{2}}{\frac{\bar{\lambda}-\underline{\lambda}}{2}}\right)\right|}\\ &=&\frac{1}{\left| T_i\left( \frac{\frac{\bar{\lambda}}{\underline{\lambda}}+1}{\frac{\bar{\lambda}}{\underline{\lambda}}-1} \right) \right|}=\frac{1}{\left| T_i\left( \frac{\alpha+1}{\alpha-1} \right) \right|}\\ &&\leq\approx2\left( \frac{\sqrt{\alpha}-1}{\sqrt{\alpha}+1} \right)^i \end{aligned}
∣∣e(0)∣∣∣∣e(i)∣∣≤==∣∣(1−τi−1λ)...(1−τ0λ)∣∣∞,[λ,λˉ]∣Ti
(0)∣1=∣∣∣∣Ti(2λˉ−λ0−2λˉ−λ)∣∣∣∣1∣∣∣∣Ti(λλˉ−1λλˉ+1)∣∣∣∣1=∣∣Ti(α−1α+1)∣∣1≤≈2(α
+1α
−1)i
其中
α
=
λ
ˉ
λ
‾
>
1
\alpha=\frac{\bar{\lambda}}{\underline{\lambda}}>1
α=λλˉ>1。
而如果做
i
i
i次单步的定常迭代，每次单步的收敛速度为
P
1
∗
(
λ
)
=
1
∣
T
1
(
α
+
1
α
−
1
)
∣
=
1
∣
α
+
1
α
−
1
∣
=
α
−
1
α
+
1
P_1^*(\lambda)=\frac{1}{\left| T_1\left( \frac{\alpha+1}{\alpha-1} \right) \right|} = \frac{1}{\left| \frac{\alpha+1}{\alpha-1} \right|} = \frac{\alpha-1}{\alpha+1}
P1∗(λ)=∣∣T1(α−1α+1)∣∣1=∣∣α−1α+1∣∣1=α+1α−1
连续做
i
i
i次以后的收敛速度为
∣
∣
e
(
i
)
∣
∣
∣
∣
e
(
0
)
∣
∣
≤
(
α
−
1
α
+
1
)
i
\frac{||\bm{e}^{(i)}||}{||\bm{e}^{(0)||}}\leq\left( \frac{\alpha-1}{\alpha+1}\right)^i
∣∣e(0)∣∣∣∣e(i)∣∣≤(α+1α−1)i。这显然要比Chebyshev迭代的收敛速度慢。
两步迭代格式
单步的Chebyshev迭代格式的缺点是 τ j ( i ) \tau_j^{(i)} τj(i)依赖于迭代次数 i i i的取值。自然地想，能不能不依赖于后者的确定？这就可以用上Chebyshev多项式的三项递推关系来解决！
对标准的Chebyshev多项式，
T
n
+
1
(
y
)
=
2
y
T
n
(
y
)
−
T
n
−
1
(
y
)
,
.
.
.
,
T
1
(
y
)
=
y
,
T
0
(
y
)
=
1
T_{n+1}(y)=2yT_n(y)-T_{n-1}(y),...,T_1(y)=y,T_0(y)=1
Tn+1(y)=2yTn(y)−Tn−1(y),...,T1(y)=y,T0(y)=1。但需要的是定义在
[
λ
‾
,
λ
ˉ
]
[\underline{\lambda}, \bar{\lambda}]
[λ,λˉ]上的Chebyshev多项式，做变换：
T
i
~
(
x
)
=
T
i
(
y
=
x
−
1
2
(
λ
ˉ
+
λ
‾
)
1
2
(
λ
ˉ
−
λ
‾
)
)
=
T
i
(
y
=
a
x
−
y
0
)
\widetilde{T_i}(x)=T_i\left(y=\frac{x-\frac{1}{2}(\bar{\lambda}+\underline{\lambda})}{\frac{1}{2}(\bar{\lambda}-\underline{\lambda})}\right)=T_i(y=ax-y_0)
Ti
(x)=Ti(y=21(λˉ−λ)x−21(λˉ+λ))=Ti(y=ax−y0)
其中记
y
0
=
λ
ˉ
+
λ
‾
λ
ˉ
−
λ
‾
,
a
=
2
λ
ˉ
−
λ
‾
y_0=\frac{\bar{\lambda}+\underline{\lambda}}{\bar{\lambda}-\underline{\lambda}}, a=\frac{2}{\bar{\lambda}-\underline{\lambda}}
y0=λˉ−λλˉ+λ,a=λˉ−λ2。按照
x
=
0
x=0
x=0这一点的
T
i
~
\widetilde{T_i}
Ti
值做归一化：
T
i
^
(
x
)
=
T
i
~
(
x
)
T
i
~
(
0
)
=
T
i
(
a
x
−
y
0
)
T
i
(
−
y
0
)
=
T
i
(
y
0
−
a
x
)
T
i
(
y
0
)
\widehat{T_i}(x)=\frac{\widetilde{T_i}(x)}{\widetilde{T_i}(0)}=\frac{T_i(ax-y_0)}{T_i(-y_0)}=\frac{T_i(y_0-ax)}{T_i(y_0)}
Ti
(x)=Ti
(0)Ti
(x)=Ti(−y0)Ti(ax−y0)=Ti(y0)Ti(y0−ax)，因此三项递推关系为
T
0
^
(
x
)
=
T
0
(
y
0
−
a
x
)
T
0
(
y
0
)
=
1
T
1
^
(
x
)
=
T
1
(
y
0
−
a
x
)
T
1
(
y
0
)
=
y
0
−
a
x
y
0
=
T
0
^
(
x
)
−
a
y
0
x
T
0
^
(
x
)
.
.
.
.
.
.
T
i
+
1
^
(
x
)
=
T
i
+
1
(
y
0
−
a
x
)
T
i
+
1
(
y
0
)
=
2
(
y
0
−
a
x
)
T
i
(
y
0
−
a
x
)
−
T
i
−
1
(
y
0
−
a
x
)
T
i
+
1
(
y
0
)
=
2
(
y
0
−
a
x
)
T
i
(
y
0
)
T
i
+
1
(
y
0
)
T
i
(
y
0
−
a
x
)
T
i
(
y
0
)
−
T
i
−
1
(
y
0
)
T
i
+
1
(
y
0
)
T
i
−
1
(
y
0
−
a
x
)
T
i
−
1
(
y
0
)
=
2
(
y
0
−
a
x
)
T
i
(
y
0
)
T
i
+
1
(
y
0
)
T
i
^
(
x
)
−
T
i
−
1
(
y
0
)
T
i
+
1
(
y
0
)
T
i
−
1
^
(
x
)
=
2
y
0
T
i
(
y
0
)
T
i
+
1
(
y
0
)
T
i
^
(
x
)
−
T
i
−
1
(
y
0
)
T
i
+
1
(
y
0
)
T
i
−
1
^
(
x
)
−
2
a
T
i
(
y
0
)
T
i
+
1
(
y
0
)
x
T
i
^
(
x
)
=
(
1
+
T
i
−
1
(
y
0
)
T
i
+
1
(
y
0
)
)
T
i
^
(
x
)
−
T
i
−
1
(
y
0
)
T
i
+
1
(
y
0
)
T
i
−
1
^
(
x
)
−
2
a
T
i
(
y
0
)
T
i
+
1
(
y
0
)
x
T
i
^
(
x
)
=
α
i
T
i
^
(
x
)
+
(
1
−
α
i
)
T
i
−
1
^
(
x
)
−
β
i
x
T
i
^
(
x
)
\widehat{T_0}(x)=\frac{T_0(y_0-ax)}{T_0(y_0)}=1\\ \widehat{T_1}(x)=\frac{T_1(y_0-ax)}{T_1(y_0)}=\frac{y_0-ax}{y_0}=\widehat{T_0}(x)-\frac{a}{y_0}x\widehat{T_0}(x)\\ ... ...\\ \begin{aligned} \widehat{T_{i+1}}(x)&=\frac{T_{i+1}(y_0-ax)}{T_{i+1}(y_0)}=\frac{2(y_0-ax)T_i(y_0-ax)-T_{i-1}(y_0-ax)}{T_{i+1}(y_0)}\\ &=\frac{2(y_0-ax)T_i(y_0)}{T_{i+1}(y_0)} \frac{T_i(y_0-ax)}{T_i(y_0)}-\frac{T_{i-1}(y_0)}{T_{i+1}(y_0)} \frac{T_{i-1}(y_0-ax)}{T_{i-1}(y_0)}\\ &=\frac{2(y_0-ax)T_i(y_0)}{T_{i+1}(y_0)} \widehat{T_i}(x) -\frac{T_{i-1}(y_0)}{T_{i+1}(y_0)} \widehat{T_{i-1}}(x)\\ &=\frac{2y_0T_i(y_0)}{T_{i+1}(y_0)} \widehat{T_i}(x) -\frac{T_{i-1}(y_0)}{T_{i+1}(y_0)} \widehat{T_{i-1}}(x)-2a\frac{T_i(y_0)}{T_{i+1}(y_0)}x\widehat{T_i}(x)\\ &=\left( 1+\frac{T_{i-1}(y_0)}{T_{i+1}(y_0)} \right) \widehat{T_i}(x) - \frac{T_{i-1}(y_0)}{T_{i+1}(y_0)}\widehat{T_{i-1}}(x) -2a\frac{T_i(y_0)}{T_{i+1}(y_0)}x\widehat{T_i}(x)\\ &=\alpha_i\widehat{T_i}(x) + (1-\alpha_i)\widehat{T_{i-1}}(x) - \beta_ix\widehat{T_i}(x) \end{aligned}
T0
(x)=T0(y0)T0(y0−ax)=1T1
(x)=T1(y0)T1(y0−ax)=y0y0−ax=T0
(x)−y0axT0
(x)......Ti+1
(x)=Ti+1(y0)Ti+1(y0−ax)=Ti+1(y0)2(y0−ax)Ti(y0−ax)−Ti−1(y0−ax)=Ti+1(y0)2(y0−ax)Ti(y0)Ti(y0)Ti(y0−ax)−Ti+1(y0)Ti−1(y0)Ti−1(y0)Ti−1(y0−ax)=Ti+1(y0)2(y0−ax)Ti(y0)Ti
(x)−Ti+1(y0)Ti−1(y0)Ti−1
(x)=Ti+1(y0)2y0Ti(y0)Ti
(x)−Ti+1(y0)Ti−1(y0)Ti−1
(x)−2aTi+1(y0)Ti(y0)xTi
(x)=(1+Ti+1(y0)Ti−1(y0))Ti
(x)−Ti+1(y0)Ti−1(y0)Ti−1
(x)−2aTi+1(y0)Ti(y0)xTi
(x)=αiTi
(x)+(1−αi)Ti−1
(x)−βixTi
(x)
其中系数
α
i
\alpha_i
αi和
β
i
\beta_i
βi的递归计算方式为
β
0
=
2
a
T
0
(
y
0
)
T
1
(
y
0
)
=
2
a
y
0
=
4
λ
ˉ
+
λ
‾
1
β
i
=
1
2
a
T
i
+
1
(
y
0
)
T
i
(
y
0
)
=
λ
ˉ
−
λ
‾
4
2
y
0
T
i
(
y
0
)
−
T
i
−
1
(
y
0
)
T
i
(
y
0
)
=
λ
ˉ
−
λ
‾
2
y
0
−
λ
ˉ
−
λ
‾
4
T
i
−
1
(
y
0
)
T
i
(
y
0
)
=
λ
ˉ
+
λ
‾
2
−
λ
ˉ
−
λ
‾
4
1
2
a
β
i
−
1
=
λ
ˉ
+
λ
‾
2
−
(
λ
ˉ
−
λ
‾
4
)
2
β
i
−
1
α
i
=
2
y
0
T
i
(
y
0
)
T
i
+
1
(
y
0
)
=
2
y
0
2
a
β
i
=
λ
ˉ
+
λ
‾
2
β
i
\beta_0=2a\frac{T_0(y_0)}{T_1(y_0)}=\frac{2a}{y_0}=\frac{4}{\bar{\lambda}+\underline{\lambda}}\\ \begin{aligned} \frac{1}{\beta_i}&=\frac{1}{2a}\frac{T_{i+1}(y_0)}{T_i(y_0)}=\frac{\bar{\lambda}-\underline{\lambda}}{4} \frac{2y_0T_i(y_0)-T_{i-1}(y_0)}{T_i(y_0)}=\frac{\bar{\lambda}-\underline{\lambda}}{2} y_0 - \frac{\bar{\lambda}-\underline{\lambda}}{4} \frac{T_{i-1}(y_0)}{T_i(y_0)}\\ &=\frac{\bar{\lambda}+\underline{\lambda}}{2} - \frac{\bar{\lambda}-\underline{\lambda}}{4} \frac{1}{2a} \beta_{i-1} = \frac{\bar{\lambda}+\underline{\lambda}}{2} - \left(\frac{\bar{\lambda} - \underline{\lambda}}{4} \right)^2 \beta_{i-1} \end{aligned} \\ \alpha_i=2y_0 \frac{T_i(y_0)}{T_{i+1}(y_0)}=\frac{2y_0}{2a}\beta_i=\frac{\bar{\lambda}+\underline{\lambda}}{2} \beta_i
β0=2aT1(y0)T0(y0)=y02a=λˉ+λ4βi1=2a1Ti(y0)Ti+1(y0)=4λˉ−λTi(y0)2y0Ti(y0)−Ti−1(y0)=2λˉ−λy0−4λˉ−λTi(y0)Ti−1(y0)=2λˉ+λ−4λˉ−λ2a1βi−1=2λˉ+λ−(4λˉ−λ)2βi−1αi=2y0Ti+1(y0)Ti(y0)=2a2y0βi=2λˉ+λβi
由此对应的两步Chebyshev迭代格式（具体的导出步骤，笔者也不确定）为：
x
(
i
+
1
)
=
α
i
x
(
i
)
+
(
1
−
α
i
)
x
(
i
−
1
)
+
β
i
r
(
i
)
\bm{x}^{(i+1)}=\alpha_i\bm{x}^{(i)} + (1-\alpha_i)\bm{x}^{(i-1)} + \beta_i\bm{r}^{(i)}
x(i+1)=αix(i)+(1−αi)x(i−1)+βir(i)
共轭梯度（Conjugate Gradient）方法
将原来的代数方程组求解问题
A
x
=
b
\bm{Ax}=\bm{b}
Ax=b转化为全空间无约束的优化问题，目标函数为
φ
=
1
2
(
A
x
,
x
)
−
(
b
,
x
)
\varphi=\frac{1}{2}(\bm{Ax},\bm{x})-(\bm{b}, \bm{x})
φ=21(Ax,x)−(b,x)
它的驻点
x
∗
\bm{x^*}
x∗满足
δ
φ
(
x
∗
)
=
A
x
∗
−
b
=
0
\delta\varphi(\bm{x^*})=\bm{Ax^*}-\bm{b}=0
δφ(x∗)=Ax∗−b=0，即为原方程组的解。目标函数有下界
φ
(
x
∗
)
=
φ
(
A
−
1
b
)
=
−
1
2
(
b
,
A
−
1
b
)
=
−
1
2
(
A
x
∗
,
x
∗
)
φ
(
x
)
−
φ
(
x
∗
)
=
1
2
(
A
(
x
−
x
∗
)
,
x
−
x
∗
)
=
1
2
∣
∣
x
−
x
∗
∣
∣
A
≥
0
\varphi(\bm{x^*})=\varphi(\bm{A}^{-1}\bm{b})=-\frac{1}{2}(\bm{b},\bm{A^{-1}b})=-\frac{1}{2}(\bm{Ax^*},\bm{x^*})\\ \varphi(\bm{x})-\varphi(\bm{x^*})=\frac{1}{2}\left(\bm{A}(\bm{x}-\bm{x^*}), \bm{x}-\bm{x^*}\right)=\frac{1}{2}||\bm{x}-\bm{x^*}||_{\bm{A}} \geq0
φ(x∗)=φ(A−1b)=−21(b,A−1b)=−21(Ax∗,x∗)φ(x)−φ(x∗)=21(A(x−x∗),x−x∗)=21∣∣x−x∗∣∣A≥0
其中由
A
\bm{A}
A内积诱导的
A
\bm{A}
A范数
∣
∣
x
∣
∣
A
2
=
(
A
x
,
x
)
=
(
x
,
x
)
A
||\bm{x}||_{\bm{A}}^2=(\bm{Ax},\bm{x})=(\bm{x},\bm{x})_{\bm{A}}
∣∣x∣∣A2=(Ax,x)=(x,x)A
CG和最速下降都属于子空间搜索方法。如果给定搜索方向
p
(
k
)
\bm{p}^{(k)}
p(k)，对此做一维极小搜索，则可以得到该方向的步长
α
k
\alpha_k
αk：
α
k
=
arg min
α
φ
(
x
(
k
)
+
α
p
(
k
)
)
\alpha_k={\underset {\alpha} {\operatorname {arg\,min} }} \space \varphi({\bm{x}^{(k)}}+\alpha\bm{p}^{(k)})
αk=αargmin φ(x(k)+αp(k))
同之前的迭代法，误差和残量之间关系为 e ( k ) = x ( k ) − x ∗ = − A − 1 r ( k ) \bm{e}^{(k)}=\bm{x}^{(k)}-\bm{x^*}=-\bm{A}^{-1}\bm{r}^{(k)} e(k)=x(k)−x∗=−A−1r(k)。
最速下降法的收敛速度
简单地取
p
(
k
)
=
r
(
k
)
\bm{p}^{(k)}=\bm{r}^{(k)}
p(k)=r(k)，即当前的残量方向。
∣
∣
e
(
k
+
1
)
∣
∣
A
2
=
∣
∣
e
(
k
)
+
α
k
r
(
k
)
∣
∣
A
2
=
(
A
e
(
k
)
+
α
k
A
r
(
k
)
,
e
(
k
)
+
α
k
r
(
k
)
)
=
∣
∣
e
(
k
)
∣
∣
A
2
+
∣
∣
α
k
r
(
k
)
∣
∣
A
2
−
2
(
r
(
k
)
,
r
(
k
)
)
=
∣
∣
e
(
k
)
∣
∣
A
2
−
α
k
(
r
(
k
)
,
r
(
k
)
)
=
∣
∣
e
(
k
)
∣
∣
A
2
−
(
r
(
k
)
,
r
(
k
)
)
(
r
(
k
)
,
r
(
k
)
)
(
r
(
k
)
,
r
(
k
)
)
A
=
∣
∣
e
(
k
)
∣
∣
A
2
(
1
−
1
(
r
(
k
)
,
r
(
k
)
)
A
(
r
(
k
)
,
r
(
k
)
)
(
r
(
k
)
,
r
(
k
)
)
A
−
1
(
r
(
k
)
,
r
(
k
)
)
)
≤
∣
∣
e
(
k
)
∣
∣
A
2
(
1
−
1
(
λ
1
+
λ
n
)
2
4
λ
1
λ
n
)
=
∣
∣
e
(
k
)
∣
∣
A
2
(
λ
1
−
λ
n
)
2
(
λ
1
+
λ
n
)
2
\begin{aligned} ||\bm{e}^{(k+1)}||_{\bm{A}}^2&=||\bm{e}^{(k)}+\alpha_k\bm{r}^{(k)}||_{\bm{A}}^2=\left(\bm{Ae^{(k)}}+\alpha_k\bm{A}\bm{r}^{(k)}, \bm{e}^{(k)}+\alpha_k\bm{r}^{(k)}\right)\\ &=||\bm{e}^{(k)}||_{\bm{A}}^2+||\alpha_k\bm{r}^{(k)}||_{\bm{A}}^2 - 2(\bm{r}^{(k)}, \bm{r}^{(k)})\\ &=||\bm{e}^{(k)}||_{\bm{A}}^2-\alpha_k(\bm{r}^{(k)}, \bm{r}^{(k)})=||\bm{e}^{(k)}||_{\bm{A}}^2-\frac{(\bm{r}^{(k)},\bm{r}^{(k)})(\bm{r}^{(k)},\bm{r}^{(k)})}{{(\bm{r}^{(k)},\bm{r}^{(k)})}_{\bm{A}}}\\ &=||\bm{e}^{(k)}||_{\bm{A}}^2\left( 1 - \frac{1}{ \frac{(\bm{r}^{(k)},\bm{r}^{(k)})_{\bm{A}}}{(\bm{r}^{(k)},\bm{r}^{(k)})} \frac{(\bm{r}^{(k)},\bm{r}^{(k)})_{\bm{A}^{-1}}}{(\bm{r}^{(k)},\bm{r}^{(k)})} } \right)\\ &\leq ||\bm{e}^{(k)}||_{\bm{A}}^2\left( 1 - \frac{1}{ \frac{(\lambda_1+\lambda_n)^2}{4\lambda_1\lambda_n} } \right)=||\bm{e}^{(k)}||_{\bm{A}}^2 \frac{(\lambda_1-\lambda_n)^2}{(\lambda_1+\lambda_n)^2} \end{aligned}
∣∣e(k+1)∣∣A2=∣∣e(k)+αkr(k)∣∣A2=(Ae(k)+αkAr(k),e(k)+αkr(k))=∣∣e(k)∣∣A2+∣∣αkr(k)∣∣A2−2(r(k),r(k))=∣∣e(k)∣∣A2−αk(r(k),r(k))=∣∣e(k)∣∣A2−(r(k),r(k))A(r(k),r(k))(r(k),r(k))=∣∣e(k)∣∣A2⎝⎛1−(r(k),r(k))(r(k),r(k))A(r(k),r(k))(r(k),r(k))A−11⎠⎞≤∣∣e(k)∣∣A2(1−4λ1λn(λ1+λn)21)=∣∣e(k)∣∣A2(λ1+λn)2(λ1−λn)2
所以收敛速度为
∣
∣
e
(
k
)
∣
∣
A
∣
∣
e
(
0
)
∣
∣
A
≤
(
λ
1
−
λ
n
λ
1
+
λ
n
)
k
\frac{||\bm{e}^{(k)}||_{\bm{A}}}{||\bm{e}^{(0)}||_{\bm{A}}} \leq \left(\frac{\lambda_1-\lambda_n}{\lambda_1+\lambda_n}\right)^k
∣∣e(0)∣∣A∣∣e(k)∣∣A≤(λ1+λnλ1−λn)k
CG方法的收敛速度
既然最速下降法是naive地选取局部残量方向作为搜索方向，自然的想法是能否同时优化2个目标（方向 p ( k ) \bm{p}^{(k)} p(k)和步长 α k \alpha_k αk）从而得到一些特殊的性质/好处。CG法希望在确定 p ( k ) \bm{p}^{(k)} p(k)方向之后所进行的一维极小搜索得到的新解 x ( k + 1 ) = x ( k ) + α k p ( k ) \bm{x}^{(k+1)}=\bm{x}^{(k)}+\alpha_k\bm{p}^{(k)} x(k+1)=x(k)+αkp(k)，同时也是加入向量 p ( k ) \bm{p}^{(k)} p(k)后张成的Krylov子空间内的最小值。为达到这个好的性质所做的一些特别的推导在此不提，后续有时间再补上。我们只要知道第 k k k步得到的解 x ( k ) \bm{x}^{(k)} x(k)就是Krylov子空间 K k ( r ( 0 ) ) \bm{K_k}(\bm{r}^{(0)}) Kk(r(0))内的最小值，并根据此得出它的收敛速度。
取 x ( 0 ) = 0 \bm{x}^{(0)}=\bm{0} x(0)=0，则 r ( 0 ) = b \bm{r}^{(0)}=\bm{b} r(0)=b，令 y = A − 1 r ( 0 ) \bm{y}=\bm{A}^{-1}\bm{r}^{(0)} y=A−1r(0)误差之比
∣ ∣ e ( k ) ∣ ∣ A ∣ ∣ e ( 0 ) ∣ ∣ A = min P k − 1 ∈ P k − 1 ∣ ∣ P k − 1 ( A ) r ( 0 ) − A − 1 r ( 0 ) ∣ ∣ A ∣ ∣ A − 1 r ( 0 ) ∣ ∣ A ≤ min P k − 1 ∈ P k − 1 sup y ≠ 0 ∣ ∣ P k − 1 ( A ) A y − y ∣ ∣ A ∣ ∣ y ∣ ∣ A = min P k ∈ P k , P k ( 0 ) = 1 sup y ≠ 0 ∣ ∣ P k ( A ) y ∣ ∣ A ∣ ∣ y ∣ ∣ A \begin{aligned} \frac{||\bm{e}^{(k)}||_{\bm{A}}}{||\bm{e}^{(0)}||_{\bm{A}}}&= {\underset {P_{k-1}\in\mathbf{P_{k-1}}} {\operatorname {min} }} \frac{||P_{k-1}(\bm{A})\bm{r}^{(0)}-\bm{A}^{-1}\bm{r}^{(0)}||_{\bm{A}}}{||\bm{A}^{-1}\bm{r}^{(0)}||_{\bm{A}}}\\ &\leq {\underset {P_{k-1}\in\mathbf{P_{k-1}}} {\operatorname {min} }} {\underset {\bm{y}\neq\bm{0}} {\operatorname {sup} }} \frac{||P_{k-1}(\bm{A})\bm{Ay}-\bm{y}||_{\bm{A}}}{||\bm{y}||_{\bm{A}}}\\ &={\underset {P_{k}\in\mathbf{P_{k}}, P_k(0)=1} {\operatorname {min} }} {\underset {\bm{y}\neq\bm{0}} {\operatorname {sup} }} \frac{||P_{k}(\bm{A})\bm{y}||_{\bm{A}}}{||\bm{y}||_{\bm{A}}}\\ \end{aligned} ∣∣e(0)∣∣A∣∣e(k)∣∣A=Pk−1∈Pk−1min∣∣A−1r(0)∣∣A∣∣Pk−1(A)r(0)−A−1r(0)∣∣A≤Pk−1∈Pk−1miny=0sup∣∣y∣∣A∣∣Pk−1(A)Ay−y∣∣A=Pk∈Pk,Pk(0)=1miny=0sup∣∣y∣∣A∣∣Pk(A)y∣∣A
设 ( λ i , y i ) (\lambda_i, \bm{y}_i) (λi,yi)为 P k ( A ) P_k(\bm{A}) Pk(A)的特征值-向量对，则 y \bm{y} y在这组基底下表示为 y = ∑ i a i y i \bm{y}=\sum_ia_i\bm{y}_i y=∑iaiyi，上式化为
∣ ∣ e ( k ) ∣ ∣ A ∣ ∣ e ( 0 ) ∣ ∣ A ≤ min P k ∈ P k , P k ( 0 ) = 1 sup y ≠ 0 ∣ ∣ ∑ i P k ( λ i ) a i y i ∣ ∣ A ∣ ∣ ∑ i a i y i ∣ ∣ A = min P k ∈ P k , P k ( 0 ) = 1 max λ ∈ σ ( A ) ∣ P k ( λ ) ∣ ≤ min P k ∈ P k , P k ( 0 ) = 1 ∣ ∣ P k ( λ ) ∣ ∣ L ∞ ( λ n λ 1 , 1 ) \begin{aligned} \frac{||\bm{e}^{(k)}||_{\bm{A}}}{||\bm{e}^{(0)}||_{\bm{A}}} &\leq {\underset {P_{k}\in\mathbf{P_{k}}, P_k(0)=1} {\operatorname {min} }} {\underset {\bm{y}\neq\bm{0}} {\operatorname {sup} }} \frac{||\sum_iP_{k}(\lambda_i)a_i\bm{y}_i||_{\bm{A}}} { ||\sum_ia_i\bm{y}_i||_{\bm{A}} } \\ &= {\underset {P_{k}\in\mathbf{P_{k}}, P_k(0)=1} {\operatorname {min} }} {\underset {\lambda\in\sigma(\bm{A})} {\operatorname {max} }} | P_k(\lambda)| \\ &\leq {\underset {P_{k}\in\mathbf{P_{k}}, P_k(0)=1} {\operatorname {min} }} || P_k(\lambda) ||_{L\,\infty\,\left(\frac{\lambda_n}{\lambda_1},1\right)} \\ \end{aligned} ∣∣e(0)∣∣A∣∣e(k)∣∣A≤Pk∈Pk,Pk(0)=1miny=0sup∣∣∑iaiyi∣∣A∣∣∑iPk(λi)aiyi∣∣A=Pk∈Pk,Pk(0)=1minλ∈σ(A)max∣Pk(λ)∣≤Pk∈Pk,Pk(0)=1min∣∣Pk(λ)∣∣L∞(λ1λn,1)
该值即为定义在
[
λ
n
λ
1
,
1
]
[\frac{\lambda_n}{\lambda_1},1]
[λ1λn,1]上且满足
0
0
0点取值为
1
1
1的所有
k
k
k次多项式中，无穷范数最小的那个。这就是Chebyshev多项式，记
α
=
λ
n
λ
1
\alpha=\frac{\lambda_n}{\lambda_1}
α=λ1λn：
T
k
^
(
x
)
=
T
k
~
(
x
)
T
k
~
(
0
)
=
T
k
(
x
−
1
2
(
1
+
α
)
1
2
(
1
−
α
)
)
T
k
(
−
1
+
α
1
−
α
)
=
P
k
∗
(
x
)
\widehat{T_k}(x)=\frac{\widetilde{T_k}(x)}{\widetilde{T_k}(0)}=\frac{T_k\left(\frac{x-\frac{1}{2}(1+\alpha)}{\frac{1}{2}(1-\alpha)} \right)}{T_k\left(-\frac{1+\alpha}{1-\alpha} \right)}=P_k^*(x)
Tk
(x)=Tk
(0)Tk
(x)=Tk(−1−α1+α)Tk(21(1−α)x−21(1+α))=Pk∗(x)
利用其逼近性质，可知收敛速度为
∣
∣
e
(
k
)
∣
∣
A
∣
∣
e
(
0
)
∣
∣
A
≤
∣
∣
P
k
∗
(
x
)
∣
∣
L
∞
(
α
,
1
)
=
1
∣
T
k
(
1
+
α
1
−
α
)
∣
=
2
(
1
−
α
1
+
α
)
k
\frac{||\bm{e}^{(k)}||_{\bm{A}}}{||\bm{e}^{(0)}||_{\bm{A}}} \leq ||P_k^*(x)||_{L\,\infty\,\left(\alpha,1\right)}=\frac{1}{\left| T_k\left( \frac{1+\alpha}{1-\alpha} \right) \right|}=2\left( \frac{1-\sqrt{\alpha}}{1+\sqrt{\alpha}} \right)^k
∣∣e(0)∣∣A∣∣e(k)∣∣A≤∣∣Pk∗(x)∣∣L∞(α,1)=∣∣Tk(1−α1+α)∣∣1=2(1+α
1−α
)k
由于 0 < α < 1 0<\alpha<1 0<α<1，故 α \sqrt{\alpha} α 比 α \alpha α更靠近1，从而CG法有比最速下降更快的收敛速度。
标签：Chebyshev,frac,bm,漫谈,Krylov,Ti,y0,left,lambda 来源： https://blog.csdn.net/weixin_43614211/article/details/122268793
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。