上一篇博文讲了高斯-若尔单的矩阵消元方法,而今天讲的LU分解正是对上一节高斯若尔单消元的拓展,逐步引导你走向线性代数的深处。
引入
我们直接从3×3的矩阵开始。假设我们已经从一个三元方程组中提取出一个3×3的矩阵,接下来我们的任务是对这个矩阵进行消元,使其变为一个上三角矩阵(Upper triangular matrix),后文称为U。
示例:
1 2 1 3 8 1 0 4 2 \begin{matrix} 1 & 2 & 1 \\ 3 & 8 & 1 \\ 0 & 4 & 2 \end{matrix} 130284112
而要使这个矩阵 A 变为上三角矩阵 U,需要乘上两个矩阵——E(3,2)与E(2,1)。(本文章后面都统一采用左乘,即行的线性组合)
E 3 , 2 表示此矩阵乘矩阵 A 能消掉 ( 3 , 2 ) 处的元素 E 2 , 1 表示次矩阵乘矩阵 A 能消掉 ( 2 , 1 ) 处的元素 E_{3,2}表示此矩阵乘矩阵A能消掉(3,2)处的元素\\E_{2,1}表示次矩阵乘矩阵A能消掉(2,1)处的元素 E3,2表示此矩阵乘矩阵A能消掉(3,2)处的元素E2,1表示次矩阵乘矩阵A能消掉(2,1)处的元素
可以得到如下式子:
E 3 , 2 E 2 , 1 A = U E_{3,2}E_{2,1}A=U E3,2E2,1A=U
根据矩阵的逆的性质可得:
E 3 , 2 E 2 , 1 E 2 , 1 − 1 E 3 , 2 − 1 A = E 2 , 1 − 1 E 3 , 2 − 1 U E_{3,2}E_{2,1}E_{2,1}^{-1}E_{3,2}^{-1}A=E_{2,1}^{-1}E_{3,2}^{-1}U E3,2E2,1E2,1−1E3,2−1A=E2,1−1E3,2−1U
即
A = E 2 , 1 − 1 E 3 , 2 − 1 U A=E_{2,1}^{-1}E_{3,2}^{-1}U A=E2,1−1E3,2−1U
这里提示另一点,矩阵AB的逆为 B − 1 A − 1 B^{-1}A^{-1} B−1A−1。
如果硬要比喻的话这个过程就像先穿袜子再穿鞋子,脱的时候就是先脱鞋子再脱袜子。用C++比喻,矩阵的逆就像对栈的特性,先进后出,后进先出
我们将右边两个消元矩阵的逆乘到一起称为L。那么,我们得到了如下两个式子
E A = U A = L U EA=U\\A=LU EA=UA=LU
而LU想要说明的就是下面的式子,A=LU的计算量是小于上面式子AE=U的,因此我们消元时倾向使用LU分解。
为什么选择LU
下面我就举例说明LU分解的合理性。
因为我们后面研究的是E与L,和A和U的关系不大,所以就只写出E矩阵:
[ 1 0 0 0 1 0 0 − 5 1 ] [ 1 0 0 − 2 1 0 0 0 1 ] = [ 1 0 0 − 2 1 0 10 − 5 0 ] \left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & -5 & 1 \end{matrix} \right] \left[ \begin{matrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right] = \left[\begin{matrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 10 & -5 & 0 \end{matrix} \right] 10001−5001 1−20010001 = 1−21001−5000
左侧矩阵为E(3,2),右侧矩阵为E(2,1),结果矩阵为E。我们可以看到,在消元的过程中出现了非1元素相乘的情况(即10的存在),如果消元矩阵是未知数,那么我们得到的将是二次项。
我们在来看看矩阵L的情况:
[ 1 0 0 2 1 0 0 0 1 ] [ 1 0 0 0 1 0 0 − 5 1 ] = [ 1 0 0 2 1 0 0 5 0 ] \left[ \begin{matrix} 1 & 0 & 0\\ 2 & 1 & 0\\ 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1 \end{matrix} \right] = \left[\begin{matrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 5 & 0 \end{matrix} \right] 120010001 10001−5001 = 120015000
左侧矩阵为E(2,1)的逆,右侧矩阵为E(3,2)的逆,结果矩阵为L。我们可以看到,这次的乘法计算中并未出现非1元素相乘的情况,即结果最多只是一次项。
这就是为什么我们说LU分解计算量小的核心原因。
下面我们就来研究为什么会如此。
首先来看看矩阵乘法中矩阵的顺序,除了矩阵乘法中本身顺序不可改变的原因外,我们等会儿还会看到想明白这一点会给你带来什么样的启发。
因为进行的是3×3的矩阵消元,而一般的方法是从上往下,从左往右进行消元,即矩阵乘法的顺序如下:
A E 2 , 1 A E 3 , 1 E 2 , 1 A − − − 我举例子时为求简洁当做矩阵 ( 3 , 1 ) 处初始就为 0 ,所以没放这一步 E 3 , 2 E 3 , 1 E 2 , 1 A A\\E_{2,1}A\\E_{3,1}E_{2,1}A ---我举例子时为求简洁当做矩阵(3,1)处初始就为0,所以没放这一步\\E_{3,2}E_{3,1}E_{2,1}A AE2,1AE3,1E2,1A−−−我举例子时为求简洁当做矩阵(3,1)处初始就为0,所以没放这一步E3,2E3,1E2,1A
由于这是消元矩阵,除了顺序,仔细想想,我们甚至能确定每个E的大致样貌。
比如E(2,1),由于其任务是要消掉矩阵A(2,1)处的值,所以考虑由矩阵 A 的x倍的第一行与1倍的第二行的线性组合得到U的第二行。
所以,可以确定,E(2,1)应该是如下样子:
提示,如果这里脑筋没转过弯,为什么能确定E(2,1)是这个样子,建议读读我的这篇博文。
[ 1 0 0 x 1 0 0 0 1 ] \left[ \begin{matrix} 1 & 0 & 0\\ x & 1 & 0\\ 0 & 0& 1 \end{matrix} \right] 1x0010001
同理,(忽略E(3,1)的情况下)E(3,2)的样子也应是如下样子:
[ 1 0 0 0 1 0 0 y 1 ] \left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & y& 1 \end{matrix} \right] 10001y001
而此时E(3,2)乘E(2,1)即为:
[ 1 0 0 0 1 0 0 y 1 ] [ 1 0 0 x 1 0 0 0 1 ] = [ 1 0 0 x 1 0 x y y 1 ] \left[ \begin{matrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & y & 1 \end{matrix} \right] \left[ \begin{matrix} 1 & 0 & 0 \\ x & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right] = \left[\begin{matrix} 1 & 0 & 0 \\ x & 1 & 0 \\ xy & y & 1 \end{matrix} \right] 10001y001 1x0010001 = 1xxy01y001
可以看到,这时的(3,1)项一定会产生二次项。
当然,你也可以直接想到有y倍的 [x 1 0] 这一线性组合的部分,所以会产生二次项。
然而,当我们取逆后,计算变为如下:
[ 1 0 0 − x 1 0 0 0 1 ] [ 1 0 0 0 1 0 0 − y 1 ] = [ 1 0 0 − x 1 0 0 − y 1 ] \left[ \begin{matrix} 1 & 0 & 0\\ -x & 1 & 0\\ 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -y & 1 \end{matrix} \right] = \left[\begin{matrix} 1 & 0 & 0 \\ -x & 1 & 0 \\ 0 & -y & 1 \end{matrix} \right] 1−x0010001 10001−y001 = 1−x001−y001
这时,我们看到,计算结果就并不存在二次项了,取逆的重点不在于改变矩阵本身,而是改变了矩阵相乘的顺序,这样就能让本身元素大致分布已知的消元矩阵相乘时错开,不产生二次项,从而达到简化计算的效果。
总结
如上,基本就是LU分解的大致内容。下一篇博文大概会围绕着LU分解的计算复杂度来讲。