线性方程的迭代方法

本文将介绍求解线性方程· \(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种常用的选择

1. \(A\) 的对角元 \(D\) (Jacobi 格式)，其具体迭代为：

\[x^{(k+1)}_i={1\over a_{ii}}\left(-\sum_{j\ne i}a_{ij}x^{(k)}_j+b_i\right)\tag{2} \]

2. \(A\) 的下三角元 \(D+L\) (Gauss-Siedel 格式)，其具体迭代为：

\[x^{(k+1)}_i={1\over a_{ii}}\left(-\sum_{j<i}a_{ij}x^{(k+1)}_j- \sum_{j>i}a_{ij}x^{(k)}_j+b_i\right)\tag{3} \]

3. 以上二者线性组合 \({1\over \omega}D+L\) (SOR 格式)，其具体迭代为：

\[x_i^{(k+1)}={\omega\over a_{ii}}\left(b_i-{\omega-1\over\omega}x_i^{(k)}-\sum_{j<i}a_{ij}x_j^{(k+1)}-\sum_{j>i}a_{ij}x_j^{(k)}\right)\tag{4} \]

共轭梯度法(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\) 构成一组正交基。

未完待续。

标签：phi,迭代,qquad,线性方程,tag,mathcal,方法,left
来源： https://www.cnblogs.com/songlun/p/15309927.html