高斯列主消元法

1. 高斯列主消元法的基本思想

列主元素消去法是为控制舍入误差而提出来的一种算法。列主元素消去法计算基本上能控制舍入误差的影响，其基本思想是：在进行第 $k(k=1,2,\dots ,n-1)$步消元时，从第$k$列的 $a_{kk}$及其以下的各元素中选取绝对值最大的元素，然后通过行变换将它交换到主元素$a_{kk}$的位置上，再进行消元。
误差来源:由于计算机表示浮点数有位数限制(IEEE 754 标准中，单精度浮点数(float)是16位，双精度浮点数(double)为32位)，在高斯消元法中，计算行乘数$m_{ij}=\frac{a_{ij}}{a_{ii}}$，如果$a_{ii}$很小，那么行乘数就可能溢出。例如:
$$\begin{array}{rl}(A,b)& =\left(\begin{array}{llll}10^{-8}& 2& 3& 1\\ -1& 3.712& 4.623& 2\\ -2& 1.072& 5.643& 3\end{array}\right),\end{array}$$
如果直接使用$10^{-8}$作为$a_{ii}$，假设计算机使用的是8位十进制尾数，那么在计算$a_{22}$时，$$a_{22}=a_{22}-m_{21}{a}_{12}=a_{22}-\frac{a_{21}}{a_{11}}{a}_{12}=3.712-\frac{-1}{10^{-8}}\times 2=0.0000000\times 10^9+0.2\times 10^9=0.2\times 10^9,$$
由于减数在规格化以后为$10^9$，导致3.712在进行浮点运算的对阶操作时阶码为9，尾数右移为机器0，这样在进行后续的消元后导致方程的解误差过大甚至是错误:经过两步消元后，有:
$$(A^{(3)},b^{(3)})=\left(\begin{array}{cccc}10^{-8}& 2& 3& 1\\ & 0.2\times 10^9& 0.3\times 10^9& 0.1\times 10^9\\ & & 0& 0\end{array}\right).$$
该方程在使用很小的数作为主元时导致方程没有唯一解，实际上该方程是有解的($detA\ne 0$)。因此需要在进行高斯消元法时选择该列的绝对值最大元素作为主元，然后在进行消元

2. 高斯列主消元法的基本步骤

其中$eps$是一个阈值，当某一列的所有列主元都小于$eps$时，该方程的解误差会很大，因此直接放弃使用计算机求解

3. 高斯列主消元法的程序设计

python"># 高斯列主消元法# 高斯列主消元法
import numpy as np
eps = 1e-20
def gaussian_elimination_with_pivoting(A, b):n = A.shape[0]# 增广矩阵Ab = np.hstack([A, b.reshape(-1, 1)])for i in range(n):# 找到列主元的行号max_row_index = np.argmax(np.abs(Ab[i:n, i])) + i# 交换当前行和列主元所在的行if i != max_row_index:Ab[[i, max_row_index]] = Ab[[max_row_index, i]]if abs(Ab[max_row_index][i]) < eps:print(f'列{i}的最大主元为{Ab[max_row_index][i]}')print("数值不稳定/无解")return None# 消元for j in range(i + 1, n):factor = Ab[j, i] / Ab[i, i]Ab[j, i:] = Ab[j, i:] - factor * Ab[i, i:]# 回带求答案x = np.zeros(n)for i in range(n - 1, -1, -1):x[i] = (Ab[i, -1] - Ab[i, i + 1: n] @ x[i + 1:]) / Ab[i, i]return x

4.计算复杂度

根据算法中的消元过程中,第$k$步运算做乘法$（n-k）（n-k+1)$次，做除法$n-k$ 次,完成$n-1$ 步消元共需做乘除法的总次数为$$\sum _{k=1}^{n-1}(n-k)(n-k+2)=\frac{n^3}{3}+\frac{n^2}{2}-\frac{5n}{6},$$

在回带带计算解时，乘除法计算次数为$$\sum _{i=1}^n(n-i+1)=\frac{n^2}{2}+\frac{n}{2}$$
因此全部乘除法运算次数为$$MD=\frac{n^3}{3}+n^2-\frac{n}{3}$$