给定 $n$ 阶方阵 $A$ 与系数向量 $\boldsymbol b \in \mathbb R^n$，求解满足

$$A \boldsymbol x = \boldsymbol b$$

的未知向量 $\boldsymbol x$ 的问题被称为线性方程式问题。
本节讨论计算机中如何求解该类问题

线性代数中，如果 $\det A \neq 0$，则可以写出解

$$\boldsymbol x = A^{-1} \boldsymbol b$$

但实际应用上，由于阶数 $n$ 可能非常巨大，使得逆矩阵计算困难。
此时需要考虑基于 Gauss 消元法的求解

如果系数矩阵是一个上三角或者下三角矩阵，那么求解该方程将极其容易
将系数矩阵 $A$ 改写为上三角矩阵 $U$ 与下三角矩阵 $L$ 的乘积，即 $A = LU$，称为 LU 分解
此时，方程可以改写为

$$L(U \boldsymbol x) = \boldsymbol b$$

我们可以很轻松地求解出 $U \boldsymbol x$，进而求解出 $\boldsymbol x$。

并且，我们可以要求下三角矩阵 $L$ 的对角线元素全为 $1$，以此来保证分解的唯一性

接下来让我们看看如何实现 LU 分解

# 无 Pivoting 的 LU 分解

假设我们已经实现了 LU 分解，也就是 $A = LU$
对 $L = (\ell_{ij})$ 进行列向量分解，对 $U = (u_{ij})$ 进行行向量分解，得到

$$A = LU = \begin{pmatrix} \boldsymbol l_1 & \boldsymbol l_2 & \cdots & \boldsymbol l_n \end{pmatrix} \begin{pmatrix} \boldsymbol u_1^T \\ \boldsymbol u_2^T \\ \vdots \\ \boldsymbol u_n^T \end{pmatrix} = \sum_{k=1}^n \boldsymbol l_k \boldsymbol u_k^T$$

我们关注求和项，各自得到

$$\begin{aligned} \boldsymbol l_k \boldsymbol u_k^T &= \begin{pmatrix} 0 \\ 0 \\ \vdots \\ \ell_{kk} \\ \vdots \\ \ell_{nk} \end{pmatrix} \begin{pmatrix} 0 & 0 & \cdots & u_{kk} & \cdots & u_{kn} \end{pmatrix} \\ &= \begin{pmatrix} 0 & 0 & \cdots & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \ell_{kk} u_{kk} & \cdots & \ell_{kk} u_{kn} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \ell_{nk} u_{kk} & \cdots & \ell_{nk} u_{kn} \end{pmatrix} \end{aligned}$$

观察不难得到，每一个第 $k$ 个求和项实际上的贡献在于右下角开始的一个边长 $\ell(k)$ 的方阵
因此，通过这样的逐级分解，就可以得到 LU 分解

具体实现方法如下：
无 Pivoting 的 LU 分解算法

• 关注 $A$ 的第一行第一列，这部分只能由单一的求和项给出，因此可以确定 $\boldsymbol \ell_1$ 与 $\boldsymbol u_1^T$ 的值。只需要将 $\boldsymbol u_1^T$ 指定为 $A$ 的第一行，再相除得到 $\boldsymbol \ell_1$
• 令 $A_1 = A - \boldsymbol \ell_1 \boldsymbol u_1^T$
• 关注 $A_1$ 的第二行第二列，这部分只能由单一的求和项给出，因此可以确定 $\boldsymbol \ell_2$ 与 $\boldsymbol u_2^T$ 的值。只需要将 $\boldsymbol u_2^T$ 指定为 $A_1$ 的第二行，再相除得到 $\boldsymbol \ell_2$
• 以此类推

以上，通过 $n-1$ 次除法，$n(n^2-1)/3$ 次乘法，$n(n^2-1)/3$ 次减法，就可以得到 LU 分解。全体的计算量为 $2n^3/3$，因此时间复杂度为 $O(n^3)$

算法实现

for(k=1; k<=n-1; k++){

    w = 1/a[k][k];

    for(i=k+1; i<=n; i++){

        a[i][k] = a[i][k] * w;

        for(j=k+1; j<=n; j++){

            a[i][j] = a[i][j] - a[i][k] * a[k][j];

        }

    }

}

这样的 LU 分解有一个极其显著的优点：不占用内存

在完成 LU 分解后，就可以快速求解方程了
令 $\boldsymbol y = U \boldsymbol x$，我们先求解 $L \boldsymbol y = \boldsymbol b$
$L$ 是一个下三角矩阵，通过每一次向矩阵的右侧移动消元，可以得到 $\boldsymbol y$ 的值。这个方法叫做 前代 (Forward Substitution)「前進消去」。

此时方程等价于

$$\begin{cases} y_1 & = b_1 \\ \ell_{21} y_1 + y_2 & = b_2 \\ \ell_{31} y_1 + \ell_{32} y_2 + y_3 & = b_3 \\ \vdots \\ \ell_{n1} y_1 + \ell_{n2} y_2 + \cdots + \ell_{n, n-1} y_{n-1} + y_n & = b_n \end{cases}$$

对第二个式子代入得到

$$y_2 = b_2 - \ell_{21} y_1$$

对第三个式子代入得到

$$y_3 = b_3 - \ell_{31} y_1 - \ell_{32} y_2$$

以此类推，我们就可以得到 $\boldsymbol y$ 的值了

算法实现

for(i=1; i<=n; i++){

    s = 0;

    for(j=1; j<=i-1; j++){

        s = s + a[i][j] * y[j];

    }

    y[i] = b[i] - s;

}

接下来求解 $U \boldsymbol x = \boldsymbol y$，$U$ 是一个上三角矩阵，通过每一次向矩阵的左侧移动消元，可以得到 $\boldsymbol x$ 的值。这个方法叫做 回代 (Backward Substitution)「後退代入」。

此时方程等价于

$$\begin{cases} u_{11} x_1 + u_{12} x_2 + \cdots + u_{1n} x_n & = y_1 \\ u_{22} x_2 + \cdots + u_{2n} x_n & = y_2 \\ \vdots \\ u_{nn} x_n & = y_n \end{cases}$$

我们需要反向代入 $x_n = y_n / u_{nn}$ 到倒数第二个式子中，得到

$$x_{n-1} = \frac{y_{n-1} - u_{n-1, n} x_n}{u_{n-1, n-1}}$$

以此类推，我们就可以得到 $\boldsymbol x$ 的值了

算法实现

for(i=n; i>=1; i--){

    s = 0;

    for(j=i+1; j<=n; j++){

        s = s + a[i][j] * x[j];

    }

    x[i] = (y[i] - s) / a[i][i];

}

# Pivoting 下的 LU 分解

LU 分解基于 Gauss 消元，而消元法中要求每一步所用的主元都不能为 $0$
如果出现了这种情况，就需要进行行交换，称为 Pivoting

回顾上述 LU 分解算法中，我们存在系数计算

$$w = \frac{1}{a_{kk}}$$

这个计算显然只能在 $a_{kk} \neq 0$ 的情况下进行，但是实际应用上可能会出现 $a_{kk} = 0$ 的情况。细分上来说可以分为两种

• 在 $a_{kk} = 0$ 的情况下，同一列的其他元素有 $a_{ik}, 1 \leq i \leq n$ 有非零元。这个时候我们可以考虑对 $A$ 执行行交换来解决问题
• 如果整个列都是 $0$，那么不仅没办法做 LU 分解，方程本身都没有唯一解。我们也就不讨论了。

记行交换矩阵为 $P_{ij}$，并简记所有行交换矩阵的乘积为 $P$，那么此时 LU 分解的形式为

$$PA = LU$$

虽然看似没有单纯对 A 执行分解，但是实际上行交换并不会对方程的解造成影响，所以在解方程的视角下是等价的

由于我们需要做除法，即使非零的除数，过小也会导致误差过大。所以原则上我们需要将 $a_{kk}$ 替换为第 $k$ 列中绝对值最大的元素来进行除法计算，这样就可以保证数值稳定性了
逐步交换行，并逐步执行 LU 分解的算法被称为 部分 Pivoting 的 LU 分解

等价地，我们也可以进行列的交换，记列交换矩阵为 $Q_{ij}$，并简记所有列交换矩阵的乘积为 $Q$，那么此时 LU 分解的形式为

$$AQ = LU$$

注意 $Q^{-1} = Q^T$，因此

$$A \boldsymbol x = \boldsymbol b \Leftrightarrow AQ Q^T \boldsymbol x = \boldsymbol b \Leftrightarrow LU Q^T \boldsymbol x = \boldsymbol b$$

从左侧逐层求解即可

如果同时允许交换行和列，那么分解的形式为

$$PAQ = LU$$

这样的分解被称为 完全 Pivoting 的 LU 分解

# LDU 分解

单纯的 LU 分解虽然可以快速求解矩阵，但是并没有体现出 $A$ 的诸多性质，尤其是正定值。
因此，通过对上三角矩阵 $U$ 执行一次分解

$$U = DU'$$

其中 $D$ 是一个对角矩阵，$U'$ 是一个上三角矩阵，并且 $U'$ 的对角线元素全为 $1$，我们就可以得到 (便于理解换回 $U$)

$$A = LDU$$

该分解被称为 LDU 分解

考虑无 Pivoting 下的分解 $A = LU$
如果 $A$ 是一个对称矩阵，那么

$$A^T = (LDU)^T = U^T D L^T$$

结合 $A^T = A$，我们可以得到

$$U = L^T$$

此时，分解形如

$$A = L D L^T$$

该分解被称为 $LD L^T$ 分解

更进一步，若 $A$ 不仅对称，还是一个正定的，这意味着 $D$ 的对角线元素全为正数，因此我们可以将 $D$ 分解为

$$D = \sqrt{D} \sqrt{D}$$

其中 $\sqrt{D}$ 是一个将 $D$ 的对角线元素开平方的对角矩阵，此时分解形如

$$A = L \sqrt{D} \sqrt{D} L^T = (L \sqrt{D})(L \sqrt{D})^T$$

令 $C = L \sqrt{D}$，则分解形如

$$A = CC^T$$

该分解被称为 Cholesky 分解

# 迭代法

上述的分解再回代的方法在数值计算中被称为直接方法，直接方法的优点是可以稳定地求解方程，缺点是当 $n$ 非常大时，计算量过大，尤其是面对一下具有很多零元素的稀疏矩阵时，直接方法的效率会大大降低，因为分解反而会引入很多非零元素

因此，针对稀疏矩阵，我们可以考虑迭代方法来求解方程
迭代方法的基本思想是先给出一个初始解 $\boldsymbol x_0$，然后通过不断迭代来得到更好的解 $\boldsymbol x_1, \boldsymbol x_2, \dots$，直到满足某个收敛条件为止

考虑线性方程组

$$A \boldsymbol x = \boldsymbol b$$

我们将通过一系列变形得到

$$\boldsymbol x = M \boldsymbol x + \boldsymbol c$$

其中 $M$ 是一个矩阵，$\boldsymbol c$ 是一个向量
记 $\boldsymbol x^{k}$ 为第 $k$ 次迭代的解，那么我们通过代入

$$\boldsymbol x^{k+1} = M \boldsymbol x^k + \boldsymbol c$$

就可以反复更新 $\boldsymbol x$ 的值，并且判断误差

$$\boldsymbol r := A \boldsymbol x^k - \boldsymbol b$$

是否足够小

线性方程组问题的迭代法一般有两种

• Gauss-Jacobi 迭代法，这个迭代方法需要占用内存
• Gauss-Seidel 迭代法，这个迭代方法不占用内存，但是由于后者的运算依赖于先前的结果，所以无法并行

我们先来看 Gauss-Jacobi 迭代法
对于矩阵 $A$，我们可以将其分解为一个一个对角矩阵 $D$，一个狭义下三角矩阵 $E$，一个狭义上三角矩阵 $F$ 的和，即

$$A = D + E + F$$

其中成分分别为

$$D = \begin{pmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn} \end{pmatrix},\quad E = \begin{pmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & 0 \end{pmatrix},\quad F = \begin{pmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{pmatrix}$$

通过移项，原方程改写为

$$D \boldsymbol x = \boldsymbol b - E \boldsymbol x - F \boldsymbol x$$

对角矩阵可以简单地求逆，因此

$$\begin{aligned} \boldsymbol x & = D^{-1} (\boldsymbol b - E \boldsymbol x - F \boldsymbol x) \\ & = -D^{-1} (E + F) \boldsymbol x + D^{-1} \boldsymbol b \end{aligned}$$

接下来就只需要不断地更新 $\boldsymbol x$ 的值了

算法实现

while(...){

    for(i=1; i<n; i++){

        double s = b[i];

            for(j=1; j<n; j++){

                if(j != i){

                    s = s - a[i][j] * x[j];

                }

            }

            xnew[i] = s / a[i][i];

    }

    for(i=1; i<n; i++){

        x[i] = xnew[i];

    }

}

这里我们每一次都需要暂时储存一下更新后的 $\boldsymbol x$ 的值到 $x_{\text{new}}$ 中，等到更新完成后再将 $x_{\text{new}}$ 的值赋给 $x$，因此该方法需要占用内存

Gauss-Seidel 迭代法考虑直接使用迭代后的值来更新 $\boldsymbol x$ 的值
对应到方程式中时，我们每次是先迭代出第一个成分 $x_1^{k+1}$ 再利用它去迭代其他成分，所以下三角矩阵 $E$ 对应着迭代后的值，也就是说

$$D \boldsymbol x^{k+1} = \boldsymbol b - E \boldsymbol x^{k+1} - F \boldsymbol x^k$$

汇总整理得到

$$\boldsymbol x^{k+1} = -(D + E)^{-1} F \boldsymbol x^k + (D + E)^{-1} \boldsymbol b$$

从这里就可以看出该算法的一个劣势是需要计算 $(D + E)^{-1}$，这并不如直接计算 $D^{-1}$ 来得简单，但是如下算法展示，他不需要占用内存

算法实现

while(...){

    for(i=1; i<n; i++){

        double s = b[i];

            for(j=1; j<n; j++){

                if(j != i){

                    s = s - a[i][j] * x[j];

                }

            }

            x[i] = s / a[i][i];

    }

}

# 共轭梯度法

共轭梯度法 (Conjugate Gradient Method) 简称 CG 方法
是一种非常有趣的方法，它告诉我们，求解线性方程组问题也可以转为优化问题

设 $A$ 为一个正定对称矩阵（这是优化问题有解的必要条件），定义二次型

$$\varphi(\boldsymbol x) = \frac{1}{2} \boldsymbol x^T A \boldsymbol x - \boldsymbol b^T \boldsymbol x$$

我们想一想如何求解 $\varphi$ 的最小值：梯度计算可得

$$\nabla \varphi(\boldsymbol x) = A \boldsymbol x - \boldsymbol b$$

因此，$\varphi$ 的最小值点 $\boldsymbol x^*$ 必须满足

$$A \boldsymbol x^* = \boldsymbol b$$

这正是我们想要求解的线性方程组问题
而且误差 $\boldsymbol r = \boldsymbol b - A \boldsymbol x = -\nabla \varphi(\boldsymbol x)$ 正是 $\varphi$ 的下降方向

最朴素的想法自然是直接沿着 $\boldsymbol r = -\nabla \varphi(\boldsymbol x)$ 的方向进行梯度下降，设步长为 $\alpha_k$，则更新公式为

$$\boldsymbol x^{k+1} = \boldsymbol x^k + \alpha_k \boldsymbol r^k$$

一维精确线搜索要求：在固定的方向 $\boldsymbol r^k$ 上，找到一个步长 $\alpha$ 使得目标函数最小。令

$$f(\alpha) = \varphi(\boldsymbol{x}^k + \alpha \boldsymbol{r}^k)$$

则条件为

$$\frac{d}{d\alpha} f(\alpha) = 0$$

利用多元函数链式法则求导

$$\begin{aligned} \frac{d}{d\alpha} \varphi(\boldsymbol{x}^k + \alpha \boldsymbol{r}^k) & = \nabla \varphi(\boldsymbol{x}^k + \alpha \boldsymbol{r}^k)^T \cdot \frac{d}{d\alpha} (\boldsymbol{x}^k + \alpha \boldsymbol{r}^k) \\ & = \nabla \varphi(\boldsymbol{x}^k + \alpha \boldsymbol{r}^k)^T \cdot \boldsymbol{r}^k \\ & = (A (\boldsymbol{x}^k + \alpha \boldsymbol{r}^k) - \boldsymbol{b})^T \cdot \boldsymbol{r}^k \\ & = (A \boldsymbol{x}^k - \boldsymbol{b})^T \cdot \boldsymbol{r}^k + \alpha (\boldsymbol{r}^k)^T A \boldsymbol{r}^k \\ & = -(\boldsymbol{r}^k)^T \cdot \boldsymbol{r}^k + \alpha (\boldsymbol{r}^k)^T A \boldsymbol{r}^k \end{aligned}$$

解出

$$\alpha_k = \frac{\boldsymbol r^{kT} \boldsymbol r^k}{\boldsymbol r^{kT} A \boldsymbol r^k}$$

但是我们回顾这个过程；我们沿着一个方向 $\boldsymbol r^k$ 在精确线搜索下，走到了函数的最小值点
此时该点的切线方向就是我们的搜索方向 $\boldsymbol r^k$，法线方向就是下一步搜索方向 $\boldsymbol r^{k+1}$，因此相邻搜索方向正交，即

$$\langle \boldsymbol r^k, \boldsymbol r^{k+1} \rangle = \boldsymbol r^{kT} \boldsymbol r^{k+1} = 0$$

如果我们的搜索范围是一个二维平面，来回两次之后我们就会回到与先前平行的方向。然而我们已经走过了在该方向上的最优点，这就会导致先前的结果被破坏，形成 Z 震荡

回到我们的起点，假设 $A\boldsymbol x = \boldsymbol b$ 的精确解为 $\boldsymbol x^*$，定义每一次迭代的误差

$$\boldsymbol e^k = \boldsymbol x^* - \boldsymbol x^k$$

将 $\boldsymbol x^k = \boldsymbol x^* - \boldsymbol e^k$ 和 $\boldsymbol b = A\boldsymbol x^*$ 代入目标函数 $\varphi(\boldsymbol x^k)$ 中展开：

$$\begin{aligned} \varphi(\boldsymbol x^k) &= \frac{1}{2} (\boldsymbol x^* - \boldsymbol e^k)^T A (\boldsymbol{x}^* - \boldsymbol{e}^k) - (\boldsymbol{x}^*)^T A^T (\boldsymbol{x}^* - \boldsymbol{e}^k) \\ &= \frac{1}{2} (\boldsymbol{x}^*)^T A \boldsymbol{x}^* - (\boldsymbol{x}^*)^T A \boldsymbol{e}^k + \frac{1}{2} (\boldsymbol{e}^k)^T A \boldsymbol{e}^k - (\boldsymbol{x}^*)^T A \boldsymbol{x}^* + (\boldsymbol{x}^*)^T A \boldsymbol{e}^k \\ &= \frac{1}{2} (\boldsymbol{e}^k)^T A \boldsymbol{e}^k - \frac{1}{2} (\boldsymbol{x}^*)^T A \boldsymbol{x}^* \end{aligned}$$

在该结果中，$\frac{1}{2} (\boldsymbol{x}^*)^T A \boldsymbol{x}^*$ 是一个完全无关的常数项，这意味着数学意义上

$$\text{Find} \min_{\boldsymbol{x}} \varphi(\boldsymbol{x}) \iff \text{Find} \min_{\boldsymbol{e}} \frac{1}{2} \boldsymbol{e}^T A \boldsymbol{e}$$

对于方阵 $A$，定义 $A$ - 范数

$$\|\boldsymbol{x}\|_A = \sqrt{\boldsymbol{x}^T A \boldsymbol{x}}$$

并可以自然诱导出内积

$$\langle \boldsymbol x, \boldsymbol y \rangle_A = \frac{1}{4} (\|\boldsymbol x + \boldsymbol y\|_A^2 - \|\boldsymbol x - \boldsymbol y\|_A^2) = \boldsymbol x^T A \boldsymbol y$$

那么我们可以说：真正应该关注的是 $\|\boldsymbol e\|_A$ 的大小

从当前点 $\boldsymbol{x}^k$ 出发，沿着给定的搜索方向 $\boldsymbol{p}^k$ 走一个确定的步长 $\alpha_k$，得到下一个点：

$$\boldsymbol{x}^{k+1} = \boldsymbol{x}^k + \alpha_k \boldsymbol{p}^k$$

残差 $\boldsymbol{r}$ 的定义得到

$$\begin{aligned} \boldsymbol{r}^{k+1} & = \boldsymbol b - A \boldsymbol{x}^{k+1} \\ & = \boldsymbol b - A (\boldsymbol{x}^k + \alpha_k \boldsymbol{p}^k) \\ & = \boldsymbol b - A \boldsymbol{x}^k - \alpha_k A \boldsymbol{p}^k \\ & = \boldsymbol r^k - \alpha_k A \boldsymbol p^k \end{aligned}$$

可以看出，沿着 $\boldsymbol{p}^k$ 走完一步后，系统中产生了一个全新的增量向量：$-\alpha_k A\boldsymbol{p}^k$。
按照 Gram-Schmidt 正交化的思路，我们把这个 $-\alpha_k A\boldsymbol{p}^k$ 在过去所有残差方向上的投影，一个一个地减掉，剥离出纯净的新维度

$$\boldsymbol r^{k+1} = -\alpha_k A \boldsymbol p^k - \sum_{i=0}^{k} \frac{\langle \boldsymbol r^i, -\alpha_k A \boldsymbol p^k \rangle}{\langle \boldsymbol r^i, \boldsymbol r^i \rangle} \boldsymbol r^i$$

如果我们已经计算出了此时的负梯度，也就是残差方向 $\boldsymbol r^k$
如果不加任何处理，直接令 $\boldsymbol{p}^k = \boldsymbol{r}^k$ 并顺着走下去，那就是最速下降法，上面分析已经说过了这会导致问题
我们对 $\{\boldsymbol r^i\}$ 执行 $A-$ 范数下的 Gram-Schmidt 正交化，得到新的搜索方向 $\boldsymbol p^k$，保证它与之前的残差方向 $\{\boldsymbol r^i\}$ 都是 $A-$ 正交的，即

$$\boldsymbol p^k = \boldsymbol r^k - \sum_{i=0}^{k-1} \frac{\langle \boldsymbol r^k, \boldsymbol p^i \rangle_A}{\langle \boldsymbol p^i, \boldsymbol p^i \rangle_A} \boldsymbol p^i$$

该算法计算量非常复杂，内存占用也极高

共轭梯度法的高明之处就在于此

由于我们在由 $\boldsymbol{p}^0, \dots, \boldsymbol{p}^{k-1}$ 张成的子空间中找到了极小值，所以当前的负梯度（残差 $\boldsymbol{r}^k$）必然垂直于这个子空间中的所有向量

$$(\boldsymbol{r}^k)^T \boldsymbol{r}^i = 0 \quad (\forall i < k)$$

根据残差更新公式 $\boldsymbol{r}^{k+1} = \boldsymbol{r}^k - \alpha_k A \boldsymbol{p}^k$ 可以得到

$$A \boldsymbol{p}^i = \frac{1}{\alpha_i} (\boldsymbol{r}^i - \boldsymbol{r}^{i+1})$$

代入到搜索方向的求和项的分子中

$$(\boldsymbol{r}^k)^T A \boldsymbol{p}^i = (\boldsymbol{r}^k)^T \left[ \frac{1}{\alpha_i} (\boldsymbol{r}^i - \boldsymbol{r}^{i+1}) \right] = \frac{1}{\alpha_i} \left[ (\boldsymbol{r}^k)^T \boldsymbol{r}^i - (\boldsymbol{r}^k)^T \boldsymbol{r}^{i+1} \right]$$

考察求和号中的历史索引 $0 \leq i < k$

• 当 $i < k - 1$ 时，$(\boldsymbol{r}^k)^T \boldsymbol{r}^i = 0$，因此分子为 $0$
• 当 $i = k - 1$ 时，$(\boldsymbol{r}^k)^T \boldsymbol{r}^{k-1} = 0$，因此剩余

$$(\boldsymbol{r}^k)^T A \boldsymbol{p}^{k-1} = - \frac{1}{\alpha_{k-1}} (\boldsymbol{r}^k)^T \boldsymbol{r}^k$$

这意味着巨大的求和项中，只有最后一个成分是非零的，因此搜索方向的更新公式可以简化为

$$\boldsymbol p^k = \boldsymbol r^k - \beta_{k-1} \boldsymbol p^{k-1}$$

接下来推导系数 $\alpha_k$ 和 $\beta_k$ 的计算公式，根据条件

$$(\boldsymbol{r}^{k+1})^T \boldsymbol{p}^k = 0$$

代入 $\boldsymbol{r}^{k+1} = \boldsymbol{r}^k - \alpha_k A \boldsymbol{p}^k$：

$$(\boldsymbol{r}^k - \alpha_k A \boldsymbol{p}^k)^T \boldsymbol{p}^k = 0 \implies \alpha_k = \frac{(\boldsymbol{r}^k)^T \boldsymbol{p}^k}{(\boldsymbol{p}^k)^T A \boldsymbol{p}^k}$$

代入 $\boldsymbol{p}^k = \boldsymbol{r}^k - \beta_{k-1} \boldsymbol{p}^{k-1}$：

$$(\boldsymbol{r}^k)^T \boldsymbol{p}^k = (\boldsymbol{r}^k)^T (\boldsymbol{r}^k - \beta_{k-1} \boldsymbol{p}^{k-1}) = (\boldsymbol{r}^k)^T \boldsymbol{r}^k + \underbrace{-\beta_{k-1} (\boldsymbol{r}^k)^T \boldsymbol{p}^{k-1}}_{=0}$$

因此

$$\alpha_k = \frac{(\boldsymbol{r}^k)^T \boldsymbol{r}^k}{(\boldsymbol{p}^k)^T A \boldsymbol{p}^k}$$

接下来看 $\beta_k$，代入求和项中的系数得到

$$\beta_k = \frac{(\boldsymbol{r}^{k+1})^T A \boldsymbol{p}^k}{(\boldsymbol{p}^k)^T A \boldsymbol{p}^k}$$

代入上述推导中剩余项 $(\boldsymbol{r}^k)^T A \boldsymbol{p}^{k-1}$ 的结果

$$\beta_k = \frac{\frac{1}{\alpha_k} (\boldsymbol{r}^{k+1})^T \boldsymbol{r}^{k+1}}{(\boldsymbol{p}^k)^T A \boldsymbol{p}^k} = \frac{(\boldsymbol{r}^{k+1})^T \boldsymbol{r}^{k+1}}{(\boldsymbol{r}^k)^T \boldsymbol{r}^k}$$

综上，共轭梯度法的算法如下

• 初始化解 $\boldsymbol x^0$，设定合适的极小值值 $\varepsilon$
• 计算初始残差 $\boldsymbol r^0 = \boldsymbol b - A \boldsymbol x^0$
• 设定初始搜索方向 $\boldsymbol p^0 = \boldsymbol r^0$
• 对 $k = 0, 1, 2, \dots$ 开始迭代
  • 计算步长 $\alpha_k = \dfrac{(\boldsymbol{r}^k)^T \boldsymbol{r}^k}{(\boldsymbol{p}^k)^T A \boldsymbol{p}^k}$
  • 更新解 $\boldsymbol x^{k+1} = \boldsymbol x^k + \alpha_k \boldsymbol p^k$
  • 更新残差 $\boldsymbol r^{k+1} = \boldsymbol r^k - \alpha_k A \boldsymbol p^k$
    • 如果 $\|\boldsymbol r^{k+1}\| < \varepsilon$，则认为已经收敛，停止迭代
  • 为了更新搜索方向，计算系数 $\beta_k = \dfrac{(\boldsymbol{r}^{k+1})^T \boldsymbol{r}^{k+1}}{(\boldsymbol{r}^k)^T \boldsymbol{r}^k}$
  • 更新搜索方向 $\boldsymbol p^{k+1} = \boldsymbol r^{k+1} - \beta_k \boldsymbol p^k$