线性方程组之高斯消元法

在上一篇文章中，我们已经描述了线性方程组和与之相关的矩阵定义。并以一个抛物线的方程为例介绍了求解线性方程组的过程。在本文中，我们介绍线性方程组的一般解法。

1. 算法过程

给定一个由m个方程构成的n元一次方程组，将其写成矩阵的形式如下，并记为$A$： $$ \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n & = & b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n & = & b_2 \\ & \vdots & \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n & = & b_m \end{cases} \Rightarrow \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \end{bmatrix} $$ 我们可以假定矩阵$A$的第一列不全为0，因为若第一列全为零，意味着未知数$x_1$对方程的解没有影响。选择一个第一列不为零的行，与第一行交换。选择的方式随意，可以是第一列中最大的那一行，也可以是第一个非零的那一行，只要对应元素不为零就可以了，我们把这个元素称为这一列的主元。然后将新矩阵中的第一行数乘$1/a_{11}$倍，将矩阵化为第一行第一列的元素为1的形式。对于第$i, 2 \le i \le m$行，我们将之加上新矩阵第一行的$-a_{i1}$倍，得到第一列只有第一行为1其余全为0的矩阵。为了方便，我们记之为矩阵$A^{(1)}$，其中上角标表示对第一行的变换操作。 $$ \begin{bmatrix} 1 & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)} \\ 0 & a_{22}^{(1)} & \cdots & a_{2n}^{(1)} & b_2^{(1)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & a_{m2}^{(1)} & \cdots & a_{mn}^{(1)} & b_m^{(1)} \end{bmatrix} $$

在新矩阵$A^{(1)}$中，我们考察第二列第二行及其右下方的各个系数$a_{ij}^{(1)}$。它们可能全为零，此时化简结束。否则，我们找到第一个不全为0的列，将之记为$j_2$。此时的矩阵A具有如下的形式： $$ \begin{bmatrix} 1 & \cdots & a_{1j_2}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)} \\ 0 & 0 & a_{2j_2}^{(1)} & \cdots & a_{2n}^{(1)} & b_2^{(1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & a_{mj_2}^{(1)} & \cdots & a_{mn}^{(1)} & b_m^{(1)} \end{bmatrix} $$ 类似的，我们取$j_2$列中的一个不为零的元素$a_{j_2i}, 1 < i \le m$所在的行，将之与第2行交换，并数乘$1/a_{j_2i}$倍，将矩阵化成第2行第$j_2$列的元素为1的形式。再对第$i, 1 \le i \le m, i \neq 2$行加上第二行的$-a_{j_22}$倍，得到第$j_2$列只有第二行为1其余全为0的矩阵$A^{(2)}$。 $$ \begin{bmatrix} 1 & \cdots & 0 & \cdots & a_{1n}^{(2)} & b_1^{(2)} \\ 0 & 0 & 1 & \cdots & a_{2n}^{(2)} & b_2^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{mn}^{(2)} & b_m^{(2)} \end{bmatrix} $$

对于剩下的矩阵，我们重复同样的操作，直到第$k$次变换后，剩余系数$a_{ij}$均为0为止。此时得到的矩阵$A^{(k)}$具有如下的形式： $$ \begin{bmatrix} a_{1,1}^{(k)} & \cdots & 0 & \cdots & 0 & a_{1,j_k+1}^{(k)} & \cdots & b_1^{(k)} \\ & & a_{2,j_2}^{(k)} & \cdots & 0 & a_{2,j_k+1}^{(k)} & \cdots & b_2^{(k)} \\ & & & \ddots & \vdots & \vdots & \vdots & \vdots \\ & & & & a_{k,j_k}^{(k)}& a_{k,j_k+1}^{(k)} & \cdots & b_k^{(k)} \\ & & & & & & & \vdots \\ & & & & & & & b_m^{(k)} \end{bmatrix} $$ 其中左下角空白表示对应元素均为0，矩阵中$a_{1,1}, a_{2,j_2}, \cdots a_{k,j_k}$均为1。我们称之为最简形式。下面，我们分情况讨论解。

(1). 如果方程数量等于未知数数量，即$m =n$，的情况。理想情况下，我们可以得到$j_i = i, 1 \le i \le n$形式的矩阵，即： $$ \begin{bmatrix} 1 & \cdots & 0 & b_1^{(n)} \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & 1 & b_n^{(n)} \end{bmatrix} $$ 此时，我们可以直接根据矩阵得出方程组的解为： $$ \begin{cases} x_1 & = & b_1^{(n)} \\ & \vdots & \\ x_n & = & b_n^{(n)} \end{cases} $$

(2). 对于$m=n$的矩阵，有时我们会得到$j_i > i$的情形。这意味着对于矩阵我们经过了$1 \le k < n$次的变换就得到了最简的形式。那么对于矩阵中第$i > k$的行，我们有等式$0x_1 + \cdots + 0x_m = b_i^{(k)}$，若要等式成立，必须有$b_i^{(k)} = 0$，否则方程组无解。

我们把方程中$j_i$列对应的未知数称为非自由变量，其它的未知数称为自由变量。所谓的自由变量的意思是对应的未知数可以取数域$F$中定义的任何数，而总能在非自由变量取特定值时保证方程成立。我们把非自由变量对应的列，即最简形式对应的第$h_u \neq j_i; u = 1, \cdots, v; i = 1, \cdots, k; u + v = n$列，移到方程的右边得到： $$ \begin{cases} x_1 & = & b_1^{(k)} - a_{1,h_1}^{(k)}x_{h_1} - \cdots - a_{1, h_u}^{(k)}x_{h_u} \\ x_{j_2} & = & b_2^{(k)} - a_{2,h_1}^{(k)}x_{h_1} - \cdots - a_{2, h_u}^{(k)}x_{h_u} \\ & \vdots & \\ x_{j_k} & = & b_k^{(k)} - a_{k,h_1}^{(k)}x_{h_1} - \cdots - a_{k, h_u}^{(k)}x_{h_u} \end{cases} $$ 我们可以说方程的解是由自由变量的线性组合再加上一个常数项得到的，我们称之为解空间。

(3). 对于$m < n$的矩阵，我们总可以得到上述的情况(2)。略。

(4). 对于$m > n$的矩阵，我们最多进行n次变换就会得到上述情况(2)。略。

2. 算法实现

上述通过对某一方程进行数乘并与其它方程求和达到消元的目的，进而求得整个方程组的解，这一求解过程就是所谓的高斯消元法。