矩阵与高斯消元：数学算法在计算机领域的应用

一、概述和基本概念

矩阵，类似于在 C++ 中我们看到的二维数组。它有两个维度，行和列。下面是一个典型的矩阵：
$$M=\left[\begin{array}{cccc}1& 2& 3& 4\\ 2& 3& 4& 5\\ 4& 4& 5& 6\\ 10& 11& 12& 13\end{array}\right]$$
我们可以对矩阵进行一些变换，前后矩阵不变，从而实现等量交换。

高斯消元，是解决 $N$ 元一次方程组的高效方法。假设我们有如下的一个方程组：
$$\{\begin{array}{l}{a}_{1,1}{x}_1+a_{1,2}{x}_2+a_{1,3}{x}_3+\cdots +a_{1,n}{x}_n=c_1\\ a_{2,1}{x}_1+a_{2,2}{x}_2+a_{2,3}{x}_3+\cdots +a_{2,n}{x}_n=c_2\\ \dots \\ a_{m,1}{x}_1+a_{m,2}{x}_2+a_{m,3}{x}_3+\cdots +a_{m,n}{x}_n=c_m\end{array}$$
那么，使用 Gauss 消元可以在 $O(n^3)$ 的时间复杂度内解决问题，相较于枚举每个 $x$ 的朴素算法（时间复杂度为 $O(kmn)$ ，其中 $k$ 为 $x$ 的值域大小， $m$ 为方程组个数） ，时间复杂度得到了优化（尤其是 $|k|\le 10^9$ 等值域较大的时候）。

二、矩阵的基本操作

在接下来的叙述里，我们将默认采用 $N$ 个矩阵 $A_1,A_2,\dots ,A_N$ 。

1. 矩阵运算

矩阵的加法：将矩阵的每个数分别相加。具体定义：设进行加法的两个矩阵分别为 $A,B$ ，结果为矩阵 $C$ ，则
$$\forall 1\le i,j\le n,C_{i,j}=A_{i,j}+B_{i,j}$$
举例：
$$\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 2& 5& 6\end{array}\right]+\left[\begin{array}{ccc}1& 3& 5\\ 2& 4& 6\\ 3& 5& 7\end{array}\right]=\left[\begin{array}{ccc}3+1=4& 4+3=7& 5+5=10\\ 1+2=3& 2+4=6& 3+6=9\\ 2+3=5& 5+5=10& 6+7=13\end{array}\right]$$
求解的时间复杂度为 $O(nm)$ ，其中 $n$ 为矩阵列数， $m$ 为矩阵行数。下同，不再赘述。

矩阵的乘法：这里非常有必要提一下。矩阵乘法的结果是目标值所在行的所有数分别乘以目标值所在列的所有数的和。设进行乘法的两个矩阵分别是 $A(n\times r),B(r\times m)$ ，得到了矩阵 $C(n\times m)$ ，则
$$\forall 1\le i,j\le n,C_{i,j}=\sum _{k=1}^r{A}_{i,k}\cdot B_{k,j}$$
较为简便的计算方法如下：
$$C_{i,j}=\sum _{k=1}^r{A}_{i,k}\times \sum _{k=1}^r{B}_{k,j}$$
求解的时间复杂度为 $O(n^3)$ 级别。
举例：
$$\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 2& 5& 6\end{array}\right]\times \left[\begin{array}{ccc}1& 3& 5\\ 2& 4& 6\\ 3& 5& 7\end{array}\right]=\left[\begin{array}{ccc}72& 144& 212\\ 36& 72& 108\\ 78& 156& 238\end{array}\right]$$

矩阵的幂：$A^N$ 表示将同一个矩阵 $A$ 连续乘以它自身 $N$ 次。

根据以上定义，我们可以写出矩阵（matrix）结构体的基本内容，这里假设输入都合法：

struct Matrix
{
	int n, m;
	vector< vector<int> > matrix;
	
	Matrix(int n, int m) : n(n), m(m)
	{
		matrix.clear();
		for (int i = 0; i <= n; ++i)
		{
			vector<int> temp(m + 1);
			temp.clear();
			matrix.push_back(temp);
		}
	}
	
	Matrix(const Matrix& E) : n(E.n), m(E.m), matrix(E.matrix) { };
	
	Matrix operator+(Matrix b) const
	{
		int p = max(n, b.n);
		int q = max(m, b.m);
		Matrix c(p, q);
		for (int i = 1; i <= p; ++i)
		{
			for (int j = 1; j <= q; ++j)
			{
				c.matrix[i][j] = matrix[i][j] + b.matrix[i][j];
			}
		}
		return c;
	}
	
	Matrix operator*(Matrix b)
	{
		Matrix c(n, b.m);
		for (int i = 1; i <= n; ++i)
		{
			for (int j = 1; j <= b.m; ++j)
			{
				for (int k = 1; k <= n; ++k)
				{
				    c.matrix[i][j] += matrix[i][k] * b.matrix[k][j];
				}
			}
		}
		return c;
	}
	
	Matrix operator*(int b)
	{
		Matrix c(n, m);
		for (int i = 1; i <= n; ++i)
		{
			for (int j = 1; j <= m; ++j)
			{
                c.matrix[i][j] = matrix[i][j] * b;
			}
		}
		return c;
	}
};

2. 矩阵运算的性质

大写字母均代表矩阵，小写字母均代表数值。为方便推导，我们假设矩阵的大小均为 KaTeX parse error: Undefined control sequence: \time at position 3: n \̲t̲i̲m̲e̲

（1） 交换律：

i) $A+B=B+A$

证明：令 $C=A+B$ 。
$\forall 1\le i,j\le n,C_{i,j}=A_{i,j}+B_{i,j}$ ，
$\Rightarrow C_{i,j}=B_{i,j}+A_{i,j}$ 。这就证明了 i 。

（2） 结合律：

ii) $(A+B)+C=A+(B+C)$
iii) $(AB)\cdot C=A\cdot (BC)$

证明：令 $D=(A+B)+C$ 。
$\forall 1\le i,j\le n,D_{i,j}=(A_{i,j}+B_{i,j})+C_{i,j}$
$\Rightarrow D_{i,j}=A_{i,j}+(B_{i,j}+C_{i,j})$ 。这就证明了 ii 。
令 $D^{\prime }=(AB)\cdot C$ 。
$\forall 1\le i,j\le n,D_{i,j}^{\prime }=(\sum _{k=1}^n{A}_{i,k}\cdot B_{k,j})C_{i,j}=\sum _{k=1}^n(\sum _{d=1}^n{A}_{i,d}\cdot B_{d,j})\cdot C_{k,j}=\sum _{k=1}^n(\sum _{d=1}^n{B}_{i,d}\cdot C_{d,j})\cdot A_{k,j}=(\sum _{k=1}^n{B}_{i,k}\cdot C_{k,j})A_{i,j}$ 。这就证明了 iii 。

重要推论：因为矩阵的乘法有结合律，所以
$$A^N=\prod _{i=1}^{N}A=\prod _{i=1}^{\frac{N}{2}}A\cdot \prod _{i=\frac{N}{2}+1}^{N}A$$
因此，我们可以使用快速幂进行简便计算，将矩阵幂的复杂度由 $O(n^4)$ 降至 $O(n^3\log n)$ 。
所以，我们可以添加一个 Matrix 类的成员函数：

Matrix pow(int b)
{
	Matrix c(n, m);
	Matrix a(n, m);
	memcpy(a.matrix, matrix, sizeof(matrix));
	for (; b; b >>= 1)
	{
		if (b & 1)
		{
			c = c * a;
		}
		a = a * a;
	}
	return c;
}

（3） 分配率：

iv) $(A+B)C=AC+BC$

证明：令 $D=(A+B)C$ 。
$\forall 1\le i,j\le n,D_{i.j}=\sum _{i=1}^n[(A_{i,k}+B_{i,k})C_{k,j}]=\sum _{i=1}^n(A_{i,k}\cdot C_{k,j}+B_{i,k}\cdot C_{k,j})=\sum _{i=1}^n{A}_{i,k}\cdot C_{k,j}+\sum _{i=1}^n{B}_{i,k}\cdot C_{k,j}$ 。这就证明了 iv 。

3. 特殊矩阵：单位矩阵

单位矩阵是一个正方形，记作 $E$ ，其中从左上角开始到右下角结束的一条对角线上的值是 $1$ ，其余都是 $0$ 。一个典型的单位矩阵如下：
$$E_4=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$
单位矩阵具有如下性质：
v) $\forall \text{ matrix }A,AE=EA=A$

证明：假设 $P=AE$ ， $E$ 的大小为 $n\times n$ 。
显然，有 $\forall 1\le i,j\le n,P_{i,j}=\sum _{k=1}^n{A}_{i,k}\cdot E_{k,j}$ 。
因为在第 $j$ 行中有且仅有第 $j$ 个元素非 $0$ ，即 $E_{j,j}=1$ ，所以 $P_{i,j}=(\sum _{k=1,k\ne j}^n{A}_{i,k}\cdot 0)+A_{i,j}=A_{i,j}$ ，所以 $AE=A$ 。
由矩阵乘法的交换律，知 $AE=EA=A$ 。这就证明了 v 。

4. 矩阵的行初等变换

（1）对换变换：交换矩阵的任意两行，矩阵不变。
（2）倍乘变换：将矩阵的某一行乘以一个常数，矩阵不变。
（3）倍加变换：在对矩阵的某一行进行倍乘变换后加到另一行的数上。
这三点将会被用在 Gauss 消元中。
实现对换变换、倍乘变换、倍加变换的例子如下：
对换变换：
$$\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 2& 5& 6\end{array}\right]⟹\left[\begin{array}{ccc}3& 4& 5\\ 2& 5& 6\\ 1& 2& 3\end{array}\right]$$
倍乘变换：
$$\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 2& 5& 6\end{array}\right]⟹\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 4& 10& 12\end{array}\right]$$
倍加变换：
$$\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 2& 5& 6\end{array}\right]⟹\left[\begin{array}{ccc}3& 4& 5\\ 1& 2& 3\\ 4& 10& 12\end{array}\right]⟹\left[\begin{array}{ccc}3& 4& 5\\ 5& 12& 15\\ 4& 10& 12\end{array}\right]$$

三、 Gauss 消元的算法原理及过程

1. 基本定义

对于以下方程组：
$$\{\begin{array}{l}{a}_{1,1}{x}_1+a_{1,2}{x}_2+a_{1,3}{x}_3+\cdots +a_{1,n}{x}_n=c_1\\ a_{2,1}{x}_1+a_{2,2}{x}_2+a_{2,3}{x}_3+\cdots +a_{2,n}{x}_n=c_2\\ \dots \\ a_{m,1}{x}_1+a_{m,2}{x}_2+a_{m,3}{x}_3+\cdots +a_{m,n}{x}_n=c_m\end{array}$$
系数矩阵：由方程的系数 $a_{i,1},\cdots ,a_{i,n}$ 组成的矩阵，记为 $M$ 。
增广矩阵：在系数矩阵的最右侧添加一列，为 $c_i$ ，记为 $M|c$ 。
Gauss 消元的算法过程在增广矩阵上进行。
例如，对于方程组：
$$\{\begin{array}{l}3x_1+2x_2=7\\ \frac{5}{3}{x}_1+\frac{7}{3}{x}_2=\frac{100}{7}\end{array}$$
它的系数矩阵、增广矩阵分别是：
$$\left[\begin{array}{cc}3& 2\\ \frac{5}{3}& \frac{7}{3}\end{array}\right],\left[\begin{array}{ccc}3& 2& 7\\ \frac{5}{3}& \frac{7}{3}& \frac{100}{7}\end{array}\right]$$

2. 补充知识

主元：在一个矩阵中，对于该矩阵的某一行（指定的），从左到右第一个非 $0$ 的数字就是主元。
行阶梯型矩阵：满足以下 $3$ 个性质的矩阵就是行阶梯型矩阵：

• $\forall 1\le i\le n,i\in \mathbf{Z}$ ，第 $i$ 行的主元永远在第 $i+1$ 行主元的左上方。
• $\forall 1\le i\le n,i\in \mathbf{Z}$ ，第 $i$ 行的主元下方均为 $0$。这一性质其实可以由上面的那一条推出。
• 所有非零行都在所有零行之上。这一性质也可以由第一条推出。

一个典型的行阶梯型矩阵如下，加粗的为主元，可以观察到以上三条性质：
$$\left[\begin{array}{cccccc}3& 10& 2& 3& 5& 1\\ 0& 0& 7& 3& 4& 2\\ 0& 0& 0& 4& 6& 4\\ 0& 0& 0& 0& 0& 10\\ 0& 0& 0& 0& 0& 0\end{array}\right]$$

简化行阶梯型矩阵：在行阶梯型矩阵的基础上满足以下性质：

• 每一个非 $0$ 行的主元一定是 $1$ 。
• 每一个主元 $1$ 是该行唯一的非 $0$ 元素。

Gauss 消元的目标就是把一个矩阵简化为行阶梯型矩阵。

3. Gauss 消元算法的主要过程及其对应在方程组中操作的实际意义

假设我们的增广矩阵 $A$ 是 $n\times m$ 的。
先写出 Gauss 消元算法的主要流程：

1. 枚举列号 $\forall 1\le c\le m\in \mathbf{Z}$ 。
2. 找到主元。假设当前是第 $c$ 列，我们在 $[1,c-1]$ 列这一范围内已经找到了 $w-1$ 个主元，根据行阶梯型矩阵的定义以及性质，下一列的主元一定在上一列的主元的右下方，即一定在第 $c$ 列，第 $[w,n]$ 行范围内。故应该在 $[w,n]$ 行范围内寻找主元。该主元应绝对值最大，找到的主元的绝对值 ∣W∣=max⁡w≤i≤n,w∈Z{∣Ai,c∣}|W|=\displaystyle\max_{w \le i \le n,w\in\bold{Z}}\{|A_{i,c}|\}∣W∣=w≤i≤n,w∈Zmax​{∣Ai,c​∣}
3. 如果我们找到的 $A_{a,b}=0$ ，说明继续操作没有意义，跳过当前列，下一列寻找。否则，将当前找到的主元所在行与第 $w$ 行做对换变换。因此，该操作并没有改变矩阵，只是将矩阵化为另一种形式，方便统计。变换后，主元所在行变成了第 $w$ 行。
4. 对第 $w$ 行做倍乘变换。$\forall 1\le i\le m,A_{w,i}\leftarrow A_{w,i}\times \frac{1}{A_{w,c}}$ ，相当于在方程组中将等式两边的值都除以 $A_{w,c}$ ，即找到的主元。
5. 消元。这是 Gauss 消元算法中最主要的步骤。主要目的是将第 $c$ 列除主元所在位置 $(w,c)$ 之外的所有元素变为 $0$ 。为确保矩阵前后相同，我们对矩阵进行倍加变换。假设我们当前的位置是 $(i,j)$ ，满足 $i\ne w$ 。根据简化行阶梯型矩阵的性质，可知 $j$ 只需要在 $w\le j\le m,j\in \mathbf{Z}$ 的范围内枚举即可，因为 $\forall 1\le j<w,A_{i,j}=0$ 恒成立（简化行阶梯型矩阵性质 2）。倍加变换的主要过程如下：当前枚举到 $(i,j)$ ，主元位置为 $(c,w)$ ，对于第 $w$ 行的所有元素，应该乘以 $\frac{A_{i,c}}{A_{w,c}}$ 。推导过程类似于解方程组的思路，假设我们有一个如下的方程组：
   $$\{\begin{array}{l}4x+2y=6\\ 4x+3y=7\end{array}$$
   那么我们一定会尝试去消去 $4x$ ，这就是加减消元法（使用了初中数学七年级下-二元一次方程组的知识）。解得： $\{\begin{array}{l}x=1\\ y=1\end{array}$ 。如果将它化为增广矩阵，那么将会等于 $\left[\begin{array}{ccc}4& 2& 6\\ 4& 3& 7\end{array}\right]$ ，主元在前面用粗体表示，如果我们尝试让第 $1$ 列的所有元素变成 $0$ ，那么我们会让第一行的所有元素乘以 $\frac{4}{4}=1$ 再减去第二行的所有元素。注意，分母的 $4$ 代表的是 $(1,1)$ 位置的元素，而分子的 $4$ 则代表的是 $(1,2)$ 位置的元素。增广矩阵会变成：$\left[\begin{array}{ccc}4& 2& 6\\ 0& 1& 1\end{array}\right]$ ，此时消元的目的已经达到。
   将这个性质推广到 $A$ ，则可以得出 $\forall (i,j)\in A$ ， $A_{i,j}$ 应该减去 $\frac{A_{i,c}}{A_{w,c}}\cdot \frac{A_{w,c}}{A_{w,j}}=A_{i,c}\cdot A_{w,j}$ 。为方便记忆，可以观察公式，行和列已经用 $w,c$ 给定，将变量 $i,j$ 按照行和列的变化插入即可。
6. 方程组的解即为矩阵的最后一列，$x_1,x_2,\dots$ 与 $A_{n,1},A_{n,2},\dots$ 一一对应。

用一些例子来更好理解 Gauss 消元算法：
$$\begin{gathered} \{\begin{array}{l}4x_1+2x_2+3x_3=15\\ 2x_1+2x_2+4x_3=19\\ 5x_1+6x_2+7x_3=29\end{array} \\ \left[\begin{array}{cccc}4& 2& 3& 15\\ 2& 2& 4& 19\\ 5& 6& 7& 29\end{array}\right]\stackrel{对换变换}{⟹}\left[\begin{array}{cccc}5& 6& 7& 29\\ 2& 2& 4& 19\\ 4& 2& 3& 15\end{array}\right]\stackrel{倍乘变换}{⟹}\left[\begin{array}{cccc}1& \frac{6}{5}& \frac{7}{5}& \frac{29}{5}\\ 2& 2& 4& 19\\ 4& 2& 3& 15\end{array}\right]\stackrel{倍加变换}{⟹} \\ \left[\begin{array}{cccc}1& \frac{6}{5}& \frac{7}{5}& \frac{29}{5}\\ 0& -\frac{2}{5}& -1& \frac{37}{5}\\ 4& 2& 3& 15\end{array}\right]\dots \end{gathered}$$

主要代码（在 Matrix 结构体中）：

void Gauss()
{
	int w = 1;
	for (int c = 1; c <= m; ++c)
	{
		int temp = w;
		for (int i = w + 1; i <= n; ++i)
			if (abs(a[temp][c]) < abs(a[i][c]))
				temp = i;
		if (abs(a[temp][c]) < eps) continue;
		for (int i = c; i <= m; ++i) swap(a[temp][i], a[w][i]);
		for (int i = m; i >= c; --i) a[w][i] /= a[w][c];
		for (int i = 1; i <= n; ++i)
		{
			if (i == w) continue;
			for (int j = m; j >= c; --j)
				a[i][j] -= a[i][c] * a[w][j];
		}
	}
}

注：本文章所有代码未经调试及编译，如有 Bug ，欢迎在评论区指出，必互关 。