几种求解线性方程组的迭代解法
共轭梯度法
共轭梯度法（Conjugate Gradient method, CG）是一种求解线性方程组的迭代解法。
基本思路
假设 \(A\) 为 \(n\) 阶实对称正定矩阵，需求解线性方程组
定义 \(\left(u, v\right) = u^T v\)，假设 \(n\) 个方向向量 \(d_1, d_2, \cdots, d_n\) 满足
则称它们关于 \(A\) 共轭。设 \(A=LL^T\)，则 \(\left\{L^T d_i\right\}\) 构成一组 \(\mathbb{R}^n\) 的正交基，\(\left\{d_i\right\}\) 线性无关。
设 \(x_0\) 为初值，真实解 \(x\) 必能表示为 \(x=x_0+\sum\limits_{i=1}^n \alpha_i d_i\) 的形式。
定义第 \(k\) 步迭代的解 \(x_k=x_0+\sum\limits_{i=1}^k \alpha_i d_i\)，误差 \(e_k=x-x_k=\sum\limits_{i=k+1}^n \alpha_i d_i\)，残差 \(r_k=b-Ax_k=Ae_k\)。
则有
即
这种方法称为共轭方向法（Conjugate Direction Method）。
为得到向量集 \(\left\{d_k\right\}\)，假设有 \(n\) 个线性无关向量 \(\left\{u_i\right\}\)，可令 \(d_1=u_1, d_k=u_k+\sum\limits_{i=1}^{k-1}\beta_{ki}d_i\)。
为使 \(\left\{d_k\right\}\) 关于 \(A\) 共轭，其中 \(\beta_{ki}\) 应满足
若 \(i\le j\)，则
由
得 \(i\le j\) 时
若取 \(u_k=r_{k-1}\)，则
特别地，不需要担心 \(u_k=r_k=0\) 导致 \(\left\{u_i\right\}\) 不再线性相关的问题，因为 \(r_k=0\) 时已收敛。
\(k=i+1\) 时，
所以
伪代码如下：
r[0] = b - A * x[0]
d[1] = r[0]
k = 0
while 未达到收敛
k += 1
alpha[k] = (r[k-1], r[k-1]) / (d[k], A * d[k])
x[k] = x[k-1] + alpha[k] * d[k]
r[k] = r[k-1] - alpha[k] * A * d[k]
d[k+1] = r[k] + (r[k], r[k]) / (r[k-1], r[k-1]) * d[k]
收敛性分析
定义向量的 A-范数 \(\lVert x\rVert_A = \sqrt{\left(x, Ax\right)}\)，则
即：在 Krylov 子空间 \(x_0+\left\{r_0, Ar_0, \cdots, A^{k-1} r_0\right\}\) 中，\(x_k\) 是使得 \(\lVert e_k\rVert_A\) 最小的点。由于 \(d_k\) 彼此正交，\(x_k\) 同时也是子空间内使得 \(\varphi\left(x\right)\) 最小的点。
设 \(A\) 的特征值为 \(\left\{\lambda_i\right\}\)，由 A 对称正定得 \(\lambda_i\gt 0\)。对应的归一化的特征向量为 \(\left\{v_i\right\}\)（相同特征值算多次），\(r_0=\sum\limits_{i} \zeta_i v_i\)，则
其中 \(P_k\) 代表一个 \(k\) 次多项式，第 \(i\) 项系数为 \(\alpha_i\)（常数项为 \(1\)）。
若 \(P_k\) 取某一确定多项式，则可得到 \(\lVert e_k\rVert_A^2\) 的某一确定上界。
设特征值最大为 \(\lambda_1\)，最小为 \(\lambda_2\)，矩阵条件数 \(\kappa=\dfrac{\lambda_1}{\lambda_2}\)。
\(k\) 阶切比雪夫多项式
若取
则 \(P_k\left(0\right)=1\) 满足常数项为 \(1\) 的条件，且
可得
即在A-范数下收敛阶不超过 \(\dfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\)。
双共轭梯度法
双共轭梯度法（Biconjugate Gradient mothod, BiCG）可解决 \(A\) 是非奇异实方阵的情况。
假设 A 为 \(n\) 阶非奇异实方阵，需求解两个线性方程组
与 \(A\) 对称正定的情形类似，假设存在两组方向向量 \(\left\{d_1, d_2, \cdots, d_n\right\}; \left\{d_1^{\star}, d_2^{\star}, \cdots, d_n^\star\right\}\) 满足
设 \(x=x_0+\sum\limits_{i=1}^n \alpha_i d_i, x^{\star} = x_0^{\star} + \sum\limits_{i=1}^n \alpha_i^{\star} d_i^{\star}\)，\(x_k, x_k^{\star},e_k, e_k^{\star}, r_k, r_k^{\star}, u_k, u_k^{\star}, \beta_{ki}, \beta_{ki}^{\star}\) 定义同上。则
仍令 \(u_k = r_{k-1}, u_k^{\star} = r_{k-1}^{\star}\)，则 \(k=i+1\) 时
\(k\neq i+1\) 时 \(\beta_{ki} = \beta_{ki}^{\star}=0\)。
正常情况下只需要一个方程组的解，不需要得到 \(x^{\star}\)，可写出伪代码如下：
r1[0] = b - A * x[0]
choose r2[0] such that (r1[0], r2[0]) ≠ 0
d1[1] = r1[0] d2[1] = r2[0]
k = 0
while 未达到收敛
k += 1
alpha[k] = (r1[k-1], r2[k-1]) / (d2[k], A * d1[k])
x1[k] = x1[k-1] + alpha[k] * d1[k]
r1[k] = r1[k-1] - alpha[k] * A * d1[k]
r2[k] = r2[k-1] - alpha[k] * A^T * d2[k]
beta[k] = (r1[k], r2[k]) / (r1[k-1], r2[k-1])
d1[k+1] = r1[k] + beta[k] * d1[k]
d2[k+1] = r2[k] + beta[k] * d2[k]
关于这个方法的收敛性，
In the general case nothing is minimised at all, which can be regarded as a theoretical weakness. In practice however, things turn out to be far less serious than we would expect, as numerical experiments have shown.
把 \(Ax=b\) 更改为 \(A^T Ax=A^T b\) 再应用共轭梯度法可得到一个收敛阶不超过 \(\dfrac{\kappa-1}{\kappa+1}\) 的算法，但是太慢了。
共轭梯度二乘法
共轭梯度二乘法（Conjugate Gradient Squared method, CGS）是对双共轭梯度法的改进。
在双共轭梯度法中，可将 \(d_k, d_k^{\star}, r_k, r_k^{\star}\) 表示为
其中 \(D_k, R_k\) 是次数等于 \(k\) 的多项式，令 \(i=k-1\) 时 \(\beta_{ki}=\beta_i\)，则 \(D_k, R_k\) 满足递推关系
设
则 \(p_k, q_k, g_k, u_k\) 满足递推关系
由 \(\left(r_k^{\star}, r_k\right) = \left(R\left(A^T\right)r_0^{\star}\right)^T R_k\left(A\right)r_0 = \left(r_0^{\star}, g_k\right), \left(d_k^{\star}, Ad_k\right) = \left(r_0^{\star}, AD_{k-1}^2\left(A\right)r_0\right)\) 可得 \(\alpha_k, \beta_k\) 的表达式
在双共轭梯度法中 \(x_k = x_{k-1} + \alpha_k d_k\)，而因为 \(\lVert R_k\left(A\right)r_0\rVert_2\) 大概率比 \(\lVert R_k^2\left(A\right)r_0\rVert_2\) 大，所以此处如果用
更新 \(x_k\) 则大概率可以获得更快的收敛速度。
伪代码如下：
q = 0
p = g = u = A - b * x
choose r such that (r, A - b * x) ≠ 0
rg = (r, g)
k = 0
while 未达到收敛
k += 1
alpha = rg / (r, A * p)
q = u - alpha * A * p
g -= alpha * A * (u + q)
x += alpha * (u + q)
beta = (r, g) / rg
rg *= beta
u = g + beta * q
p = u + beta * (q + beta * p)
稳定双共轭梯度法
稳定双共轭梯度法（Biconjugate Gradient Stabilized method, BiCGSTAB）是对共轭梯度二乘法的改进。
基本思路
观察双共轭梯度法，可以发现实际上只需要满足
这个弱一些的条件即可。
设
其中 \(D_k, R_k, D_k^{\star}\) 为次数等于 \(k\) 的多项式。
由
可得
由条件得
且
因此对任意次数小于 \(k\) 的多项式 \(S\)，
即 \(i\lt k-1\) 时，\(\beta_{ki}=0\)，可令 \(\beta_i = \beta_{ki}\)。
因此 \(D_k^{\star}\) 实际上可以取任意的 \(k\) 次多项式，而在共轭梯度二乘法中 \(D_k^{\star}=D_k\) 并没有很好地利用这一性质。
如果令 \(D_k^{\star}\left(A\right) = D_{k-1}^{\star}\left(I-\omega_k A\right)\)，则 \(D_k, R_k, D_k^{\star}\) 满足递推关系
设
则 \(p_k, g_k, s_k\) 满足递推关系
其中
类似共轭梯度二乘法中的操作，如果令 \(x_k = x_{k-1} + \alpha_k p_{k-1} + \omega_k s_k\) 则可以保证 \(g_k = r_0 - Ax_k\) 以获得更快的收敛速度。
如果选取 \(\omega_k=\dfrac{\left(s_k, As_k\right)}{\left(As_k, As_k\right)}\)，则可以令 \(\lVert g_k\rVert_2\) 尽可能小，以提高算法的稳定性。
伪代码如下：
g = p = b - A * x
choose r such that (r, A - b * x) ≠ 0
while 未达到收敛
k += 1
alpha = (r, g) / (r, A * p)
s = g - alpha * A * p
omega = (s, A * s) / (A * s, A * s)
g = s - omega * A * s
x += alpha * p + omega * s
beta = (r, g) / (omega * (r, A * p))
p = g + beta * (p - omega * A * p)
带有预处理的稳定双共轭梯度法
预处理过程（preconditioning）可用于改善矩阵的收敛性。
假设用某种预处理方法得出 \(A\approx K_1K_2\)，其中 \(K_1, K_2\) 可实现快速求逆，对方程 \(\left(K_1^{-1}AK_2^{-1}\right)x=K_1^{-1}b\) 应用稳定双共轭梯度法，并设
则 \(p_k', g_k', s_k'\) 满足递推关系
其中
如果选取 \(\omega_k=\dfrac{\left(s_k, AK_2^{-1}K_1^{-1}s_k\right)}{\left(AK_2^{-1}K_1^{-1}s_k, AK_2^{-1}K_1^{-1}s_k\right)}\)，则可以令 \(\lVert g_k\rVert_2\) 尽可能小，以提高算法的稳定性。
伪代码如下：
g = p = b - A * x
choose r such that (r, A - b * x) ≠ 0
while 未达到收敛
k += 1
p0 = K2 ^ (-1) * K1 ^ (-1) * p
alpha = (r, g) / (r, A * p0)
s = g - alpha * A * p0
s0 = K2 ^ (-1) * K1 ^ (-1) * s
omega = (s, A * K2 ^ (-1) * K1 ^ (-1) * s0) / (A * K2 ^ (-1) * K1 ^ (-1) * s0, A * K2 ^ (-1) * K1 ^ (-1) * s0)
g = s - omega * A * s0
x += alpha * p0 + omega * s0
beta = (r, g) / (omega * (r, A * p0))
p = g + beta * (p - omega * A * p0)
广义最小残量法
广义最小残量法（Generalized Minimal Residual method, GMRES）是一种求解线性方程组的更简易的解法。
假设 A 为 \(n\) 阶非奇异实方阵，需求解线性方程组
设 \(x_0\) 为初值，\(r_0 = b - Ax_0\)，前 \(k\) 次迭代中对向量集
正交化可得到 \(k\) 个正交归一的向量
设 \(x_k = x_0 + \sum\limits_{i=1}^k \alpha_i q_i\)，则 \(r_k = r_0 - \sum\limits_{i=1}^k \alpha_i Aq_i\)，假设
则 \(AQ_k = Q_{k+1}H_k\)，其中 \(H_k\) 为 \(\left(k+1\right)\times k\) 的上 Hessenberg 矩阵（若 \(k+1\) 个向量线性相关则令最后一行为 \(0\)，并且在这轮迭代后结束即可）。
残差的二范数
其中 \(e_1\) 为 \(n+1\) 阶向量，第一个分量为 \(1\)，其余为 \(0\)。
则只需求解一个最小二乘问题即可得到 \(\alpha\) 和 \(x_k\)。
收敛性：假设 \(A\) 可对角化为 \(A=X\Lambda X^{-1}\)，第 \(k\) 轮残差 \(r_k\)，\(P_k\) 表示复数域常数项为 \(1\)、度数不超过 \(k\) 的多项式的集合，则类似共轭梯度法
如果 \(A\) 的对称部分 \(M = \dfrac{1}{2}\left(A + A^T\right)\) 为正定矩阵，则 \(A\) 的所有特征值的实部均大于 \(0\)。此时根据复杂的推导可以得到
当 \(A\) 本身即对称正定时，直接取 \(P_k\) 为切比雪夫多项式即可得到
广义最小残量法每步的运算次数与之前的迭代次数成正比，为平衡时间与收敛性可能需要定期重构。
