Linear Algebra Done Left - A Taste of Iterative Methods for Linear Systems

这个系列的第二篇从上学期鸽到现在。。。

线性代数最基本的问题就是解线性系统了，但是这也是个很不$\text{trivial}$的问题。

求解线性系统，主要有两类方法，直接求解（分解）或者迭代方法。对于大型稀疏矩阵而言，两种方法各有优势，都有不少深入的文献和代码可以学习。这次简要介绍求解线性系统的迭代方法，直接方法留到这个系列的第三篇（又不知道会鸽到什么时候了，大悲）。

迭代法是一个很大的话题，也有一些专著，比如$\text{Yousef Saad}$($\text{GMRES}$方法的发明者)的书$\text{Iterative Methods for Sparse Linear Systems}$。这次不仅参考了这本书，同时也参考了徐树方的数值线性代数。当然，由于内容丰富不可能一下子学完，这篇小水文也会随着我的继续学习长期更新写完共轭梯度就咕咕。

$\text{1. 古典迭代方法}$

$\text{1.1. 单步线性定常迭代}$

这些是最为基本的方法，基本上都满足下面这样的形式:

$x_{n+1}=Mx_n + g$

典型代表是$\text{Jacobi}$迭代和$\text{Gauss-Seidel}$迭代，这两种迭代方法在此不详细叙述了，但是下面会给出代码。

重点是对于一般的$M$的讨论。在$M$满足什么条件的时候，求解过程收敛？收敛多快？误差如何？这些都是值得讨论的问题。

$\text{1.1.1. 收敛性和误差估计}$

收敛的充要条件是$\rho(M)<1$。这个很容易理解，但是实际中并不好用，因为估计谱半径是很难的。所以我们给出一个判断收敛的充分条件以及一个误差估计:

$\text{Thm.}$若$M$的范数为$q$，则$q<1$时收敛，且有误差估计:

$$||x_k - x_*||\le \frac{q^k}{1-q}||x_1-x_0||$$

稍微放缩一下还是很好证明这个结论的。

利用这个结论，我们就可以得出一个对于一般矩阵判断$Jacobi$迭代和$G-S$迭代是否收敛的条件，但是证明和结论都挺复杂的，这里不列出了。我们对一些更特殊的矩阵给出更简单也更好的结论。

$\text{Thm.}$若$A$对称且对角元均为正数，则$\text{Jacobi}$迭代法收敛的充要条件是$A$和$2D-A$都正定。

$\text{Proof.}$设$\text{Jacobi}$迭代的迭代矩阵为$B$，由于$A$的对角元均为正，我们取其对角线矩阵$D$，则可以证明

$$B=D^{-\frac{1}{2}}(I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}})D^{\frac{1}{2}}$$

由这个式子我们可以知道$B$相似于一个对称矩阵，所以特征值均为实数。要证明$B$的谱半径小于$1$，一个间接的方法是证明$I-B$和$I+B$都是正定的（这是个充要条件，是比较容易证明的），这又可以和$A$和$2D-A$的正定性联系起来。

$\text{Thm.}$若$A$是严格对角占优的或者不可约对角占优的，则两种迭代法都收敛。

稍微解释一下可约的概念，可约的意思就是存在排列矩阵$P$使得$PAP^T=\begin{bmatrix}A_{11}&0\\A_{12}&A_{22}\end{bmatrix}$的形式。

$\text{Proof.}$

首先，严格对角占优矩阵和不可约对角占优矩阵都是可逆的。第一个证明是经典习题，第二个证明略去（

由对角占优性，又可以直接推出正定，证法是给$A$加上$\lambda I(\lambda>0)$然后证明不可逆。

最后来考虑两种迭代法的收敛性，使用反证法来证明迭代矩阵所有特征值的绝对值均小于$1$即可。

$\text{1.1.2. 收敛速度}$

设第$k$轮迭代的误差是$y_k$，则容易知道$||y_k||\le ||M^k||||y_0||$，这个结论很好算，把$x_k$和精确解$x_*$都代入定义即可得。

容易知道迭代的误差至少是指数衰减的，所以我们可以通过对数来定义量化收敛速度。定义$R_k(M)=\frac{-\ln||M^k||}{k}$为$k$次迭代的平均收敛速度。

这样，我们比较$R_k$就可以知道什么样的迭代法收敛得快，当然这是考虑$k$轮迭代的。对$R_k$取极限就可以得到渐近收敛速度$R_{\infty}$，这是我们关注的重点。

$\text{Thm.}$ $R_{\infty}=-\ln \rho(M)$

可以当作矩阵分析或者数学分析的简单习题看起来很简单，虽然我不见得能写出来，但是感性理解起来这确实很合理

$\text{1.2. SOR迭代}$

中文名为超松弛迭代法，说白了，就是做$G-S$迭代的时候，把$G-S$的计算结果和上一轮迭代结果做因子为参数$\omega(\omega>1)$的线性插值。$\omega < 1$时叫低松弛迭代，$\omega = 1$的时候自然就是$G-S$迭代。$\omega$的学名叫做松弛因子。

但是因子大于1了还能叫插值吗，不管了总之就是这个意思

$\text{1.2.1. 最佳松弛因子的选取}$

对于固定的$\omega$，每一步的迭代矩阵是固定的，所以前面讨论过的一般的单步线性定常迭代的结论还是可以套到这里来。如果希望迭代收敛，我们就需要选取比较好的参数$\omega$来让迭代矩阵的谱半径小于$1$。

$\text{Thm.}$ $\text{SOR}$迭代法收敛的必要条件是$0<\omega < 2$。

为什么这里又可以小于1了，算了也不管了，反正名字都是人瞎起的，意思就是那么个意思

$\text{Proof.}$

关键的一步是推出迭代矩阵的行列式值为$|(1-\omega)^n|$，迭代矩阵谱半径小于$1$，则其行列式必定也小于$1$，解不等式即可。对于特殊矩阵可以有更好的结果。

$\text{Thm.}$若$A$是严格对角占优的或不可约对角占优的，且$\omega\in (0,1)$，则$\text{SOR}$迭代收敛。

证明类似上面。

$\text{Thm.}$若$A$是实对称正定矩阵，则$\text{SOR}$迭代收敛的充要条件为$\omega\in (0,2)$。

$\text{1.2.2. 最佳松弛因子}$

对于不一样的$\omega$，迭代收敛的速度肯定也有区别，我们就需要考虑选取什么样的$\omega$可以使得渐近收敛速度最快。

具体问题具体分析吧。。。

$\text{1.3. Coding \& Example}$

下面的代码取自我目前正在开发的数值线性代数库$\text{NLAX}$（目前就叫这么个难听名字吧）。实现了上面提到的三种古典迭代方法。

class Solver
{
    protected:
        SolverStatus status = Undefined;
    public:
        SolverStatus info() const { return status; }

};
class IterativeSolver : public Solver
{
    protected:
        int maxRound = -1;
        Real conv_eps = EPS;
        uint nRounds = 0;
    public:
        void setEpsilon(Real eps) { conv_eps = eps; }
        void setMaxRound(int _maxRound) { maxRound = _maxRound < 0 ? -1 : _maxRound; }
        uint rounds() const { return nRounds; }
        virtual void solve(const Matrix& A, const Vector& b, Vector& ret) = 0;
        virtual Vector solve(const Matrix& A, const Vector& b)
        {
            Vector ret;
            solve(A, b, ret);
            return ret;
        }
};

class JacobiSolver final : public IterativeSolver
{
    public:
        void solve(const Matrix& A, const Vector& b, Vector& ret) override
        {
            assert(A.isSquare() && A.cols() == b.dim() && ret.dim() == b.dim());
            uint n = b.dim();
            Vector Dinv = A.diag().inv();
            Vector last;
            for(nRounds = 1; maxRound < 0 || nRounds < maxRound; nRounds++)
            {
                last = ret;
                for(int i = 0; i < n; i++)
                    ret[i] = (-A.row(i).dot(last) + b[i]) * Dinv[i] + last[i];
                if((ret - last).L1Norm() < conv_eps) break;
            }
        }
};

class GaussSeidelSolver final : public IterativeSolver
{
    public:
        void solve(const Matrix& A, const Vector& b, Vector& ret) override
        {
            assert(A.isSquare() && A.cols() == b.dim() && ret.dim() == b.dim());
            uint n = b.dim();
            Vector Dinv = A.diag().inv();
            Vector last;
            for(nRounds = 1; maxRound < 0 || nRounds < maxRound; nRounds++)
            {
                last = ret;
                for(int i = 0; i < n; i++)
                {
                    ret[i] = b[i];
                    for(int j = 0; j < i; j++)
                        ret[i] -= A(i, j) * ret[j];
                    for(int j = i + 1; j < n; j++)
                        ret[i] -= A(i, j) * last[j];
                    ret[i] *= Dinv[i];
                }
                if((ret - last).L1Norm() < conv_eps) break;
            }
        }
};

class SORSolver final: public IterativeSolver
{
    private:
        Real w = 0.5;
    public:
        SORSolver() = default;
        void setw(Real _w) { w = _w; }
        void solve(const Matrix& A, const Vector& b, Vector& ret) override
        {
            assert(A.isSquare() && A.cols() == b.dim() && ret.dim() == b.dim());
            uint n = b.dim();
            Vector Dinv = A.diag().inv();
            Vector last;
            for (nRounds = 1; maxRound < 0 || nRounds < maxRound; nRounds++)
            {
                last = ret;
                for (int i = 0; i < n; i++)
                {
                    ret[i] = w * Dinv[i] * b[i] + (1.0 - w) * last[i];
                    for (int j = 0; j < i; j++)
                        ret[i] -= w * A(i, j) * ret[j] * Dinv[i];
                    for (int j = i + 1; j < n; j++)
                        ret[i] -= w * A(i, j) * last[j] * Dinv[i];
                }
                if ((ret - last).L1Norm() < conv_eps) break;
            }
        }
};

当然，其实这里没必要用虚函数，用$\text{CRTP}$在编译期确定接口实现可能更好。但我不是来修仙$\text{C++}$黑魔法的，我是来学数值算法的，所以暂时不管这个了。

不过也得说明一点是我目前并不清楚怎样的实现最快，上面的代码并不是最优的写法。

然后做一个实际例子，这个例子来自徐树方《数值线性代数》第二版书第$136$面的第一个上机习题——离散化求解两点边值问题的数值解。

对这个问题跑了一下三种迭代算法，并且比较了$\text{SOR}$迭代在不同松弛因子下的收敛速度，所有的向量范数全部采用$L_1$范数，与微分方程的精确解进行对比，结果如下。

Comparing Jacobi solver, G-S solver and SOR solver on a tri-diagonal linear system
The matrix is 100 * 100.

1. Iterate for 20 rounds and compare the result.
Errors:
 Jacobi: 24.097786
 Gauss-Seidel: 23.323033
 SOR: 24.072686
Solution errors:
 Jacobi: 29.285151
 Gauss-Seidel: 0.214910
 SOR: 0.636983

2. Iterate to a precision of 1e-5 and compare the rounds taken.
Rounds:
 Jacobi: 27714
 Gauss-Seidel: 8495
 SOR: 21938
Errors:
 Jacobi: 57.799333
 Gauss-Seidel: 57.789314
 SOR: 57.769156
Solution errors:
 Jacobi: 0.000020
 Gauss-Seidel: 0.000010
 SOR: 0.000030

3. Compare the effect of SOR with different omegas.
(1) omega = 0.100000, take 105329 rounds
 Error: 57.607890 Solution error: 0.000191
(2) omega = 0.200000, take 56648 rounds
 Error: 57.708677 Solution error: 0.000090
(3) omega = 0.300000, take 38252 rounds
 Error: 57.742277 Solution error: 0.000057
(4) omega = 0.400000, take 28295 rounds
 Error: 57.759080 Solution error: 0.000040
(5) omega = 0.500000, take 22179 rounds
 Error: 57.769156 Solution error: 0.000030
(6) omega = 0.600000, take 17747 rounds
 Error: 57.775877 Solution error: 0.000023
(7) omega = 0.700000, take 14728 rounds
 Error: 57.780677 Solution error: 0.000019
(8) omega = 0.800000, take 12032 rounds
 Error: 57.784270 Solution error: 0.000015
(9) omega = 0.900000, take 10175 rounds
 Error: 57.787070 Solution error: 0.000012
(10) omega = 1.000000, take 8538 rounds
 Error: 57.789316 Solution error: 0.000010
(11) omega = 1.100000, take 7072 rounds
 Error: 57.791148 Solution error: 0.000008
(12) omega = 1.200000, take 5980 rounds
 Error: 57.792674 Solution error: 0.000007
(13) omega = 1.300000, take 4911 rounds
 Error: 57.793966 Solution error: 0.000005
(14) omega = 1.400000, take 4013 rounds
 Error: 57.795073 Solution error: 0.000004
(15) omega = 1.500000, take 3187 rounds
 Error: 57.796041 Solution error: 0.000003
(16) omega = 1.600000, take 2464 rounds
 Error: 57.796878 Solution error: 0.000003
(17) omega = 1.700000, take 1791 rounds
 Error: 57.797625 Solution error: 0.000002
(18) omega = 1.800000, take 1179 rounds
 Error: 57.798294 Solution error: 0.000001
(19) omega = 1.900000, take 564 rounds
 Error: 57.798909 Solution error: 0.000001

可以发现，$G-S$迭代方法确实比$\text{Jacobi}$方法快很多，如果松弛因子没选好，$\text{SOR}$的收敛速度很慢，但是选好了就非常快。经过我更精细的测试，$\omega$取$1.939$左右的时候收敛最快，在这之后收敛速度急剧变慢，在$\omega=2$的测试中，求解器进入了死循环，和理论分析一样没有收敛。

$\text{2. 对一般的投影方法的讨论}$

通过讨论一般的投影方法，能够让我们对于接下来的更多方法有更好的理解。

$\text{3. Krylov子空间方法}$

$\text{Krylov}$子空间有着形如$p(A)v$的基底，原理其实就是使用$p(A)$去近似$A^{-1}$。

我们定义$m$维的$\text{Krylov}$子空间为$\mathcal{K}_m(A,r_0)=<r_0, Ar_0,..., A^{m-1}r_0>$

$r_0$是初始猜测$x_0$的残向量。$\text{Krylov}$子空间方法其实就是在上述的投影法中选取$\mathcal{K}$为$\mathcal{K}_m$，不同方法对于$\mathcal{L}$的选取以及预优的策略则各有各的不同。按照前面一节的理论，选取$\mathcal{L}$为$\mathcal{K}_m$和$A\mathcal{K}_m$自然是一种可行的方案。这就是$\text{Saad}$书中讲到的第一类$\text{Krylov}$子空间方法。但是接下来不会讨论那么多，主要还是讲共轭梯度法，其他的视情况补充。

$\text{3.1. 剧痛版共轭梯度法介绍}$

这大概是最著名的求解对称正定系统的算法了。徐树方的书上单独用了一个章节介绍这个算法。这个算法十分重要，我们从多个角度对其进行理解，这里先按照$\text{Saad}$的思路来。

之所以标题取成这样，是为了和网上流传已久的“无痛版共轭梯度法介绍”反着来。无痛版主要从优化角度切入，这也是徐树方书上的讲法。但是$\text{Saad}$作为数值计算大师，他的思路似乎和正常人不太一样学数学的有正常人吗，这里只能尽力保留其中（我认为）的主干思路。

不过对于共轭梯度的理解角度有很多，所以后面如果我没鸽肯定也会参考无痛版，从优化角度再次理解这个算法。

$\text{3.1.1. Arnodi和Lanczos——构造Krylov子空间的正交基}$

$\text{Saad}$的书上讲法是由一般到特殊，按照他的讲法，共轭梯度法其实就是许多更一般化的算法的特例，这样的顺序不太友好，逻辑大概是这么个逻辑：

$\text{Arnoldi算法}\stackrel{A\text{对称}}{\longrightarrow}\text{Lanczos算法}\stackrel{\text{LU分解}}{\longrightarrow}\text{D-Lanczos算法}\stackrel{\text{简化}}{\longrightarrow}\text{共轭梯度法} $

书上说，共轭梯度法还可以视作是$A$对称时$\text{DIOM}$方法的变体，将从$\text{IOM}$方法推出$\text{DIOM}$方法的过程应用于$\text{Lanczos}$算法就得到了共轭梯度法。然而我都不会

上面扯了很多，但是显然没必要全都懂。。。我很怀疑这个老逼登写书的时候到底有没有想让人看下去这种讲法大概让人还没看到共轭梯度就已经被前面的各种一般化的复杂算法给劝退了。这里简单梳理一下各个方法之间的逻辑，最终综合各个方法中得到的结果，得到共轭梯度法。这里当然不会详述所有的算法，但是会抽出其中重要的过程和结论来得到共轭梯度的结果。

$\text{Arnoldi}$方法可以构造生成$\text{Krylov}$子空间的一组单位正交基$\{v_1, v_2, ...,v_n\}$。大致过程是这样的:

void Arnoldi(const Matrix& A, const Vector& v, Vector V[], Vector AV[], Vector w[], Matrix& H, int m)
    {
        assert(A.isSquare());
        V[0] = v.normalized();
        AV[0] = A * V[0];
        Vector sum(A.rows());
        for(int j = 0; j < m; j++)
        {
            for(int i = 0; i <= j; i++)
                H(i, j) = V[i].dot(AV[j]);
            sum += H(j, j) * V[j];
            w[j] = AV[j] - sum;
            H(j + 1, j) = w[j].L2Norm();
            if(isZero(H(j + 1, j))) return ;
            if(j < m - 1) 
            {
                V[j + 1] = w[j] / H(j + 1, j);
                AV[j + 1] = A * V[j + 1];
            }
        }
    }

随手写的，没打算实装在$\text{NLAX}$里，仅作为演示用途所以写得很不讲究。

上述过程迭代$m$步，生成了$\mathcal{K}_m(A, v_1)$的一组单位正交基。正交性是由构造过程保证的，这和$\text{Gram-Schimidt}$正交化很相似，用归纳法很容易验证。利用归纳法同样很容易验证生成的向量都是$p(A)v_1$的形式，这就保证了生成的$v_1, v_2, ..., v_m$是$\text{Krylov}$子空间的基。

注意，上述代码中下标从$0$开始，而理论中的下标都从$1$开始。

记矩阵$V_m = (v_1, v_2, ... v_m)$。$\text{Arnoldi}$方法给了我们重要的结论：

$$AV_m = V_mH_m + w_me_m^T = V_{m + 1}\bar{H}_m$$

$$V_m^TAV_m = H_m$$

第一个式子是由构造过程得到的，左乘$V_m^T$并利用正交性就得到了第二个式子。

这里面的$H_m$是个上海森堡矩阵应该不需要我强调断句是上 海森堡而不是上海 森堡吧，对于上海森堡矩阵这里不多说，讲特征值的时候再说好了，又一个大坑，因为我们并没怎么利用它的性质。

对于$\text{Lanczos}$方法，它的本质就是将$\text{Arnoldi}$算法用于一个实对称矩阵，这种情况下，矩阵的上海森堡形一定是一个对称三对角线矩阵$T_m$。采用$\text{Saad}$书上的记号来标记三条对角线上的值（实际上真正不同的只有两条）：

$$\alpha_j = h_{jj}, \beta_j = h_{j-1, j} = h_{j, j - 1}$$

这个结论很容易验证，$T_m = V_m^TAV_m$，$A$对称则$T_m$对称，$T_m$同时也是上海森堡形，那么它必然是三对角线矩阵。

于是我们可以简化上面的$\text{Arnoldi}$算法得到$\text{Lanczos}$算法。其实只要把内层循环去掉变成两步运算就可以了。

beta[0] = 0;
for(int j = 0; j < m; j++)
{
    if(j) w[j] = A * v[j] - beta[j] * v[j - 1];
    else w[j] = A * v[j];
    alpha[j] = w[j].dot(v[j]);
    w[j] = w[j] - alpha[j] * v[j];
    beta[j + 1] = w[j].L2Norm();
    if(isZero(beta[j + 1])) return;
    v[j + 1] = w[j] / beta[j + 1];
}

$\text{3.1.2. FOM，IOM和DIOM迭代——一般化的求解}$

现在回到解方程组的问题上来。按照投影法的思路，我们把$\mathcal{K}(A,r_0)$作为$\mathcal{K}$和$\mathcal{L}$，我们我们实际上就是要求解$x_0+\mathcal{K}_m$当中的一个$x_m$，使得

$$b-Ax_m\perp \mathcal{K}_m$$

根据$\text{Arnoldi}$方法的结果，我们知道$v_1 = \frac{r_0}{||r_0||_2}$，不妨令$\beta = ||r_0||_2$

这样$V_m^T r_0 = \beta V_m^Tv_1=\beta e_1$。我们知道$V_m$是$\mathcal{K}$的一个基底， 所以$\mathcal{K}$中的向量一定可以表达成$V_m y$的形式。所以把$x_m$写成

$$x_m = x_0+V_m y_m$$

根据$V_m^T AV_m = H_m$，结合上述结论我们可以算出来（读者自证不难）

$$y_m = \beta H_m^{-1} e_1$$

基于$\text{Arnoldi}$的迭代构造，我们在此基础上就能构造迭代求解了，这就是$\text{FOM}$算法，这里就不写了。

下面一个重要的推论证明，对于$\text{FOM}$迭代，每一步的残向量相互正交。

$\text{Prop. }b-Ax_m = -\beta_{m+1}e_m^Ty_mv_{m+1}$

$\text{Proof.}$

$$b-Ax_m = b-A(x_0 + V_my_m)=r_0 - AV_my_m=\beta v_1-V_mH_my_m-\beta_me_m^Ty_mv_{m+1}$$

$$=\beta v_1-V_m\beta e_1-\beta_me_m^Ty_mv_{m+1}=-\beta_me_m^Ty_mv_{m+1}$$

$\text{Q.E.D.}$

$v_i$是相互正交的，自然残向量也是相互正交的。

在$\text{FOM}$的基础上，我们不必知道$\text{IOM}$算法和$\text{DIOM}$算法具体有什么新花样。只需要知道，$\text{D}$代表“直接”。“直接”的意思，就是做$\text{LU}$分解，接下来我们会把这个做法用到$\text{Lanczos}$算法里。

不管$H_m$长成什么样，我们都能对它做$\text{LU}$分解，不管$U$长成什么样，上海森堡矩阵的结构决定了$L$肯定很好，除去主对角线全是$1$以外，它只有一条次对角线$(l_{21}, l_{32}, ..., l_{n, n-1})$非零。

所以我们写出来$x_m = x_0 + \beta V_mU_m^{-1}L_m^{-1}e_1$并且令$P_m = V_mU_m^{-1},z_m = \beta L_m^{-1}e_1$就有$x_m = x_0+P_m z_m$

由于$L_m$的结构特殊，可以写出$z_m = \begin{bmatrix} z_{m-1} \\ \zeta_m \end{bmatrix}$，其中$\zeta_m = -l_{m, m-1}\zeta_{m-1}$。如果我们定义$P_m = (p_1, p_2, ..., p_m)$，展开一下就得到$x_m = x_{m-1} + \zeta_m p_m$，这就是$\text{DIOM}$的迭代式了。当然这个不止对$\text{IOM}$是成立的。

$p_m$好不好算取决于$P_m$好不好算，而在$A$对称的情形下，$P_m$就非常好算了，这样离共轭梯度法已经很近了。

$\text{3.1.3. 从D-Lanczos算法到共轭梯度法}$

书接前两回，我们已经知道了$H_m$是对称三对角矩阵，还对它写出了$\text{LU}$分解。其中的$U_m$矩阵自然也只有两条对角线，不妨记$\eta_i = u_{ii}, \beta_i = u_{i - 1, i}$

这样我们要迭代求$p_m$就非常简单了，$p_m = \eta_m^{-1}[v_m - \beta_mp_{m - 1}]$

其中$\eta_m = \alpha_m - \frac{\beta_m^2}{\eta_{m-1}}$

我们称其为$\text{D-Lanczos}$算法，这就是把$\text{LU}$分解以及迭代用于$\text{Lanczos}$算法的成果。

最后我们再加上一个推论，就得到了最终的共轭梯度法，这也是这个算法的名字的由来，所谓共轭，就是$p_iAp_j = 0(i\neq j)$。。

$\text{Prop.}$所有$p_i$关于$A$-共轭。

$\text{Proof.}$

$P_m^TAP_m = U_m^{-T}V_m^TAV_mU_m^{-1}=U_m^{-T}T_mU_m^{-1}=U_m^{-T}L_m$

注意$U_m^{-T}L_m$首先是下三角阵，然后它是对称的，所以必然是对角阵。

$\text{Q.E.D.}$

前面已经推出了残向量的正交性，下面就来做最后的化简。利用正交性和共轭性，我们使用更简单的方式计算迭代。

根据前面我们有$x_{j+1}=x_j+\zeta_jp_j$，因此$r_{j+1}=r_j - \zeta_jAp_j$

由正交性，我们有$r_j^T(r_j-\zeta_jAp_j)=0$从而解出

$$\zeta_j=\frac{r_j^Tr_j}{r_j^TAp_j}$$

前面注意到$p_j$可以由$v$和$p_{j-1}$线性递推，由$r$和$v$共线我们就知道，如果不考虑长度只考虑方向（注意我们只需要共轭性），$p_{j+1}$也可以由$r_{j+1}$和$p_j$线性递推。即

$$p_{j+1} = r_{j+1} + \gamma_jp_j$$

首先把这个式子代入$\zeta$，就求出$\zeta = \frac{r_j^Tr_j}{p_j^TAp_j}$，这是最好的形式了。

然后类似于上面的做法，利用共轭性，可以求出

$p_{j+1}=r_{j+1} + \frac{||r_{j+1}||_2^2}{||r_{j}||_2^2}p_j$

所以就得到下面最经典的共轭梯度求解器了。

$\text{3.1.4. Coding \& Example}$

下面的求解器是根据$\text{Saad}$书上的伪代码写成的，是最经典的形式，书上其它的变体这里就不管了。

class CGSolver : public IterativeSolver
{
    public:
        void solve(const Matrix& A, const Vector& b, Vector& ret) override
        {
            assert(A.isSymmetric());
            assert(A.isSymmetric() && A.cols() == ret.dim());
            Vector r(b - A * ret);
            Vector p(r);
            Real alpha, beta, rNormSqrCache = r.L2NormSqr();
            Vector Ap = A * p;
            for(nRounds = 1; maxRound < 0 || nRounds < maxRound; nRounds++)
            {
                alpha = rNormSqrCache / p.dot(Ap);
                ret.saxpy(p, alpha);
                if(alpha * p.L1Norm() < conv_eps) break;
                r.saxpy(Ap, -alpha);
                beta = r.L2NormSqr() / rNormSqrCache;
                rNormSqrCache = r.L2NormSqr();
                p.scadd(r, beta);
                Ap = A * p;
            }
        }
};

代码并不长，但这背后的故事实在是不容易（

同样说明的一点是我并不知道怎样的写法是最好的，按照我的理解，上面主要使用向量的$\text{saxpy}$操作和$\text{scaleAndAdd}$操作，如果这些线性代数基本操作能够做好，那么速度应当是不慢的。

同样做一个测试例子，这里选用徐树方书上第$\text{160}$面的$\text{Hilbert}$矩阵进行测试，并对比$\text{Gauss-Seidel}$迭代法和共轭梯度法，$\text{Jacobi}$迭代法似乎无法收敛，因此不测试，$\text{SOR}$迭代法调参难调，不稳定，因此也不测试。

测试内容是对$200$阶的矩阵进行迭代求解，包括迭代$\text{20}$轮比较求解精度和求解至精度为$10^{-5}$比较迭代轮数。最后为了测试共轭梯度法速度，对其在$2048$阶的矩阵上的表现进行测试。

Comparing classical solvers and conjugate gradient solver on Hilbert matrix
The matrix is 200 * 200.

1. Iterate for 20 rounds and compare the result.
Solution errors:
 Gauss-Seidel: 0.098972
 CG: 0.000001

2. Iterate to a precision of 1e-5 and compare the rounds taken.
Rounds:
 Gauss-Seidel: 30349
 CG: 25
Solution errors:
 Gauss-Seidel: 0.000001
 CG: 0.000000

Now test CG solver on 2048 * 2048 Hilbert matrix.
Start timing...
Stop timing. Taking 924.225000 ms, 51 rounds.
 The solution error is 0.000000

不得不说共轭梯度法的效果比我想得要好太多，$2048$阶的稠密矩阵也能在一秒以内快速求解，这还是在没有任何优化技巧，没有并行也没有预优的情况下，这种效率在以前只学过高斯消元的时候简直是无法想象的。上面的代码优化空间非常大，如果能进行优化，效果一定更优秀。对于稠密矩阵尚且如此，对于单步计算量更小的稀疏矩阵，能够求解的规模更大，表现肯定也会更好，用途更广泛。

$\text{3.1.5. 收敛性和收敛速度}$

证明涉及到数值分析中的切比雪夫多项式，所以这里略去其实是懒得写了，手上都有求解器了谁还看收敛性分析啊

上结论！

$$||x_* - x_m||_A\le 2(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1})^m||x_*-x_0||_A$$

$\kappa$是矩阵在$2$范数下的条件数，也就是绝对值最大特征值和绝对值最小特征值的比值。

所以可以看到，只要是对称正定矩阵，就可以应用共轭梯度法，并且必定收敛。

$\text{Update on 2023.5.7:}$看业内人士说$\text{Saad}$的这本书被认为是这个领域的圣经。。。你们的圣经就这么劝退异教徒是吧

$\text{3.2. 预优}$

根据前面的收敛性分析，如果矩阵比较病态，收敛得会很慢。如果我们能降低矩阵的条件数，那么就能加速收敛。所谓预优，其实就是将求解$Ax=b$转化为求解$M^{-1}A = M^{-1}b$，如果$M$选得好，那么收敛速度自然也就会更快。当然，如果预优做得不好那么反而可能会起到反作用跑得更慢预劣。实际上，我们并不指望$M$和$M^{-1}A$都能显式计算出来，但是我们希望对$M^{-1}$和$A$的操作不能太慢，这样才有可能提速。

预优同样是个极其庞杂的问题，这里也只能介绍一些简单的预优手段，而且实现同样有难度，预优的实现做得不好很可能就是反效果。

这个时候自然会去看看$\text{Eigen}$实现了哪些，来获得重要的参考。。。

答案是$\text{Eigen}$实现了三种预优器，最简单的对角线预优，$\text{ILUT}$预优和不预优。

$\text{3.2.1. ILUT预优}$

$\text{LU}$自然是$\text{LU}$分解，$\text{I}$代表着“不完全”($\text{Incomplete}$)。把$A$写成$A=LU-R$，其中$R$代表残量矩阵。我们希望$R$具有一些特殊的性质，方便我们求出一个好用的近似$LU$来做预优。

$\text{ILU}$这个方法下有太多变种，这里自然就选$\text{Eigen}$用了也说好的$\text{ILUT}$预优。$\text{T}$是$\text{Threshold}$的缩写。

还是那句话，预优要做好并不容易，实际上，为了给出高效的实现，$\text{Saad}$在书中甚至专门用一节来讲了实现细节。

$\text{3.2.1.1. ILU(0)}$

在讨论$\text{ILUT}$之前，首先提一下最简单的$\text{ILU}$方法，也就是$\text{ILU(0)}$。其中$0$代表没有非零项填充。

简单来说，我们对于一种特殊的情况考虑$L,U$和$A$的稀疏模式——有限差分中常见的五点差分格式矩阵。如果我们取$A$的下三角的稀疏模式作为$L$的稀疏模式，取$A$的上三角的稀疏模式作为$U$的稀疏模式，那么$A$和$LU$的稀疏模式大体上就比较接近。但是$LU$中还会多出一些额外的对角线，忽略这些多出来的对角线就行了，也就是把它们全都算在$R$里。

这就是$ILU(0)$的思想：寻找$L,U$使得$A-LU$在$A$的稀疏模式的位置上的元素均为$0$。显然这样的$L,U$并不唯一，因此书上给出了一种常用的构造性算法，这是更一般的算法（此处略去）在$\text{ILU(0)}$情况下的特例。

for(int i = 2; i <= n; i++) {
    for(int k = 1; k <= i - 1; k++) {
        if(a[i][k] == 0) continue;
        a[i][k] /= a[k][k];
        for(int j = k +1; k <= n; j++) {
            if(a[i][j] == 0) continue;
            a[i][j] -= a[i][k] * a[k][j];
        }
    }
}

实际上，更一般的算法也就是把上面代码中的两个非$0$判断都换成判断当前位置是否属于集合$P$，其中$P$代表了我们希望把哪些元素变成$0$。在$\text{ILU(0)}$情况下，$P$就对应于$A$的稀疏模式。经过这个过程，$A$存储的元素就变成了$L$和$U$存在一块儿的模式（因为一个上三角一个下三角）。

至于这样的方法为什么是有效的，讨论很长，此处略去，可以大致上认为这就是对矩阵中的一部分元素做了高斯消元/$\text{LU}$分解。

$\text{3.2.1.2. ILUT原理}$

前面的$\text{ILU(0)}$忽视数值而只考虑稀疏模式，对于实际问题来说这是很要命的。而$\text{ILUT}$一定程度上弥补了这个缺陷。

$\text{ILUT}$算法大概是这样的

for(int i = 1; i <= n; i++) {
    Vector w = A.Row(i);
    for(int k = 1; k <= i - 1; k++) {
        if(w[k] == 0) continue;
        w[k] /= a[k][k];
        DropElement(w[k]);
        if(w[k] != 0) w -= w[k] * U.Row(k);
    }
    Drop(w);
    for(int j = 1; j < i; j++)
        L[i][j] = w[j];
    for(int j = i; j <= n; j++)
        U[i][j] = w[j];
}

其中DropElement(w[k])表示以某种规则将w[k]选择性地置为$0$，Drop同理，规则与DropELement可能不同。

在$\text{ILUT}(p,\tau)$算法中，这个规则是这么选定的：

1. 在DropElement(x)中，如果$|x|<\tau||A_i||$，我们将$x$置为$0$。

2. 在Drop中，首先对每个元素应用上述的同样规则，然后只在$L$和$U$中各自保留除主对角线外最大的$p$个元素。

$\text{3.2.1.3. 实现细节}$

如果算法实现得不好，那么预优这一步很可能做得非常慢，所以需要仔细考虑实现。在$\text{ILUT}$这一步，主要有$3$个地方，如果做得不好可能导致低效。

1. w -= w[k] * U.Row(k)这一步

2. 选择最大的$p$个元素

3. 在第二重循环中以升序按列访问$L$

第一个倒不是啥大麻烦，对w使用稀疏存储或者稠密存储都行。或者可以使用稠密存储，然后记下非零的位置以便后续清零。

对于第二个，使用堆排序或者快排都可以。

第三个问题。。。我还真没想明白这为什么是个问题（