课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第二讲，这一讲主要讨论了线性方程组的求解以及LU分解。

Chapter 2 线性系统和LU分解

线性系统的可求解性

线性方程组的形式如下：

解的情形一共以下三种：

  1. 唯一解：

  2. 无解：

  3. 无穷多组解：

这一讲我们只讨论$A$是可逆方阵的情形。注意求解线性方程组并不需要计算$A^{-1}$，这是因为$A^{-1}$的计算成本太大。

消元法

课程中老师演示了消元法求解线性方程组的过程，其中只涉及到三种操作：

  1. 对调矩阵中某两行。
  2. 用非零常数乘以矩阵的某一行。
  3. 将矩阵的某一个行乘以常数$c$后加到另一行。

这三种操作对应了三种矩阵：

置换

假设$e_i $是$m$维向量，且$e_i$第$i$个元素为$1$，其余元素都为$0$，定义如下置换映射（该映射为双射）：

那么第一种操作对应的矩阵为：

伸缩

第二种操作对应的矩阵为：

消除

对第$k$行乘以$c$加到第$l$行对应的矩阵为：

注意以上三种矩阵都可可逆。

消元法可以用左乘这些矩阵来代替，用$E$表示上述三种矩阵，所以消元的过程可表示为

最后我们的结果是

从而

因此

高斯消元法

将消元法稍作总结即可得到高斯消元法，分为两个步骤，分别是Forward Substitution和Back Substitution。

Forward Substitution

首先得到增广矩阵：

利用伸缩操作将第一行第一列的元素变为$1$（这里假设这个步骤可行），得到

用第三类变化将第一列的其余元素都变为$0​$可得

对其余每列重复该操作，最终得到

Forward Substitution最终结果得到的是上三角矩阵，接下的步骤是Back Substitution。

Back Substitution

Back Substitution的步骤是将上三角矩阵变成单位阵，首先对最后一列进行操作：

重复此步骤最终得到

最后一列即为方程组的解，整个过程的算法如下：

接着分析算法时间复杂度，注意每一步操作都是对一整行进行，所以每一步操作的时间复杂度为$O(n)$，因为一共有$O(n^2)$个步骤，所以算法的时间复杂度为$O(n^3)​$。

分析

高斯消元法会带来一定的问题，首先考虑下例

显然之前的算法无法处理上述情形。接着考虑下例

利用之前的算法，我们会得到

因此我们会得到$\frac 1 \epsilon ​$的项，如果$\epsilon ​$非常小，那么$\frac 1 \epsilon ​$会非常大，这就会产生误差。

解决上述问题的方法是选主元：即交换行和列以避免除以很小的数或零，选主元有两个主要策略：

1.部分主元：交换行，使得当前列绝对值最大的元素出现在对角线位置。
2.全局主元：交换行列，使得绝对值最大的元素出现在对角线位置。

如果我们有一系列线性方程组：

因为对$A$的操作相同，所以可以将上述方程组合并为

这样可以减少很多重复操作。

$LU$分解

首先考虑对上三角矩阵进行高斯消元法：

由之前的算法可知，一共要循环$n$次，注意到这里为上三角矩阵，所以每一轮迭代的时间复杂度实际上为$O(n)$（因为不需考虑$0$元素），所以总共的时间复杂度为$O(n^2)$，同理可知，对下三角矩阵使用高斯消元法的时间复杂度也为$O(n^2)$。这个重要的结果是$LU$分解的出发点。

回顾高斯消元法，我们知道Forward Substitution的过程实际上是左乘一些列初等矩阵，得到

其中$U=E_k…E_1 A​$是上三角矩阵，那么

其中

注意$E_i​$为下三角矩阵，所以$L​$也为下三角矩阵，因此$A=LU​$将$A​$分解为下三角矩阵和上三角矩阵的乘积，所以原方程组即为

结合之前的讨论可知，我们可以把原问题拆成两个子问题：

因为$L$是下三角矩阵，$U$是上三角矩阵，所以上述算法的时间复杂度为$O(n^2)$。

实现$LU$

作为一个小技巧，我们可以把$L,U​$存储到一个矩阵中来减少空间开销，原理很简单，构造对角矩阵$S​$，使得$LS​$的对角元为$1​$，那么

因为$S​$为对角矩阵，所以$LS,S^{-1}U​$仍然分别为下三角矩阵和上三角矩阵，最终存储结果如下：

在不交换行和列来选择主元的条件下，$LU$分解伪代码如下：

// Takes as input an n-by -n matrix A[i,j]
// Edits A in place to obtain the compact LU factorization described above
for pivot from 1 to n {
	pivotValue = A[pivot , pivot ]; // Bad assumption that this is nonzero !
    for eliminationRow from ( pivot +1) to n { // Eliminate values beneath the pivot
        // How much to scale the pivot row to eliminate the value in the current
        // row; notice we 're not scaling the pivot row to 1, so we divide by the
        // pivot value
        scale = A[ eliminationRow , pivot ] / pivotValue ;
        // Since we / subtract / scale times the pivot row from the current row
        // during Gaussian elimination , we / add/ it in the inverse operation
        // stored in L
        A[ eliminationRow , pivot ] = scale ;
        // Now , as in Gaussian elimination , perform the row operation on the rest
        // of the row: this will become U
        for eliminationCol from ( pivot +1) to n
        	A[ eliminationRow , eliminationCol ] -= A[pivot , eliminationCol ] * scale ;
     }
}

代码解释：

因为

所以$L$中存储的是$E_i^{-1}$，因为$E_{i}$对应的是减法（消去某行的元素），所以$E_i^{-1}$对应的是加法，这一步对应了

scale = A[ eliminationRow , pivot ] / pivotValue ;
A[ eliminationRow , pivot ] = scale ;

然后

对应的代码为

for eliminationCol from ( pivot +1) to n
        A[ eliminationRow , eliminationCol ] -= A[pivot , eliminationCol ] * scale ;