5.5_对分法和反迭代法

5.5 对分法和反迭代法

对分法 (Bisection) 的基本思想是利用惯性定理来计算所需的部分特征值.

定义5.1 设 AA 为对称矩阵, 则其惯性定义为

Inertia(A)=(ν,ζ,π)\operatorname {I n e r t i a} (A) = (\nu , \zeta , \pi)

其中 ν,ζ,π\nu, \zeta, \pi 分别表示 AA 的负特征值,零特征值和正特征值的个数。

定理 5.11 (Sylvester 惯性定理) 设 ARn×nA \in \mathbb{R}^{n \times n} 是对称矩阵, XRn×nX \in \mathbb{R}^{n \times n} 非奇异, 则 XAXX^{\top}AXAA 有相同的惯性.

利用LU分解可得 AzI=LDLTA - zI = LDL^{\mathsf{T}} ,其中 LL 为非奇异下三角矩阵, DD 为对角阵,则

Inertia(AzI)=Inertia(D).\operatorname {I n e r t i a} (A - z I) = \operatorname {I n e r t i a} (D).

由于 DD 是对角矩阵, 所以 Inertia(D)\operatorname{Inertia}(D) 很容易计算.

αR\alpha \in \mathbb{R} , 记 Negcount (A,α)(A, \alpha) 为小于 α\alphaAA 的特征值的个数, 即

Negcount(A,α)=#(λ(A)<α).\operatorname {N e g c o u n t} (A, \alpha) = \# (\lambda (A) < \alpha).

α1<α2\alpha_{1} < \alpha_{2} , 则 AA 在区间 [α1,α2)[\alpha_{1}, \alpha_{2}) 中的特征值个数为

Negcount(A,α2)Negcount(A,α1).\operatorname {N e g c o u n t} (A, \alpha_ {2}) - \operatorname {N e g c o u n t} (A, \alpha_ {1}).

如果 α2α1<tol\alpha_{2} - \alpha_{1} < tol (其中 tol1tol \ll 1 为事先给定的阈值), 且 AA[α1,α2)[\alpha_{1}, \alpha_{2}) 中有特征值, 则我们可将 [α1,α2)[\alpha_{1}, \alpha_{2}) 中的任意一个值作为 AA 在该区间中的特征值的近似. 由此我们可以给出下面的对分法.

算法5.8.对分法:计算 AA[a,b)[a,b) 中的所有特征值

1: Let toltol be a given threshold
2: compute na=Negcount(A,a)n_a = \mathrm{Negcount}(A, a)
3: compute nb=Negcount(A,b)n_b = \mathrm{Negcount}(A, b)
4: if na=nbn_a = n_b then return %此时[a,b)中没有A的特征值
5: put (a,na,b,nb)(a, n_a, b, n_b) onto worklist
6: % worklist 中的元素是 “四元素对”, 即由四个数组成的数对
7: while worklist not empty do
8: remove (low,nlow,up,nup)(low, n_{low}, up, n_{up}) from the worklist
9: %\% (low, nlow,up,nupn_{low}, up, n_{up} ) 是 worklist 中的任意一个元素
10: if (uplow)<tol(up - low) < tol then
11: print "There are nupnlown_{up} - n_{low} eigenvalues in [low, up)"
12: else
13: compute mid=(low+up)/2\text{mid} = (\text{low} + \text{up}) / 2
14: compute nmid=Negcount(A,mid)n_{mid} = \mathrm{Negcount}(A, mid)
15: if (nmid>nlow)(n_{mid} > n_{low}) then

16: put (low, nlown_{low} , mid, nmidn_{mid} ) onto worklist
17: end if
18: if (nup>nmid)(n_{up} > n_{mid}) then
19: put (mid, nmidn_{mid} , up, nupn_{up} ) onto worklist
20: end if
21: end if
22: end while

显然,对分法的主要运算量集中在计算 Negcount(A,z)\mathrm{Negcount}(A,z) 。通常是事先将 AA 转化成对称三对角矩阵,这样计算 AzIA - zILDLT\mathrm{LDL}^{\mathrm{T}} 分解就非常简单:

AzI=[a1zb1b1bn1bn1anz]=[1l1ln11][d1dn][1l1ln11]LDLT.\begin{array}{l} A - z I = \left[ \begin{array}{c c c c} a _ {1} - z & b _ {1} & & \\ b _ {1} & \ddots & \ddots & \\ & \ddots & \ddots & b _ {n - 1} \\ & & b _ {n - 1} & a _ {n} - z \end{array} \right] \\ = \left[ \begin{array}{c c c c} 1 & & & \\ l _ {1} & \ddots & & \\ & \ddots & \ddots & \\ & & l _ {n - 1} & 1 \end{array} \right] \left[ \begin{array}{c c c c} d _ {1} & & & \\ & \ddots & & \\ & & \ddots & \\ & & & d _ {n} \end{array} \right] \left[ \begin{array}{c c c c} 1 & l _ {1} & & \\ & \ddots & \ddots & \\ & & \ddots & l _ {n - 1} \\ & & & 1 \end{array} \right] \triangleq L D L ^ {\mathsf {T}}. \\ \end{array}

利用待定系数法, 可以得到下面的递推公式

d1=a1z,di=(aiz)bi12di1,i=2,3,,n.(5.8)d _ {1} = a _ {1} - z, \quad d _ {i} = \left(a _ {i} - z\right) - \frac {b _ {i - 1} ^ {2}}{d _ {i - 1}}, \quad i = 2, 3, \dots , n. \tag {5.8}

用上面的公式计算 did_{i} 的运算量约为 4n4n

注意这里没有选主元, 但针对对称三对角矩阵, 该算法是非常稳定的, 即使当 did_{i} 有可能很小时, 算法依然很稳定.

定理5.12 [30]利用公式(5.5)计算所得的 did_{i} 与精确计算 A^\hat{A}d^i\hat{d}_i 有相同的符号,故有相同的惯性.这里 A^\hat{A}AA 非常接近,即

A^(i,i)=ai,A^(i,i+1)=bi(1+εi),\hat {A} (i, i) = a _ {i}, \quad \hat {A} (i, i + 1) = b _ {i} (1 + \varepsilon_ {i}),

其中 εi2.5ε+O(ε2)|\varepsilon_i| \leq 2.5\varepsilon + \mathcal{O}(\varepsilon^2) , 这里的 ε\varepsilon 为机器精度.

由于单独调用一次 Negcount 的运算量为 4n4n , 故计算 kk 个特征值的总运算量约为 O(kn)\mathcal{O}(kn) .

当特征值计算出来后, 我们可以使用带位移的逆迭代来计算对应的特征向量. 通常只需迭代 1 至 2 次即可, 由于 AA 是三对角矩阵, 故计算每个特征向量的运算量为 O(n)\mathcal{O}(n) . 整个合起来就构成对分法和逆迭代.

当特征值紧靠在一起时,计算出来的特征向量可能会失去正交性,此时需要进行再正交化,可通过MGS的QR分解来实现.