## Krylov 子空間法──線性方程的數值解法 (一)：Arnoldi 與 Lanczos 算法

$A$ 為一 $n\times n$ 階複矩陣，$\mathbf{b}\in\mathbb{C}^n$ 為一非零向量。向量序列 $\{\mathbf{b},A\mathbf{b},\ldots,A^{j-1}\mathbf{b}\}$ 稱為 Krylov 序列，此序列所生成的子空間稱為 Krylov 子空間 (見“Krylov 子空間法”)，記為

$\displaystyle \mathcal{K}_j(A,\mathbf{b})=\text{span}\{\mathbf{b},A\mathbf{b},\ldots,A^{j-1}\mathbf{b}\}$

$A^k\mathbf{b}=c_0\mathbf{b}+c_1A\mathbf{b}+\cdots+c_{k-1}A^{k-1}\mathbf{b}$

$\displaystyle m(t)=t^k-c_{k-1}t^{k-1}-\cdots-c_1t-c_0$

$\displaystyle A\left(\frac{A^{k-1}\mathbf{b}-c_{k-1}A^{k-2}\mathbf{b}-\cdots-c_1\mathbf{b}}{c_0}\right)=\mathbf{b}$

$\displaystyle \mathbf{x}_\ast=\frac{1}{c_0}\left(A^{k-1}\mathbf{b}-c_{k-1}A^{k-2}\mathbf{b}-\cdots-c_1\mathbf{b}\right)$

$\displaystyle \mathbf{y}_j=A^j\mathbf{v}-\sum_{i=1}^j(\mathbf{q}_i^\ast A^j\mathbf{v})\mathbf{q}_i$

\displaystyle \begin{aligned} \mathcal{K}_{j+1}(A,\mathbf{v})&=\hbox{span}\{\mathbf{v},A\mathbf{v},A^2\mathbf{v},\ldots,A^j\mathbf{v}\}\\ &=\hbox{span}\{\mathbf{q}_1,A\mathbf{q}_1,A^2\mathbf{q}_1,\ldots,A^j\mathbf{q}_1\}\\ &=\hbox{span}\{\mathbf{q}_1,d_1\mathbf{q}_1+d_2\mathbf{q}_2,A(d_1\mathbf{q}_1+d_2\mathbf{q}_2),\ldots,A^{j-1}(d_1\mathbf{q}_1+d_2\mathbf{q}_2)\}\\ &=\hbox{span}\{\mathbf{q}_1,\mathbf{q}_2,A\mathbf{q}_2,\ldots,A^{j-1}\mathbf{q}_2\}\\ &=\cdots\\ &=\hbox{span}\{\mathbf{q}_1,\mathbf{q}_2,\ldots,\mathbf{q}_j,A\mathbf{q}_j\}. \end{aligned}

$\displaystyle \mathbf{u}_j=A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i$

$\mathbf{u}_j=\mathbf{0}$，停止計算，表示 $\mathcal{K}_{j+1}(A,\mathbf{v})=\mathcal{K}_j(A,\mathbf{v})$；否則令 $\mathbf{q}_{j+1}=\mathbf{u}_j/\Vert\mathbf{u}_j\Vert$，並重複上述步驟繼續擴充 Krylov 子空間。因為 $\mathbf{q}_{j+1}\perp\hbox{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\}$

$\displaystyle \mathbf{q}_{j+1}^\ast\mathbf{u}_j=\mathbf{q}_{j+1}^\ast\left(A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\right)=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j$

\displaystyle\begin{aligned} A\mathbf{q}_j&=\mathbf{u}_j+\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\ &=\Vert\mathbf{u}_j\Vert\mathbf{q}_{j+1}+\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\ &=\sum_{i=1}^{j+1}(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i\\ &=\sum_{i=1}^{j+1}h_{ij}\mathbf{q}_i, \end{aligned}

Arnoldi-MGS 算法

1. $\mathbf{q}_1\leftarrow\mathbf{v}/\Vert\mathbf{v}\Vert$
2. 對於 $j=1,2,\ldots,n$
3.     $\mathbf{u}\leftarrow A\mathbf{q}_j$
4.     對於 $i=1,\ldots,j$
5.         $h_{ij}\leftarrow\mathbf{q}_i^\ast\mathbf{u}$
6.         $\mathbf{u}\leftarrow\mathbf{u}-h_{ij}\mathbf{q}_i$
7.     $h_{j+1,j}\leftarrow\Vert\mathbf{u}\Vert$
8.     若 $h_{j+1,j}=0$，停止計算。
9.     $\mathbf{q}_{j+1}\leftarrow\mathbf{u}/h_{j+1,j}$

Arnoldi-MGS 算法可用簡明的矩陣運算表示。令 $n\times j$ 階矩陣 $Q_j=\begin{bmatrix} \mathbf{q}_1&\cdots&\mathbf{q}_j \end{bmatrix}$ 以及 $j\times j$ 階上 Hessenberg 矩陣 $H_j=[h_{il}]$，即 $h_{il}=0$$i>l+1$ (見“特殊矩陣 (19)：Hessenberg 矩陣”)。寫出單範正交基底的關鍵式

$\displaystyle A\mathbf{q}_l=\begin{bmatrix} \mathbf{q}_1&\cdots&\mathbf{q}_l \end{bmatrix}\begin{bmatrix} h_{1l}\\ \vdots\\ h_{ll} \end{bmatrix}+\mathbf{q}_{l+1}h_{l+1,l},~~l=1,\ldots,j$

\displaystyle\begin{aligned} AQ_j&=\begin{bmatrix} A\mathbf{q}_1&A\mathbf{q}_2&\cdots &A\mathbf{q}_j \end{bmatrix}\\ &=\begin{bmatrix} \mathbf{q}_1&\mathbf{q}_2&\cdots&\mathbf{q}_j \end{bmatrix}\begin{bmatrix} h_{11}&h_{12}& h_{13}&\cdots&h_{1j}\\ h_{21}&h_{22}&h_{23}&\cdots&h_{2j}\\ 0&h_{32}&h_{33}&\cdots&h_{3j}\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\cdots&0&h_{j,j-1}&h_{jj} \end{bmatrix}+\begin{bmatrix} \mathbf{0}&\cdots&\mathbf{0}&\mathbf{q}_{j+1}h_{j+1,j} \end{bmatrix}\\ &=Q_jH_j+\mathbf{q}_{j+1}h_{j+1,j}\mathbf{e}_j^T, \end{aligned}

$AQ_j=Q_{j+1}\tilde{H_j}$

$\displaystyle \tilde{H}_j=\begin{bmatrix} H_j\\ h_{j+1,j}\mathbf{e}_j^T \end{bmatrix}=\begin{bmatrix} h_{11}&h_{12}& h_{13}&\cdots&h_{1j}\\ h_{21}&h_{22}&h_{23}&\cdots&h_{2j}\\ 0&h_{32}&h_{33}&\cdots&h_{3j}\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\cdots&0&h_{j,j-1}&h_{jj}\\ 0&\cdots&\cdots&0&h_{j+1,j} \end{bmatrix}$

$Q_j^\ast AQ_j=Q_j^\ast Q_{j+1}\tilde{H_j}=H_j$

$AQ_k=Q_kH_k$

$A$ 是 Hermitian 矩陣，即 $A^\ast=A$，則

$H_j^\ast=(Q_j^\ast AQ_j)^\ast=Q_j^\ast AQ_j=H_j$

$\displaystyle T_j=\begin{bmatrix} \alpha_1&\beta_1& & & \\ \beta_1&\alpha_2&\beta_2&&\\ & \beta_2&\alpha_3&\ddots&\\ && \ddots&\ddots&\beta_{j-1}\\ &&&\beta_{j-1}&\alpha_j \end{bmatrix}$

$\displaystyle \mathbf{u}_j=A\mathbf{q}_j-\sum_{i=1}^j(\mathbf{q}_i^\ast A\mathbf{q}_j)\mathbf{q}_i=A\mathbf{q}_j-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}$

\begin{aligned} \Vert\mathbf{u}_j\Vert&=\mathbf{q}_{j+1}^\ast\mathbf{u}_j=\mathbf{q}_{j+1}^\ast\left(A\mathbf{q}_j-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}\right)\\ &=\mathbf{q}_{j+1}^\ast A\mathbf{q}_j=(\mathbf{q}_j^\ast A\mathbf{q}_{j+1})^\ast=\bar\beta_j,\end{aligned}

$A\mathbf{q}_j=\beta_{j-1}\mathbf{q}_{j-1}+\alpha_j\mathbf{q}_j+\beta_j\mathbf{q}_{j+1}$

1. $\mathbf{q}_1\leftarrow\mathbf{v}/\Vert\mathbf{v}\Vert$$\mathbf{q}_0\leftarrow\mathbf{0}$$\beta_0\leftarrow 0$
2. 對於 $j=1,2,\ldots,n$
3.     $\mathbf{u}\leftarrow A\mathbf{q}_j$
4.     $\alpha_j\leftarrow \mathbf{q}_j^\ast\mathbf{u}$
5.     $\mathbf{u}\leftarrow\mathbf{u}-\alpha_j\mathbf{q}_j-\beta_{j-1}\mathbf{q}_{j-1}$
6.     $\beta_j\leftarrow\Vert\mathbf{u}\Vert$
7.     若 $\beta_j=0$，停止計算。
8.     $\mathbf{q}_{j+1}\leftarrow\mathbf{u}/\beta_j$

$AQ_j=Q_jT_j+\mathbf{q}_{j+1}\beta_j\mathbf{e}_j^T$

$AQ_j=Q_{j+1}\tilde{T}_j$

$AQ_k=Q_kT_k$

$A$ 為大尺寸矩陣，我們可以採用迭代法求解。定常迭代法是較古老且相對簡單的線性方程迭代算法，但僅保證有限的 $A$ 矩陣具備收斂性，如對角佔優矩陣與正定矩陣 (見“定常迭代法──線性方程迭代解法”)。即便 Krylov 序列未直接提供線性方程 $A\mathbf{x}=\mathbf{b}$ 的精確解 $\mathbf{x}_\ast=A^{-1}\mathbf{b}$，但至少我們知道精確解 $\mathbf{x}_\ast$ 屬於 Krylov 子空間 $\mathcal{K}_k(A,\mathbf{b})$，其中 $k$$\mathbf{b}$ 相對 $A$ 的最小多項式的次數。下文將介紹通過 Krylov 序列的線性方程的迭代解法，統稱為 Krylov 子空間法，適用於任何型態的可逆矩陣，是目前數值線性代數中最成功的一類算法。

This entry was posted in 線性代數專欄, 數值線性代數 and tagged , , , , , , , , . Bookmark the permalink.

### 3 Responses to Krylov 子空間法──線性方程的數值解法 (一)：Arnoldi 與 Lanczos 算法

1. Ran says:

请问老师 subspace iteration 和 Lanczos method 有什么异同呢？

• ccjou says:

我還沒寫完，後續會介紹。簡單講，Krylov subspace methods由兩個子空間($U_j$$V_j$)定義：
1. $x_j-x_0\in U_j$,
2. $b-Ax_j\perp V_j$.

兩類方法的最大差異在於$V_j$

2. bingkunblog says:

周老师您好~能否讲下LSQR 算法？