【NA】高斯消元法

前言.

• 科学计算与工程应用中，很多问题话归到最终都是线性方程组的形式，例如前面《函数最佳逼近》《最小二乘拟合》中推导出的正规方程，此外还有差分法求解常微分方程、偏微分方程等。
• 在天文学中，Kepler第一定律指出行星的轨道是椭圆，如果想确定一颗小行星绕太阳运行的轨道，可在轨道平面内建立以太阳为原点的直角坐标系，轨道方程可以设为$$a_1{x}^2+2\cdot a_{2}xy+a_3{y}^2+2\cdot a_{4}x+2\cdot a_{5}y+1=0$$
• 如果要确定轨道方程中的参数，可以在5个不同的时间点对小行星进行观测，获取其不同时刻的坐标 $(x_i,y_i)|i=1,2,...,5$，于是轨道参数的求解转化为求解一个线性方程组 $Ax=b.$
• 后面讨论的线性方程组大致可以按照系数矩阵分为两类，一类是低阶稠密矩阵，阶数通常不超过150，非零元素很多；另一类是高阶稀疏矩阵，阶数较高但零元素较多。两类线性方程组也分别对应两种求解策略，对于低阶稠密阵，通常直接进行求解，例如高斯消元法及其诸多改进算法；对于高阶稀疏阵，通常使用迭代过程逐步逼近线性方程组的解。

下面的讨论中线性方程组 $Ax=b$ 的系数矩阵 $A$ 都是可逆矩阵，因此方程组具有唯一解。( $b=\vec{0}$ 时零解是唯一解)

高斯消元法.

• 基本思想就是手工解方程组时的消元策略，不过计算机执行时，首先将系数矩阵通过初等变换转化为上三角矩阵，而后对待求变量从后向前 $x_n\to x_{n-1}\to ...\to x_1$ 回代计算。
• 在计算机中使用高斯消元法求解线性方程组 $Ax=b$ 时只需要存储增广矩阵 $[Ab]$，而后消元过程中只需要对其中的元素直接更新即可。得到上三角矩阵后，可以将解直接存储在增广矩阵最右列。
• 对于下面的增广矩阵$$[Ab]=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& b_1\\ a_{21}& a_{22}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}& b_n\end{array}\right]$$
• 欲将其转化成上三角矩阵，就是要将系数矩阵主对角线下方的元素通过初等变换化为 $0$，以第一列为例，对于主对角线下方的元素 $a_{k1},k=2,3,...,n$ ，只需要做变换 $$a_{k1}=a_{k1}-\frac{a_{k1}}{a_{11}}\cdot a_{11}$$即第 $k$ 行减去 $\frac{a_{k1}}{a_{11}}$ 倍的第 $1$ 行。
• 将其一般化得到，对于第 $t$ 次高斯消元更新，我们更新的目标是第 $t$ 列主对角线以下的元素，即 $a_{kt},k=t+1,t+2,...,n$，令 l k t = − a k t a t t ,   t ∈ { 1 , 2 , . . . , n } , k = t + 1 , t + 2 , . . . , n l_{kt}=-\frac{a_{kt}}{a_{tt}},~t\in\{1,2,...,n\},k=t+1,t+2,...,n lkt​=−att​akt​​, t∈{1,2,...,n},k=t+1,t+2,...,n
• 因此对下面第 $t+1,t+2,...,n$ 行所作变换为 a k = a k + l k t ⋅ a t ,   t ∈ { 1 , 2 , . . . , n } , k = t + 1 , t + 2 , . . . , n a_{k}=a_k+l_{kt}·a_t,~t\in\{1,2,...,n\},k=t+1,t+2,...,n ak​=ak​+lkt​⋅at​, t∈{1,2,...,n},k=t+1,t+2,...,n
• 最终得到的增广矩阵如下$$[Ab]=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& b_1\\ 0& a_{22}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& & ⋮& ⋮\\ 0& 0& \cdots & a_{nn}& b_n\end{array}\right]$$注意此处增广矩阵内元素的值已经发生了变化，据此可得$$x_n=\frac{b_n}{a_{nn}}$$$$x_k=\frac{b_k-{\sum }_{j=k+1}^n{a}_{kj}{x}_j}{a_{kk}},k=n-1,n-2,...,1$$对 $x_k$ 的求解过程即为回代求解。
• 高斯消元法大致分为两步：消元与回代。在消元过程，$l_{kt}$ 的计算需要进行 $n-1$ 次乘法运算，而系数矩阵元素的乘法次数可以这样考虑，第一次消元时，可以直接将第一列主对角线以下的元素赋为0，而需要计算 $(n-1{)}^2$ 个元素的更新值，依次类推，整个系数矩阵的乘法次数为$$\sum _{i=1}^{n-1}{i}^2=\frac{(n-1)\cdot n\cdot (2n-1)}{6}\approx \frac{n^3}{3}$$后面 $b$ 更新所需的乘法次数以及回代过程的乘法次数可以类似计算，最终可以判断出高斯消元法大概的时间复杂度为 $O(n^3)$，相较于仅有理论意义的克莱默法则 $O(n!)$ 无疑是一个极大的提升。
• 高斯消元法的特点是，按顺序对第 $1,2,...,n$ 列在对角线以下的元素进行校园，我们将每次消元时作为 $l_{kt}$ 分母的主对角线元素称为主元素。一般高斯消元法可能会遇到如下两个问题，①如果某个主元素为零，则消元过程无法进行；②如果某个主元素绝对值过小，可能会导致消元过程数值不精确。于是我们尝试对一般高斯消元法做出改进。
• 下图来自Anany P162
  

列主元消去法.

• 为了避免上部分最后提到的两个问题，列主元消去法在每一步进行消元时，先遍历当前第 $t$ 列对角线下方的元素，从中选出绝对值最大的数，假设它处于第 $q$ 行，于是将 $t,q$ 两行交换，就能在一定程度上提升高斯消元法的数值稳定性。
• 时间复杂度和高斯消元法相同为 $O(n^3).$
• 下图来自Anany P163
  

Gauss-Jordan消元法.

• 高斯消元法将系数矩阵 $A$ 转化为上三角矩阵，而高斯若当消元法则将其直接转化为单位矩阵 $E$，并且其中也引入了选取列主元素的策略。
• 在第 $t$ 步消元时，选择第 $t$ 列对角线下方绝对值最大的元素，将其所在行与第 $t$ 行交换。令$$l_{kt}=-\frac{a_{kt}}{a_{tt}},k=1,2,...,t-1,t+1,...,n$$$$l_{tt}=\frac{1}{a_{tt}}$$对每一行都做初等变换$$a_k=a_k+l_{kt}\cdot a_t,k=1,2,...,n$$最后进行对角元素的单位化，即$$a_k=l_{tt}\cdot a_k$$
• 可以发现高斯若当消元法求解线性方程组的乘法次数大约是 $\frac{n^3}{2}$，比一般高斯消元法以及列主元消去法略大，但高斯若当消元法可以很方便用于求出可逆矩阵 $A$ 的逆矩阵 $A^{-1}.$
• 根据线性代数的知识，对分块矩阵 $M=[AE]$ 做初等变换，当分块矩阵中左块 $A$ 被转化为单位矩阵时，右块 $E$ 即被转化为 $A$ 的逆矩阵 $A^{-1}.$

无关紧要.

• 看过Anany这本书的话，在第一段伪代码的后文中，Anany指出此段代码最内层循环效率低下，并且在第二段伪代码中做出了优化。
• 第一段伪代码中效率低下的原因在于，在最内层以 $k$ 为循环变量的循环体中，$A[j,i]\ast A[i,i]$ 是循环不变量，也就是说对于不同的 $k$ 是不会变化的，很自然我们能想到对于它的计算发生了无用的重复计算，因此将其外提到 $k$ 循环外，以临时变量 $temp$ 来标识。
• 编译原理中循环优化会完成这一工作，即循环不变量的外提，现代的高性能编译器在编译优化时，即使我们写出了上面这种无用计算的代码，它也会自动将其转换为第二种使用临时变量的高性能代码。但知道Compiler为我们默默付出了努力，还是有些必要的。
• 《循环优化》