CS205A Lecture 2 Linear systems and LU

课程主页: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$分解伪代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// 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}$对应的是加法,这一步对应了

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

然后

对应的代码为

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

本文标题:CS205A Lecture 2 Linear systems and LU

文章作者:Doraemonzzz

发布时间:2019年03月16日 - 22:19:00

最后更新:2019年03月19日 - 00:54:21

原始链接:http://doraemonzzz.com/2019/03/16/CS205A Lecture 2 Linear systems and LU/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。