fatcat 枝头不见绿,蓄势待春风

线性代数:3.LU分解

2018-09-15
fatcat22

什么是LU分解

将一个矩阵 $A$ 表示成 $A = LU$ 的形式,其中 $L$ 是下三角矩阵, $U$ 是上三角矩阵。这就是所谓的LU分解。例如:

上面的例子中 $A$ 是方阵。如果 $A$ 不是方阵呢?这种情况下,我们希望 $L$ 和 $U$ 中有一个是方阵,并且是三角矩阵;而另一个尽量做成三角矩阵。比如这样的形式:

为什么要做LU分解

LU分解其实就是一个运算优化。无论是求解行列式,还是线性方程组等,完全可以直接使用 $A$ 本身来求解。但使用 $LU$ 矩阵的计算量非常小,小到相对来说是可以忽略的。更何况实际编程中,我们还可以将 $LU$ 矩阵保存,即可以用来求行列式,也可以用来求解线性方程组,可谓“一次计算,到处方便”。下面就分别举例看一下。

利用LU分解求行列式

行列式有一个性质,即如果 $A=LU$ ,则 $detA = det(LU) = (detL)(detU)$ 。
行列式还有一个性质,即对于对角矩阵 $B$ ,其行列式的值对角元素的乘积。
因此对于文章开头的例子中的矩阵 $A$ ,其行列式 $detA = (1 * 1 * 1)(1 * (-3) * 4) = -12$

利用LU分解求解线性方程组

如果要求 $Ax = y$ 中的向量 $x$ ,将 $A$ 进行 $LU$ 分解以后,则 $LUx = y$ 。因此求解 $x$ 的问题变成了先求 $Lz = y$ 中的向量 $z$ ,再求 $Ux = z$ 中的向量 $x$ 。乍一看LU分解没有简化问题,反而使问题复杂化了。其实不然,因为 $L$ 和 $U$ 都是三角矩阵,而求解三角矩阵构成的方程组是非常简单的。仍然以文章开头的例子举例,要求

的解 $x$ 。

因为 $A=LU$ ,令 $Ux=z$ ,所以 $Lz = y$ 。因此我们先求向量 $z$ :

由此可以得到一个很简单的线性方程组:

所以 $z = (1, -3, 0)^T$

有了 $z$ 向量以后,接下来我们求解 $Ux = z$ ,即

我们也可以得到一个很简单的线性方程组:

由此可得到线性方程组的解 $x = (-1, 1, 0)^T$。
从上面的过程中可以看出来,进行LU分解以后,整个求解过程都非常简单。

使用python库进行LU分解

虽然这篇文章主要介绍LU分解,但我目前并不打算学习和介绍使用笔算进行LU分解。取而代之的是使用python库求解。

import sympy

A = sympy.Matrix([[2, 4, 2], [1, 5, 2], [4, -1, 9]])
L, U, P = A.LUdecomposition()
#output:
#Matrix([
#[  1,  0, 0],
#[1/2,  1, 0],
#[  2, -3, 1]])
#U
#Matrix([
#[2, 4, 2],
#[0, 3, 1],
#[0, 0, 8]])
#P
#[]
import scipy.linalg

A = scipy.matrix([[2, 4, 2], [1, 5, 2], [4, -1, 9]])
P, L, U = scipy.linalg.lu(A)
#output:
#P:
#array([[ 0.,  0.,  1.],
#       [ 0.,  1.,  0.],
#       [ 1.,  0.,  0.]])
#L:
#array([[ 1.        ,  0.        ,  0.        ],
#       [ 0.25      ,  1.        ,  0.        ],
#       [ 0.5       ,  0.85714286,  1.        ]])
#U:
#array([[ 4.        , -1.        ,  9.        ],
#       [ 0.        ,  5.25      , -0.25      ],
#       [ 0.        ,  0.        , -2.28571429]])

上面代码中,$P$ 称为置换矩阵。

总结

这篇文章介绍了什么是LU分解,更主要的是,为什么要进行LU分解:进行分解后,可以利用分解结果多次求解,从而增加程序效率。


参考资料

  1. 《程序员的数学 - 线性代数》

Similar Posts

Comments