Arnoldi迭代

目标：得到一组标准正交基${\mathbf{v}}_1,\cdots ,{\mathbf{v}}_m$，使得
K m ( A , r 0 ) = span ⁡ { v 1 , ⋯   , v m } \mathcal{K}_m\left(\mathbf{A},\mathbf{r}_0\right)=\operatorname{span}\left\{\mathbf{v}_1,\cdots,\mathbf{v}_m\right\} Km​(A,r0​)=span{v1​,⋯,vm​}
显然${\mathbf{v}}_1=\frac{{\mathbf{r}}_0}{\Vert {\mathbf{r}}_0\Vert }$

假设已经得到${\mathcal{K}}_j\left(\mathbf{A},{\mathbf{r}}_0\right)=\mathrm{span}\left({\mathbf{v}}_1,\cdots ,{\mathbf{v}}_j\right)$
现在要求${\mathbf{v}}_{j+1}$
$${\mathbf{v}}_j\in {\mathcal{K}}_j\left(\mathbf{A},{\mathbf{r}}_0\right)\Rightarrow \mathbf{A}{\mathbf{v}}_j\in \mathbf{A}{\mathcal{K}}_j\left(\mathbf{A},{\mathbf{r}}_0\right)\subset {\mathcal{K}}_{j+1}\left(\mathbf{A},{\mathbf{r}}_0\right)$$
那么
$$\mathbf{A}{\mathbf{v}}_j=\sum _{i=0}^{j+1}{h}_{ij}{\mathbf{v}}_i$$

$${\mathbf{v}}_i^T\mathbf{A}{\mathbf{v}}_j={\mathbf{h}}_{ij}\left(i=0,1,\cdots ,j\right)$$
于是
$$h_{j+1,j}{\mathbf{v}}_{j+1}=\mathbf{A}{\mathbf{v}}_j-\sum _{i=0}^j{h}_{ij}{\mathbf{v}}_j$$
所以
$$\{\begin{array}{l}{\mathbf{v}}_{j+1}=\frac{\mathbf{A}{\mathbf{v}}_j-{\sum }_{i=0}^j{h}_{ij}{\mathbf{v}}_j}{\Vert \mathbf{A}{\mathbf{v}}_j-{\sum }_{i=0}^j{h}_{ij}{\mathbf{v}}_j\Vert }\\ h_{j+1,j}=\Vert \mathbf{A}{\mathbf{v}}_j-{\sum }_{i=0}^j{h}_{ij}{\mathbf{v}}_j\Vert \end{array}$$

如果某一轮迭代中，$h_{j+1,j}=0$,则可以提前停止
可以证明Arnoldi迭代提前停止，当且仅当$\dim {\mathcal{K}}_k<k$

令$q_m\left(\mathbf{A}\right)$是$m$阶多项式
${\mathbf{v}}_1=q_0\left(\mathbf{A}\right){\mathbf{v}}_1$
假设当$i\le j$时，${\mathbf{v}}_j=q_{j-1}\left(\mathbf{A}\right){\mathbf{v}}_1$成立，则
$h_{j+1,j}{\mathbf{v}}_{j+1}=\mathbf{A}{\mathbf{v}}_j-{\sum }_{i=0}^j{h}_{ij}{\mathbf{v}}_j=q_{j+1}\left(\mathbf{A}\right){\mathbf{v}}_1$
也就是说，只要没有提前停止，就有 K m ( A , r 0 ) = span ⁡ { v 1 , ⋯   , v m } \mathcal{K}_m\left(\mathbf{A},\mathbf{r}_0\right)=\operatorname{span}\left\{\mathbf{v}_1,\cdots,\mathbf{v}_m\right\} Km​(A,r0​)=span{v1​,⋯,vm​}

Arnoldi迭代最后会产生一个Hessenberg矩阵$\bar{\mathbf{H}}\in {\mathbb{R}}^{\left(k+1\right)\times k}$
使得
$$\begin{array}{rl}\mathbf{A}{\mathbf{V}}_m& ={\mathbf{V}}_{m+1}{\bar{\mathbf{H}}}_m\\ & ={\mathbf{V}}_m{\mathbf{H}}_m+{\mathbf{v}}_{m+1}\left(h_{m+1,m}{\mathbf{e}}_m^T\right)\\ & ={\mathbf{V}}_m{\mathbf{H}}_m+{\mathbf{w}}_m{\mathbf{e}}_m^T\\ {\mathbf{V}}_m^T\mathbf{A}{\mathbf{V}}_m& ={\mathbf{H}}_m\end{array}$$
其中${\mathbf{H}}_m$是${\bar{\mathbf{H}}}_m$的前$m$行
注意${\mathbf{V}}_{k+1}^T{\mathbf{V}}_{k+1}=\mathbf{I}$,但是${\mathbf{V}}_{k+1}{\mathbf{V}}_{k+1}^T$不一定等于单位矩阵

修改版

主要区别就是$h_{ij}=\left(w_j,v_i\right)$，因为${\mathbf{v}}_i^T{\mathbf{v}}_j=0$,所以数值上是一样的

Householder变换

其中3-5行是产生Householder的反射矩阵

其实相当于对$\mathbf{v},\mathbf{A}{\mathbf{v}}_1,\mathbf{A}{\mathbf{v}}_2,\cdots ,\mathbf{A}{\mathbf{v}}_m$做householder正交化或者householder版QR分解

令
$${\mathbf{Q}}_j={\mathbf{P}}_j{\mathbf{P}}_{j-1}\cdots {\mathbf{P}}_1$$
由第8行
$${\mathbf{Q}}_j\mathbf{A}{\mathbf{v}}_j={\mathbf{z}}_{j+1}$$
由第6行
$${\mathbf{h}}_j={\mathbf{P}}_{j+1}{\mathbf{z}}_{j+1}={\mathbf{P}}_{j+1}{\mathbf{Q}}_j\mathbf{A}{\mathbf{v}}_j={\mathbf{Q}}_{j+1}\mathbf{A}{\mathbf{v}}_j$$
注意到${\mathbf{h}}_j$的$j+1,\cdots ,n$分量都是$0$，所以${\mathbf{P}}_i{\mathbf{h}}_j={\mathbf{h}}_j\left(i\ge j+2\right)$,于是
$${\mathbf{h}}_j={\mathbf{P}}_m\cdots {\mathbf{P}}_{j+2}{\mathbf{h}}_j={\mathbf{Q}}_m\mathbf{A}{\mathbf{v}}_j$$
于是
$${\mathbf{Q}}_m\left(\mathbf{v},\mathbf{A}{\mathbf{v}}_1,\cdots ,\mathbf{A}{\mathbf{v}}_m\right)=\left({\mathbf{h}}_0,\cdots ,{\mathbf{h}}_m\right)$$
$\left({\mathbf{h}}_0,\cdots ,{\mathbf{h}}_m\right)$是$n\times \left(m+1\right)$的矩阵，${\mathbf{Q}}_m$是标准正交矩阵
令${\bar{\mathbf{H}}}_m$为$\left({\mathbf{h}}_1,\cdots ,{\mathbf{h}}_m\right)$的前$m+1$行
$$\mathbf{A}{\mathbf{v}}_j={\mathbf{Q}}_{j+1}^T{\mathbf{h}}_j={\mathbf{Q}}_{j+1}^T\sum _{i=1}^{j+1}{h}_{ij}{\mathbf{e}}_i=\sum _{i=1}^{j+1}{h}_{ij}{\mathbf{Q}}_{j+1}^T{\mathbf{e}}_i$$
注意到${\mathbf{P}}_k{\mathbf{e}}_i={\mathbf{e}}_i\left(i<k\right)$
$${\mathbf{Q}}_{j+1}^T{\mathbf{e}}_i={\mathbf{P}}_1\cdots {\mathbf{P}}_{j+1}{\mathbf{e}}_i={\mathbf{v}}_i\left(i\le j+1\right)$$
于是
$$\mathbf{A}{\mathbf{v}}_j=\sum _{i=1}^{j+1}{h}_{ij}{\mathbf{v}}_i$$
所以
$$\mathbf{A}{\mathbf{V}}_m={\mathbf{V}}_{m+1}{\bar{\mathbf{H}}}_m$$

广义极小残量法

广义极小残量法（Generalized Minimal RESidual，GMRES）
考虑大型线性方程组
$$\mathbf{Ax}=\mathbf{b}$$

考虑Krylov子空间
$${\mathcal{K}}_m\left(\mathbf{A},{\mathbf{r}}_0\right)=\mathrm{span}\left\{{\mathbf{r}}_0,\mathbf{A}{\mathbf{r}}_0,\cdots ,{\mathbf{A}}^{m-1}{\mathbf{r}}_0\right\}$$
设${\mathbf{x}}^{(0)}$为起点
GMRES考虑解$\underset{\mathbf{x}\in {\mathbf{x}}^{(0)}+{\mathcal{K}}_m\left(\mathbf{A},{\mathbf{r}}_0\right)}{\min }\Vert \mathbf{Ax}-\mathbf{b}\Vert$

设${\mathbf{v}}_1,\cdots ,{\mathbf{v}}_m$为一组标准正交基
${\mathbf{V}}_m=\left({\mathbf{v}}_1,\cdots ,{\mathbf{v}}_m\right)\in {\mathbb{R}}^{n\times m}$
则${\mathbf{x}}^{(m)}\in {\mathbf{x}}^{(0)}+{\mathcal{K}}_m\left(\mathbf{A},{\mathbf{r}}_0\right)$可以写作${\mathbf{x}}^{(m)}={\mathbf{x}}^{(0)}-{\mathbf{V}}_m\mathbf{y}$,其中$\mathbf{y}\in {\mathbb{R}}^m$
由Arnoldi迭代
$$\begin{array}{rl}\mathbf{b}-\mathbf{Ax}& =\mathbf{b}-\mathbf{A}\left({\mathbf{x}}^{(0)}+{\mathbf{V}}_m\mathbf{y}\right)\\ & ={\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}\\ & =\beta {\mathbf{v}}_1-{\mathbf{V}}_{m+1}{\bar{\mathbf{H}}}_m\mathbf{y}\\ & ={\mathbf{V}}_{m+1}\left(\beta {\mathbf{e}}_1-{\bar{\mathbf{H}}}_m\mathbf{y}\right)\end{array}$$
其中${\mathbf{r}}_0=\mathbf{b}-\mathbf{A}{\mathbf{x}}^{(0)},\beta =\Vert {\mathbf{r}}_0\Vert ,{\mathbf{v}}_1=\frac{{\mathbf{r}}_0}{\beta }$

于是
$$\underset{\mathbf{x}\in {\mathbf{x}}^{(0)}+{\mathcal{K}}_m\left(\mathbf{A},{\mathbf{r}}_0\right)}{\min }\Vert \mathbf{Ax}-\mathbf{b}\Vert =\underset{\mathbf{y}\in {\mathbb{R}}^m}{\min }\Vert {\mathbf{r}}_0-\mathbf{A}{\mathbf{V}}_m\mathbf{y}\Vert =\underset{\mathbf{y}\in {\mathbb{R}}^m}{\min }\Vert \beta {\mathbf{e}}_1-{\bar{\mathbf{H}}}_m\mathbf{y}\Vert$$

Arnoldi-Householder

设${\mathbf{y}}_m=\left(\begin{array}{c}{\eta }_1\\ ⋮\\ {\eta }_m\end{array}\right)$
则${\mathbf{x}}^{(m)}={\mathbf{x}}^{(0)}+{\eta }_1{\mathbf{v}}_1+\cdots +{\eta }_m{\mathbf{v}}_m$
因为
$${\mathbf{v}}_j={\mathbf{P}}_1\cdots {\mathbf{P}}_j{\mathbf{e}}_j$$
所以
$$\begin{array}{rl}{\mathbf{x}}^{(m)}& ={\mathbf{x}}^{(0)}+{\eta }_1{\mathbf{v}}_1+\cdots +{\eta }_m{\mathbf{v}}_m\\ & ={\mathbf{x}}^{(0)}+{\eta }_1{\mathbf{P}}_1{\mathbf{e}}_1+\cdots +{\eta }_m{\mathbf{P}}_1\cdots {\mathbf{P}}_m{\mathbf{e}}_m\\ & ={\mathbf{x}}^{(0)}+{\mathbf{P}}_1\left({\eta }_1{\mathbf{e}}_1+{\mathbf{P}}_2\left({\eta }_2{\mathbf{e}}_2+\cdots {\mathbf{P}}_{m-1}\left({\eta }_{m-1}{\mathbf{e}}_{m-1}+{\mathbf{P}}_m{\eta }_m{\mathbf{e}}_m\right)\right)\right)\end{array}$$
最后就可以减小代价

进一步化简

使用Givens旋转变换，把Hessenberg矩阵${\bar{\mathbf{H}}}_m$化成上三角矩阵加一行$0$
设${\Omega }_i$为第$i$个旋转矩阵
$$\begin{gathered} {\mathbf{Q}}_m={\Omega }_m{\Omega }_{m-1}\cdots {\Omega }_1 \\ {\bar{\mathbf{R}}}_m={\mathbf{Q}}_m{\bar{\mathbf{H}}}_m \\ {\bar{\mathbf{g}}}_m={\mathbf{Q}}_m\left(\beta {\mathbf{e}}_1\right)={\left({\gamma }_1,\cdots ,{\gamma }_{m+1}\right)}^T \end{gathered}$$
于是
$$\underset{\mathbf{y}\in {\mathbb{R}}^m}{\min }\Vert \beta {\mathbf{e}}_1-{\bar{\mathbf{H}}}_m\mathbf{y}\Vert =\underset{\mathbf{y}\in {\mathbb{R}}^m}{\min }\Vert {\bar{\mathbf{g}}}_m-{\bar{\mathbf{R}}}_m\mathbf{y}\Vert$$
令${\mathbf{R}}_m$为${\bar{\mathbf{R}}}_m$删掉最后一行，${\mathbf{g}}_m=\left(\begin{array}{c}{\gamma }_1\\ ⋮\\ {\gamma }_m\end{array}\right)$
$$\begin{array}{rl}\Vert {\bar{\mathbf{g}}}_m-{\bar{\mathbf{R}}}_m\mathbf{y}{\Vert }^2& =\Vert {\mathbf{g}}_m-{\mathbf{R}}_m\mathbf{y}{\Vert }^2+{\left|{\gamma }_{m+1}\right|}^2\end{array}$$
$$\mathbf{A}{\mathbf{V}}_m={\mathbf{V}}_{m+1}{\bar{\mathbf{H}}}_m={\mathbf{V}}_{m+1}{\mathbf{Q}}_m^T{\bar{\mathbf{R}}}_m\Rightarrow \mathrm{rank}\left(\mathbf{A}{\mathbf{V}}_m\right)=\mathrm{rank}\left({\bar{\mathbf{R}}}_m\right)$$
所以只要$r_{ii}\ne 0$,$\mathbf{A}{\mathbf{V}}_m$非奇异
因为${\mathbf{R}}_m$是上三角矩阵，可以通过回代解，${\mathbf{y}}_m={\mathbf{R}}_m^{-1}{\mathbf{g}}_m$
$$\underset{\mathbf{y}\in {\mathbb{R}}^m}{\min }\Vert {\bar{\mathbf{g}}}_m-{\bar{\mathbf{R}}}_m\mathbf{y}\Vert =\left|{\gamma }_{m+1}\right|$$

重启

收敛性

假设$\mathbf{A}$可对角化，即$\mathbf{A}=\mathbf{X}\Lambda {\mathbf{X}}^{-1}$
其中$\Lambda =\mathrm{diag}\left({\lambda }_1,\cdots ,{\lambda }_n\right)$是$\mathbf{A}$的特征值组成的对角矩阵
令
$${ϵ}^{(m)}=\underset{p\in {\mathbb{P}}_m,p\left(0\right)=1}{\min }\underset{i=1,\cdots ,n}{\max }\left|p\left({\lambda }_i\right)\right|$$
则
$$\Vert {\mathbf{r}}_m\Vert \le {\kappa }_2\left(\mathbf{x}\right){ϵ}^{(m)}\Vert {\mathbf{r}}_0\Vert$$
证明：因为$\mathbf{x}\in {\mathcal{K}}_m$,有$\mathbf{b}-\mathbf{Ax}=p\left(\mathbf{A}\right){\mathbf{r}}_0$
$$\Vert \mathbf{b}-\mathbf{Ax}\Vert =\Vert \mathbf{X}p\left(\Lambda \right){\mathbf{X}}^{-1}{\mathbf{r}}_0\Vert \le \Vert \mathbf{X}{\Vert }_2\Vert {\mathbf{X}}^{-1}{\Vert }_2\Vert {\mathbf{r}}_0\Vert \Vert p\left(\Lambda \right)\Vert \le {\kappa }_2\left(\mathbf{X}\right){ϵ}^{(m)}\Vert {\mathbf{r}}_0\Vert$$

推论

假设$\mathbf{A}$可对角化，即$\mathbf{A}=\mathbf{X}\Lambda {\mathbf{X}}^{-1}$
其中$\Lambda =\mathrm{diag}\left({\lambda }_1,\cdots ,{\lambda }_n\right)$是$\mathbf{A}$的特征值组成的对角矩阵
假设$\mathbf{A}$的特征值落在椭圆$E\left(c,d,a\right)$中（$c$为中心，$d$为焦距长，$a$为长半轴）
则
$$\Vert {\mathbf{r}}_m\Vert \le {\kappa }_2\left(\mathbf{X}\right)\frac{C_m\left(\frac{a}{d}\right)}{\left|C_m\left(\frac{c}{d}\right)\right|}\Vert {\mathbf{r}}_0\Vert$$
其中$C$为切比雪夫多项式

证明：
摸了

参考

Iterative methods for sparse linear systems