线性方程的迭代方法
本文将介绍求解线性方程· \(Ax = b\) 的几种方法。 其中 A 是一个大型矩阵，通常以算子的形式给出，例如在偏微分方程中。这类问题的规模太大，以至于像LU分解之类的直接方法受内存限制无法使用，故需要使用迭代求解。
静态（Stationary）方法
采用如下迭代格式：
\[x_{k+1}=-M^{-1}(A-M)x_k+M^{-1}b\tag{1} \]M是一个相对于A更便于求逆的矩阵。可以很方便的验证上述格式在\(Ax=b\) 是不动点。M有3种常用的选择
- \(A\) 的对角元 \(D\) (Jacobi 格式)，其具体迭代为：
- \(A\) 的下三角元 \(D+L\) (Gauss-Siedel 格式)，其具体迭代为：
- 以上二者线性组合 \({1\over \omega}D+L\) (SOR 格式)，其具体迭代为：
共轭梯度法(Conjugate gradient, CG)
当\(A\) 对称正定时，可将原问题化为最小化 \(\phi(x)={1\over 2}x^TAx-b^Tx\). 从一个点 \(x\)开始，沿一个方向 \(p\) 搜索 \(\phi(x+ap)\) 的最小值，得到
\[a=\arg\min\phi(x+ap)={p^Tr \over p^TAp }, \quad\text{where }r=b-Ax\tag{5} \]再把 \(x+ap\) 当做新的起点，换个方向重复以上步骤就可以不断减小 \(r\) 直到找到解。剩下就是选择搜索的方向。注意到\(\nabla\phi(x)=-r\), 即 \(r\) 是 \(\phi(x)\) 下降最快的方向，故可以选取 \(r\) 作为每次搜索的方向。但实际上有更好的办法。记 \(P_{i-1}=[p_0, \cdots, p_{i-1}]\in\mathbb{R}^{n\cross i}\) 为前 \(i\) 次的搜索方向，则 \(x_i=x_0+P_{i-1}y_{i-1}\), 其中 \(y_{i-1}\) 为每次迭代的 \(a\) 。则在第 \(i\) 次迭代中
\[\begin{align} \phi(x_i+ap_i)&=\phi(x_0+P_{i-1}y_{i-1}+ap_i)\\ &=\phi(x_0+P_{i-1}y_{i-1})+ap_i^TA(x_0+P_{i-1}y_{i-1})+{1\over2}a^2p_i^TAp_i-ab^Tp_i\\ &=\phi(x_0+P_{i-1}y_{i-1})-ap_i^Tr_0+ay_{i-1}^TP_{i-1}^TAp_i+{a^2\over2}p_i^TAp_i \tag{6} \end{align} \]如果 \(P_{i-1}^TAp_i=0\) , 即 \(\forall j<i,p_j^TAp_i=0\) ，则上式变为
\[\phi(x_i+ap_i)=\phi(x_i)-ap_i^Tr_0+{a^2\over2}p_i^TAp_i\tag{7} \]注意到现在 \(a\) 和 \(y_{i-1}\) 已经分开了。则现在有
\[\underset{y\in\mathbb{R}^{i+1}}{\arg\min}\space\phi(x_0+P_iy)=[y_{i-1},a]^T \tag{8} \]设 \(\mathcal{K}_i=\left<p_0,\cdots,p_i\right>\) ,其中尖括号表示向量张成的空间，\(\mathcal{K}_i+x_0=\{x_0+x|x\in\mathcal{K}_i\}\) 。则按这种方法迭代的 \(x_i\) 满足
\[x_i=\underset{x\in x_0+\mathcal{K}_{i-1}}{\arg\min}\space\phi(x)\tag{9} \]由此，我们选择 \(p_i\) 为 \(r_i\) 减去其在 \(A\mathcal{K}_{i-1}=\{Ax|x\in\mathcal{K}_{i-1}\}\) 上的投影。记 \((A\mathcal{K}_{i-1})^\perp=\{x|\forall x'\in\mathcal{K}_{i-1},x^TAx'=0\}\), 则有
\[r_i=\gamma_ip_i+z_i,\quad\text{where }z_i\in A\mathcal{K}_{i-1},\space p_i\in(A\mathcal{K}_{i-1})^\perp\tag{10} \]到此我们已经得到了一个算法。但其在每步中需要使用 Gram-Schmidt 方法来求 \(p_i\) 。下面来对其进行优化。先证几个命题。
Proposition 1: \(\mathcal{K}_k=\left<r_0,Ar_0,\cdots,A^kr_0\right>\)
**Proof: ** 使用归纳法。\(k=1\) 时显然成立，设命题对于 \(k\) 成立，注意到
\[r_i-r_0=-AP_{i-1}y_{i-1}\in A\mathcal{K}_{i-1}=\left<Ar_0,\cdots,A^ir_0\right>\tag{11} \]又 \(z_i\in A\mathcal{K}_{i-1}\), 则
\[p_i={1\over\gamma_i}[r_0-(z_i+r_0-r_i)]\in\left<r_0,Ar_0,\cdots,A^ir_0\right>\\ \mathcal{K}_i=\left<\mathcal{K}_{i-1},p_i\right>=\left<r_0,\cdots,A^ir_0\right>\quad\quad\quad\quad\square\tag{12} \]Proposition 2: \(P_{i-1}^Tr_i=0\)
**Proof: ** 由前得知，\(x_i=x_0+P_{i-1}y_{i-1}\) ,
\[\begin{align} y_{i-1}&=\arg\min\phi(x_i)\\ &=\arg\min\left\{\phi(x_0)+x_0^TAP_{i-1}y_{i-1}+{1\over2}y_{i-1}^T P_{i-1}^TAP_{i-1}y_{i-1}-b^TP_{i-1}y_{i-1}\right\}\\ &=(P_{i-1}^TAP_{i-1})^{-1}P_{i-1}^Tr_0\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad(13)\\ P_{i-1}^Tr_i&=P_{i-1}^T(r_0-AP_{i-1}y_{i-1})=0\qquad\qquad\square \end{align} \]**Proposition 3: ** \(\forall j\ne i,\space r_i^Tr_j=0\)
**Proof: ** 由命题二得 \(\forall x\in\mathcal{K}_{i-1},\space x^Tr_i=0\) ，故只需证明 \(r_i\in \mathcal{K}_i\) 。使用归纳法，当 \(i=1\) 时显然，若对 \(i-1\) 成立，则 \(r_i=r_{i-1}-a_{i-1}Ap_{i-1}\) ，又由命题一得 \(A\mathcal{K}_{i-1}\sub\mathcal{K}_i\)，即 \(Ap_{i-1}\in\mathcal{K}_i\)，则 \(r_i\in\mathcal{K}_i\qquad\square\)
由命题三可知，\(\mathcal{K}_i=\left<r_0,\cdots,r_i\right>\), 且 \(r_0,\cdots,r_i\) 构成一组正交基。
未完待续。