---
url: /computational-physics/lesson-3/index.md
---
## Part 2 高斯消元法

我们希望求解的线性系统可以统一表达为，

$$
A \cdot \mathbf{x} = \mathbf{b} ;,
\tag{a}
$$

其中 $A \in \mathbb{C}^{n \times n}$ 是一个一般的 $n \times n$ 的复矩阵，$\mathbf{b} \in \mathbb{C}^n$ 是 $\mathbb{C}^n$ 中已知的复矢量，而 $\mathbf{x} \in \mathbb{C}^n$ 则是待求的矢量。我们一般将假设 $A$ 是不奇异的，否则方程一般无解。

数值上求解这类方程其实与我们在初中时学习的步骤基本一致。这实际上是一个古老的方法，称为 **高斯消元法 (Gauss Elimination Method, GEM)**。我们首先写出这些方程，然后对它们进行一系列的所谓的初等变换——即将某个方程乘以一定的数然后将其与另一个方程相加减——其目的是消去某个变量前面的系数。通过这样的变换，将原先的方程化为如下形式的方程：

$$
U \cdot \mathbf{x} = \mathbf{c} ;, \quad U = \begin{bmatrix} u\_{11} & u\_{12} & \cdots & u\_{1n} \ 0 & u\_{22} & \cdots & u\_{2n} \ \vdots & & \ddots & \vdots \ 0 & \cdots & 0 & u\_{nn} \end{bmatrix} ;,
\tag{b}
$$

其中 $U$ 是一个上三角矩阵。如果矩阵 $A$ 是非奇异的，那么上三角矩阵 $U$ 一定也是如此，这意味着所有的对角元 $u\_{ii} \ne 0$。因此这个上三角系统可以运用所谓 **反代 (back-substitution)** 的方法来求解，即首先从最后一个方程解出 $x\_n$，反代到倒数第二个方程解出 $x\_{n-1}$；然后将这两者反代到倒数第三个方程解出 $x\_{n-2}$，等等。用公式来表达就是，

$$
x\_i = \frac{c\_i - \sum\_{k=i+1}^n u\_{ik}x\_k}{u\_{ii}} ;, \quad i = n, n-1, \cdots, 1 ;.
\tag{c}
$$

为了要将线性系统 (a) 化为上三角的形式 (b)，我们将原先的矩阵 $A$ 和 $\mathbf{b}$ 合并为一个 $n \times (n+1)$ 的矩阵，它称为原来线性系统 (a) 的 **增广矩阵 (augmented matrix)**，我们将用 $(A, \mathbf{b})$ 来标记它，

$$
(A, \mathbf{b}) = \begin{bmatrix} a\_{11} & \cdots & a\_{1n} & b\_1 \ \vdots & & \vdots & \vdots \ a\_{n1} & \cdots & a\_{nn} & b\_n \end{bmatrix} ;,
$$

我们随后的线性变换都是直接对增广矩阵 $(A, \mathbf{b})$ 来进行操作的。最终的目的是将其变换为上三角的形式 (b)。

我们首先去寻找一个 $a\_{r1} \ne 0$。对所有的 $1 \le r \le n$ 来说，这样的矩阵元总是存在的，否则矩阵 $A$ 将是奇异的，与我们的初始假设矛盾。如果我们找到了这样的一个矩阵元，我们首先可以做的是将第 1 行与第 $r$ 行进行互换。这可以通过一个交换矩阵来实现。假设通过交换行 1 和行 $r$ 将 $(A, \mathbf{b})$ 变为 $(\bar{A}, \bar{\mathbf{b}})$，那么这个变换可以表达为，

$$
(\bar{A}, \bar{\mathbf{b}}) = P^{(1r)} \cdot (A, \mathbf{b}) ;,
$$

其中 $n \times n$ 的 **置换矩阵 (permutation matrix)** $P^{(1r)}$ 的矩阵元为，

$$
\begin{cases}
\[P^{(1r)}]*{1r} = \[P^{(1r)}]*{r1} = 1 ;, \quad \[P^{(1r)}]*{11} = \[P^{(1r)}]*{rr} = 0 ;, \\\\
\[P^{(1r)}]*{ij} = \delta*{ij} ;, \quad (ij) \ne (1r) \ne (r1) ;.
\end{cases}
$$

也就是说，除了第 1 和第 $r$ 行、第 1 和第 $r$ 列这四个矩阵元之外，其他的矩阵元与单位矩阵的相应矩阵元完全一样；\\

而在这四个矩阵元的位置，它看起来就像 $\begin{pmatrix} 0 & 1 \ 1 & 0 \end{pmatrix}$ 一样。

大家很容易验证，这个矩阵是一个非奇异的矩阵并且它的逆矩阵就是它本身。当我们用它左乘上矩阵 $(A, \mathbf{b})$ 时，其作用是交换该矩阵的第 1 行和第 $r$ 行。经过这个变换，$(A, \mathbf{b})$ 变换到了 $(\bar{A}, \bar{\mathbf{b}})$，用数学表达式来写就是：$(\bar{A}, \bar{\mathbf{b}}) = P^{(1r)} \cdot (A, \mathbf{b})$。

另外一类初等变换的矩阵的相貌这样的，

$$
G^{(1)} = \begin{bmatrix}
1 & 0 & \cdots & 0 \\
-l\_{21} & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
-l\_{n1} & 0 & \cdots & 1
\end{bmatrix} ;,
$$

它是一个下三角矩阵，并且仅仅在一列（具体到这个例子是第一列）与单位矩阵不同，其他地方与单位矩阵相同。这样的矩阵在数学上称为 **Frobenius 矩阵**。它的作用是，只要我们适当地选取 $l\_{r1}$ 的数值，就可以将增广矩阵 $(\bar{A}, \bar{\mathbf{b}})$ 中第一列中除了第一个元素 $\bar{a}\_{11}$ 之外的其他元素统统变换为零。

这个具体的选择是 $l\_{r1} = \bar{a}*{r1}/\bar{a}*{11}, \quad r=2, \cdots, n$，其中我们假定了 $\bar{a}*{11} \ne 0$。容易证明，这个矩阵也是非奇异的，事实上它的逆矩阵也是一个 Frobenius 矩阵，只不过其中的 $-l*{r1}$ 都要换成 $+l\_{r1}$。

因此，经过上述的两个变换，我们总可以将原先的最一般的增广矩阵 $(A, \mathbf{b})$ 变换为如下的形式，

$$
(A', \mathbf{b}') = \begin{pmatrix}
a'*{11} & a'*{12} & \cdots & a'*{1n} & b'*1 \\
0 & a'*{22} & \cdots & a'*{2n} & b'*2 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & a'*{n2} & \cdots & a'\_{nn} & b'\_n
\end{pmatrix} ;, \quad (A', \mathbf{b}') = G^{(1)}(\bar{A}, \bar{\mathbf{b}}) = G^{(1)}P^{(1r)}(A, \mathbf{b}) ;, \tag{d}
$$

其中的第一个矩阵 $P^{(1r)}$ 的作用是，首先选择一个不为零的元素 $a\_{r1}$ 并且将它调到第一行，当然如果 $a\_{11}$ 本身就不等于零，这一步原则上也可以省略；第二个矩阵则将矩阵 $(\bar{A}, \bar{\mathbf{b}}) = P^{(1r)}(A, \mathbf{b})$ 的第一列从第二个元素开始往下的矩阵元全部清零。

经过上述两个步骤原先的增广矩阵 $(A, \mathbf{b})$ 被约化了。具体来说，第一个待求的变量 $x\_1$ 仅仅出现在第一个方程之中，后面的 $(n-1)$ 个方程仅仅包含其余的 $(n-1)$ 个变量：$x\_2, \cdots, x\_n$。

显然，我们可以对剩余的这 $(n-1)$ 个变量的增广矩阵——也就是公式 (d) 中由绿色标记出来的矩阵——继续利用上面提及的那两类变换，即将两行对掉（注意，由于这时的对换不牵涉第一行，因此不会破坏整个矩阵 (d) 的结构，也就是说，除了 $a'\_{11}$ 之外，第一列都是零）以及将某行乘以一个数与另一行相加，我们将可以将其变为仅仅包含变量 $(x\_3, \cdots, x\_n)$ 的增广矩阵。

这个过程可以一直迭代地做下去。最终的结果就是我们将增广矩阵变换为了我们希望的上三角的形式 (b) 并进一步利用反代的方法 (c) 获得线性系统的解。

前面提及的寻找的不为零的矩阵元 $\bar{a}*{11} = a*{r1}$ 有个专门的名称，叫做 **支点元 (pivot element)**，或者简称 **支点 (pivot)**。

而这个过程称为 **支点遴选 (pivot selection)**。虽然任意的不为零的元素都可以作为支点，但是直觉告诉我们它应当尽可能地远离零，因此一个自然的选择是
$$
\bar{a}*{11} = \max\_r |a*{r1}| ;.
$$

的确，数学上可以证明这样的选择造成的误差比起随便胡乱选择的支点来说是比较小的。这个选择一般称为 **部分支点遴选 (partial pivoting)**。与之相比，我们还可以进行所谓的 **完全支点遴选 (complete pivoting)**。在完全支点遴选中，我们不仅仅局限于第一列，而是在所有矩阵元中挑选模最大的作为支点，即，

$$
\bar{a}*{11} = \max*{r,s} |a\_{rs}| ;.
$$

随后我们将矩阵的第 1 行与第 $r$ 行，第 1 列与第 $s$ 列对换（这相当于将原先的解 $\mathbf{x}$ 的第 1 个分量与第 $s$ 个分量进行了一次对换）。

> 如果愿意，这两个操作也可以用矩阵来表达：令 $P^{(1r)} \cdot (A, \mathbf{b}) = (A'', \mathbf{b}'')$，那么 $(\bar{A}, \bar{\mathbf{b}}) = (A'' \cdot P^{(1s)}, \bar{\mathbf{b}})$。同时记住，$x\_1$ 与 $x\_s$ 进行了对换。

这样一来，原先的增广矩阵 $(A, \mathbf{b})$ 就变为了一个新的增广矩阵 $(\bar{A}, \bar{\mathbf{b}})$。然后可以进一步利用 $G^{(1)}$ 变换将其化为 (d) 的形式。我们可以将其简记为，
$$
(A', \mathbf{b}') = \begin{bmatrix} a'\_{11} & | & \mathbf{a}'^T & | & b'\_1 \ \hline 0 & | & \bar{A} & | & \bar{\mathbf{b}} \end{bmatrix} ;,
$$

其中 $\mathbf{a}'$ 和 $\bar{\mathbf{b}}$ 都是具有 $(n-1)$ 个分量的矢量， $\bar{A}$ 则是一个 $(n-1) \times (n-1)$ 的复矩阵。下面的步骤就是对 $(n-1)$ 阶的增广矩阵 $(\bar{A}, \bar{\mathbf{b}})$ 继续重复上面的步骤就可以将其最终化为上三角形式。

在第一步的支点遴选的过程中的置换矩阵 $P^{(1r)}$ 的上标 $r$ 其实并不是必须的，它只是临时从各个行中计算出来的一个具有最大模的矩阵元的行指标而已。因此为了下面描述的方便，我们将这一步的置换矩阵记为 $P^{(1)}$，即隐去其上标 $r$。

并且我们将约化后的矩阵记为 $(A^{(1)}, \mathbf{b}^{(1)})$。原始的增广矩阵则给它一个上标 0，即 $(A^{(0)}, \mathbf{b}^{(0)}) \equiv (A, \mathbf{b})$。这样一来第一步的约化可以表达为：
$$
(A^{(1)}, \mathbf{b}^{(1)}) = G^{(1)}P^{(1)}(A^{(0)}, \mathbf{b}^{(0)}) ;.
$$

利用这个记号，我们可以比较方便地将整个约化过程表述为，

$$
\begin{cases}
(A^{(j)}, \mathbf{b}^{(j)}) = G^{(j)}P^{(j)}(A^{(j-1)}, \mathbf{b}^{(j-1)}) ;, \quad j=1,2, \cdots, (n-1) ;. \\\\
(A^{(n-1)}, \mathbf{b}^{(n-1)}) \equiv (U, \mathbf{c})
\end{cases}
$$

或者更为明确地写出，

$$
(U, \mathbf{c}) = G^{(n-1)}P^{(n-1)}G^{(n-2)}P^{(n-2)} \cdots G^{(1)}P^{(1)}(A, \mathbf{b}) ;.
$$

这个过程就是高斯消元法的完整过程。经过这个约化，原先的线性系统被成功约化为一个上三角线性系统。

高斯消元法并不是总能够一直进行下去的。在消元的过程之中，有可能会遇到支点为零的情形。这时候除非我们进行支点的重新选择，否则算法没有办法继续。因此，可以直接进行消元处理的矩阵比起一般的非奇异矩阵来说，需要更多的条件。

下面我们列出进行部分支点遴选的高斯消元法的算法基本步骤。

> \[!important]
>
> **Algorithm 1 部分支点遴选的高斯消元法**
>
> > **Require**:
> >
> > 设 $A \in \mathbb{C}^{n \times n}$ 为一方阵。我们通过消元法化为上三角矩阵，然后利用反代法求解之。
> >
> > 【计算量】：大约 $2n^3/3$.
> >
> > 1. **for** $i = 1, \cdots, n$ **do**
> >
> > 2. 寻找一个不为零的支点元，令：
> >    $$
> >    \bar{a}*{ii} = \max*{i \le r \le n} |a\_{ri}| ;.
> >    $$
> >
> > 3. 将第 $i$ 行与相应的第 $r$ 行互换，这等价于左乘上一个置换矩阵 $P^{(i)}$；
> >
> > 4. 利用一般的 Frobenius 矩阵 (e) $G^{(i)}$ 左乘上面得到的矩阵，从而将 $\bar{a}\_{ii}$ 以下的矩阵元消为零；
> >
> > 5. **end for**
> >
> > 6. 这样一来矩阵已经化为上三角矩阵，可以进而利用反代法给出最后的解。

利用 GEM 进行消元并求解线性方程的计算量可以按照下列方法进行计算。消元的过程需要的计算量为 $2(n-1)n(n+1)/3 + n(n-1)$，随后求解两个三角线性系统还需要 $n^2$ 的计算量。因此整个过程需要计算量大约为 $(2n^3/3 + 2n^2)$，其领头阶的计算量为 $2n^3/3$。这对于量级为 $n \sim O(100)$ 的矩阵是没有任何问题的。

关于线性系统的稳定性与误差分析，我们需要引入矩阵的 **条件数** 的概念。给定一个非奇异的方阵 $A \in \mathbb{C}^{n \times n}$，它的条件数定义为，

$$
K(A) = ||A|| \cdot ||A^{-1}|| ;.
$$

其中的模 $||\cdot||$ 可以是任意良好定义的模。例如，如果是我们的前面定义的 $p-$ 模，我们会将相应矩阵的模记为 $K\_p(A)$。一个矩阵的条件数一般来说依赖于模的选取。不过奇异矩阵的条件数总是趋于无穷大，而一个接近奇异的矩阵则具有非常大的条件数。虽然任何的矩阵模都可以用于条件数的定义，但是通常我们总是采取服从乘法的矩阵模。对于这类矩阵模我们有，

$$
1 = ||AA^{-1}|| \le ||A|| \cdot ||A^{-1}|| = K(A) ;.
$$

因此由服从乘法的矩阵模所定义的条件数总是大于 1 的。比较常用的是采用 $p$-模定义的矩阵模。这时候矩阵的条件数记为 $K\_p(A)$。其中最为常用的是欧氏模定义的 $K\_2(A)$。我们有，

$$
K\_2(A) = \frac{\sigma\_1(A)}{\sigma\_n(A)} ;,
$$

其中 $\sigma\_1(A)$ 和 $\sigma\_n(A)$ 分别是矩阵 $A$ 的最大和最小的 **奇异值** 。如果局限于对称正定矩阵来说，我们有，

$$
K\_2(A) = \frac{\lambda\_{\max}}{\lambda\_{\min}} ;,
$$

其中 $\lambda\_{\max}$ 和 $\lambda\_{\min}$ 分别是矩阵 $A$ 的最大、最小本征值。

一般来说，如果考虑舍入误差等因素，我们的初始数据——这包括矩阵 $A$ 和方程的右边矢量 $\mathbf{b}$——都是具有误差的。因此，我们真正求解出来解可以写成下面方程的形式，

$$
(A + \delta A)(\mathbf{x} + \delta \mathbf{x}) = (\mathbf{b} + \delta \mathbf{b}) ;.
$$

其中 $\delta A$, $\delta \mathbf{b}$ 衡量了初始数据的误差，而 $\delta \mathbf{x}$ 则显示了真实求出的解 $(\mathbf{x} + \delta \mathbf{x})$ 对原始解 $\mathbf{x}$ 的偏差。那么我们有如下的结果，

**定理1**

> 对于上面定义的线性方程我们有，
>
> $$
> \frac{||\delta \mathbf{x}||}{||\mathbf{x}||} \le \frac{K(A)}{1 - K(A)||\delta A||/||A||} \left\[ \frac{||\delta \mathbf{b}||}{||\mathbf{b}||} + \frac{||\delta A||}{||A||} \right] ;.
> $$

特别值得注意的是，如果矩阵本身没有误差，即 $\delta A = 0$，那么我们一定有
$$
\frac{||\delta \mathbf{x}||}{||\mathbf{x}||} \le K(A) \frac{||\delta \mathbf{b}||}{||\mathbf{b}||}
$$
因此对于一个大的条件数的矩阵来说，方程右边的相对误差会被放大。正因为如此，我们在求解线性方程组的过程之中，必须尽可能地保持矩阵的条件数不会变大。

这也是为什么幺正变换（相应的实矩阵为正交变换）在矩阵的变换过程之中起了非常重要的意义。可以证明它们是保持一个矩阵的条件数的变换。因此在变换的过程之中不会将矩阵的条件数变差。