<复习向>线性代数之PLU分解

Tags: ,

行列式的求解

从行列式的定义出发去求行列式,是一个简单但低效的方法。而实际上,解决数学问题的途径往往有多种。这里,我将介绍其中一种比较常见的快速解法:PLU分解

PLU的LU

要理解PLU,得先搞懂LU分解。(这里分享一个外教的讲解视频,简单好理解:https://www.youtube.com/watch?v=UlWcofkUDDU 能翻墙的同学就直接看吧。)

LU分别代表:Lower Triangular Matrix 和 Upper Triangular Matrix,即下三角矩阵和上三角矩阵。

下面手动演示下LU分解过程:

设A:

\[ A = \left[ \begin{matrix} 9&6&0\\ 6&5&4\\ 3&4&10\\ \end{matrix} \right] \]

要把A分解成LU,第一步是先用高斯消元法,把A变成阶梯型矩阵U:

  • \( R2 -= R1 * a_{10}/a_{00} = R1 * 6/9 \)

\[ A_{0} = \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 3&4&10\\ \end{matrix} \right] = E_{0}A = \left[ \begin{matrix} 1&0&0\\ -6/9&1&0\\ 0&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 9&6&0\\ 6&5&4\\ 3&4&10\\ \end{matrix} \right] \]

  • \( R3 -= R1 * a_{20}/a_{00} = R1 * 3/9 \)

\[ A_{1} = \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 0&2&10\\ \end{matrix} \right] = E_{1}A_{0} = \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ -3/9&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 3&4&10\\ \end{matrix} \right] \]

  • \( R3 -= R2 * a_{21}/a_{11} = R1 * 2/1 \)

\[ U = A_{2} = \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 0&0&2\\ \end{matrix} \right] = E_{2}A_{1} = \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&-2/1&1\\ \end{matrix} \right] \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 0&2&10\\ \end{matrix} \right] \]

  • 因此得:\( U = E_{2}A_{1} = E_{2}E_{1}A_{0} = E_{2}E_{1}E_{0}A \)

  • 再调换下,得到:\( A = E_{0}^{-1}E_{1}^{-1}E_{2}^{-1}U \)

  • 所以, \( L = E_{0}^{-1}E_{1}^{-1}E_{2}^{-1} \)

要得到最终的L,需要算3个\(E_{x}\)矩阵的逆矩阵,看似麻烦,其实很简单,因为\(E_{x}\)有这样的性质:

\[ \left[ \begin{matrix} 1&0&0\\ a&1&0\\ 0&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ -a&1&0\\ 0&0&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{matrix} \right] \]

\[ \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ a&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ -a&0&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{matrix} \right] \]

\[ \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&a&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&-a&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{matrix} \right] \]

所以:

\[ L = E_{0}^{-1}E_{1}^{-1}E_{2}^{-1} = \left[ \begin{matrix} 1&0&0\\ 6/9&1&0\\ 0&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 3/9&0&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ 0&1&0\\ 0&2/1&1\\ \end{matrix} \right] \]

只要搞定右边的3矩阵乘法运算,就能得到L。而又因为:

\[ \left[ \begin{matrix} 1&0&0\\ a_{1}&1&0\\ b_{1}&c_{1}&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&0&0\\ a_{2}&1&0\\ b_{2}&c_{2}&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0&0\\ a_{1}+a_{2}&1&0\\ b_{1}+a_{2}c_{2}+b_{2}&c_{1}+c_{2}&1\\ \end{matrix} \right] \]

所以,L的结果可以迅速得到:

\[ L = E_{0}^{-1}E_{1}^{-1}E_{2}^{-1} = \left[ \begin{matrix} 1&0&0\\ 6/9&1&0\\ 3/9&2/1&1\\ \end{matrix} \right] \]

于是,A的LU分解完成了:

\[ A = LU = \left[ \begin{matrix} 1&0&0\\ 6/9&1&0\\ 3/9&2/1&1\\ \end{matrix} \right] \left[ \begin{matrix} 9&6&0\\ 0&1&4\\ 0&0&2\\ \end{matrix} \right] \]

PLU的P

这里的P,指的是Permutation Matrix,置换矩阵。

何谓置换矩阵?其实就是经过一系列初等变换的单位矩阵,且元素\( a_{ij} = 0 or 1 \)。

置换矩阵的作用,是用来交换某个矩阵的行(列)顺序。

比如:

\[ P = \left[ \begin{matrix} 0&1&0\\ 0&0&1\\ 1&0&0\\ \end{matrix} \right]\ \ A = \left[ \begin{matrix} 3&4&0\\ 1&2&9\\ 0&5&6\\ \end{matrix} \right] \]

\[ PA = \left[ \begin{matrix} 0&1&0\\ 0&0&1\\ 1&0&0\\ \end{matrix} \right] \left[ \begin{matrix} 3&4&0\\ 1&2&9\\ 0&5&6\\ \end{matrix} \right] = \left[ \begin{matrix} 1&2&9\\ 0&5&6\\ 3&4&0\\ \end{matrix} \right] \]

\[ AP = \left[ \begin{matrix} 3&4&0\\ 1&2&9\\ 0&5&6\\ \end{matrix} \right] \left[ \begin{matrix} 0&1&0\\ 0&0&1\\ 1&0&0\\ \end{matrix} \right] = \left[ \begin{matrix} 0&3&4\\ 9&1&2\\ 6&0&5\\ \end{matrix} \right] \]

从这个例子就可以看出,P左乘A时,改变了A的行的顺序;P右乘A时,改变了A的列的顺序。

PA = LU?

为什么要先对A做P置换后,再做LU分解?这是因为不这样做的话,LU会不稳定(stability)。

举个例子:

\[ A = \left[ \begin{matrix} 10^{-20}&1\\ 1&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ 10^{20}&1\\ \end{matrix} \right] \left[ \begin{matrix} 10^{-20}&1\\ 0&1-10^{20}\\ \end{matrix} \right] = L_{0}U_{0} \]

直接分解后得到的L、U矩阵,出现了大数。所以这是不能接受的。

而神奇的是,对A做一些P置换后,再来LU分解,是可以变稳定的:

\[ PA = \left[ \begin{matrix} 0&1\\ 1&0\\ \end{matrix} \right] \left[ \begin{matrix} 10^{-20}&1\\ 1&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&1\\ 10^{-20}&1\\ \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ 10^{-20}&1\\ \end{matrix} \right] \left[ \begin{matrix} 1&1\\ 0&1-10^{-20}\\ \end{matrix} \right] = L_{1}U_{1} \]

L、U中没有出现大数,于是认为这样的分解是稳定的。

PA的P,需要对A做一些检测后才可以得到,策略就是:沿着对角线从左上角到右下角遍历A,并检测当前列的最大元素在下方的哪一行(当前行上方的行保持不变),找到后就将当前行和目标行交换,并记录下一个\(P_{i}\)。最后按顺序算\(P_{i}\)的乘积就得到了P。

A的行列式

<复习向>线性代数之矩阵与行列式(1)中,已经提到了一条行列式公式:

\[ det(AB) = det(A)det(B) \]

而,\( PA = LU \)又可以变成 \( A = P^{-1}LU \),所以:

\[ det(A) = det(P^{-1}LU) = det(P^{-1})det(L)det(U) \]

可以进一步将这个式子简化:

  • L、U矩阵分别是下三角矩阵和上三角矩阵,它们的行列式等于对角线上元素的乘积
  • L矩阵对角线上的元素都为1

于是有:

\[ det(A) = det(P^{-1})u_{11}u_{22}\cdots u_{nn} \]

因为: \( PP^{-1} = PP^{T} = 1 \),\( det(P^{T}) = det(P) \),所以问题变成了求det(P)。

P怎么求?首先,P相当于多个\(E_i\)矩阵的乘积,而又有\( det(E_i)=-1 \) (行列式的基本性质:交换行列式的两行,行列式变号),所以有:

\[ P = E_t\cdots E_2E_1 \implies \det(P) = \prod^t_{i=1}\det(E_i)=(-1)^t \]

于是:

\[ det(A) = det(P^{-1}LU) = det(P)u_{11}u_{22}\cdots u_{nn} = (-1)^{t}u_{11}u_{22}\cdots u_{nn} \]

总结

矩阵的分解(factorization)有很多种,PA=LU只是其中一种,但此类分解法都离不开高斯消元这把大杀器。理解好高斯消元是关键。

P.S. 已经有人证明了,任何方阵都存在它的PLU分解:http://arxiv.org/pdf/math/0506382v1.pdf

(未经授权禁止转载)
Written on August 31, 2015

博主将十分感谢对本文章的任意金额的打赏^_^