线性代数 --- LU分解的编程实现Crout‘s method

松下J27

已于 2023-12-26 15:37:54 修改

阅读量449

点赞数 1

分类专栏： Linear Algebra 文章标签：算法线性代数 lu分解

于 2021-09-27 19:08:42 首次发布

本文链接：https://blog.csdn.net/daduzimama/article/details/120515110

版权

Linear Algebra 专栏收录该内容

68 篇文章 55 订阅

订阅专栏

对于矩阵A而言(A必须是方阵)，假设A可以被分解成两个矩阵L和U的乘积，则有：

$L\cdot U=A$

其中L为下三角矩阵(除了主对角线上的元素和主对角线以下的元素之外，其他元素都是0的矩阵)，U为上三角矩阵(除了主对角线上的元素和主对角线以上的元素之外，其他元素都是0的矩阵)。以4x4矩阵为例，假设我们令下三角矩阵L和上三角矩阵U分别为：

$L=\begin{bmatrix} \alpha _{11}& 0 &0 & 0\\ \alpha _{21} & \alpha _{22} & 0 &0 \\ \alpha _{31}& \alpha _{32} & \alpha _{33} &0 \\ \alpha _{41}& \alpha _{42} & \alpha _{43} & \alpha _{44} \end{bmatrix}$ $U= \begin{bmatrix} \beta _{11} & \beta _{12} &\beta _{13} &\beta _{14} \\ 0 & \beta _{22}& \beta _{23} & \beta _{24}\\ 0 & 0& \beta _{33}&\beta _{34} \\ 0 & 0& 0& \beta _{44} \end{bmatrix}$

则LU分解可被表示为如下形式，其中L和U中的元素均为未知数(除0外)，矩阵A中的元素为已知数：

$\begin{bmatrix} \alpha _{11}& 0 &0 & 0\\ \alpha _{21} & \alpha _{22} & 0 &0 \\ \alpha _{31}& \alpha _{32} & \alpha _{33} &0 \\ \alpha _{41}& \alpha _{42} & \alpha _{43} & \alpha _{44} \end{bmatrix} \begin{bmatrix} \beta _{11} & \beta _{12} &\beta _{13} &\beta _{14} \\ 0 & \beta _{22}& \beta _{23} & \beta _{24}\\ 0 & 0& \beta _{33}&\beta _{34} \\ 0 & 0& 0& \beta _{44} \end{bmatrix}= \begin{bmatrix} a_{11}& a_{12} &a_{13} & a_{14}\\ a_{21}&a_{22} & a_{23}&a_{24} \\ a_{31}& a_{32} &a_{33} &a_{34} \\ a_{41} &a_{42} &a_{43} & a_{44} \end{bmatrix}$

下面按照矩阵乘法的规律逐一写出矩阵A中每一个元素的等式。

对矩阵A中的第一行元素有：

$a_{11}=\alpha _{11}*\beta _{11}$ （1）

$a_{21}=\alpha _{21}*\beta _{11}$ （2）

$a_{31}=\alpha _{31}*\beta _{11}$ （3）

$a_{41}=\alpha _{41}*\beta _{11}$ （4）

对矩阵A中的第二行元素有：

$a_{12}=\alpha _{11}*\beta _{12}$ （5）

$a_{22}=\alpha _{21}*\beta _{12}+\alpha _{22}*\beta _{22}$ （6）

$a_{32}=\alpha _{31}*\beta _{12}+\alpha _{32}*\beta _{22}$ （7）

$a_{42}=\alpha _{41}*\beta _{12}+\alpha _{42}*\beta _{22}$ （8）

对矩阵A中的第三行元素有：

$a_{13}=\alpha _{11}*\beta _{13}$ （9）

$a_{23}=\alpha _{21}*\beta _{13}+\alpha _{22}*\beta _{23}$ （10）

$a_{33}=\alpha _{31}*\beta _{13}+\alpha _{32}*\beta _{23}+\alpha _{33}*\beta _{33}$ （11）

$a_{43}=\alpha _{41}*\beta _{13}+\alpha _{42}*\beta _{23}+\alpha _{43}*\beta _{33}$ （12）

对矩阵A中的第四行元素有：

$a_{14}=\alpha _{11}*\beta _{14}$ （13）

$a_{24}=\alpha _{21}*\beta _{14}+\alpha _{22}*\beta _{24}$ （14）

$a_{34}=\alpha _{31}*\beta _{14}+\alpha _{32}*\beta _{24}+\alpha _{33}*\beta _{34}$ （15）

$a_{44}=\alpha _{41}*\beta _{14}+\alpha _{42}*\beta _{24}+\alpha _{43}*\beta _{34}+\alpha _{44}*\beta _{44}$ （16）

如果矩阵A的尺寸为NxN，则共有NxN个方程和(N+1)xN个未知数。就本例而言，N=4，则共有16个方程和20个未知数。未知数的个数大于方程的个数，则方程组的解不确定。但如果我们暂时只考虑没有行交换的高斯消元过程的话，则最终得到的L矩阵的对角线上一定是1。

这样一来，我们就令 $\alpha _{11},\alpha_{22},\alpha_{33},\alpha_{44}$ 都为1，得到：

$\alpha _{11}=1$ ， $\alpha _{22}=1$ ， $\alpha _{33}=1$ ， $\alpha _{44}=1$

共16个方程，16个未知数。

后续我们还会知道，因为crout method得算法是一种“in place transform”，也就是会破坏原矩阵A，并把最终结果保存到矩阵A中。因此，我们这里可以暂时忽略矩阵L中的0和对角线上的1，同样也忽略矩阵U中的0。得到如下矩阵：

$\begin{bmatrix} \beta _{11}& \beta _{12} & \beta _{13} & \beta _{14}\\ \alpha _{21} & \beta _{22} & \beta _{23} & \beta _{24} \\ \alpha _{31}& \alpha _{32} & \beta _{33} & \beta _{34} \\ \alpha _{41}& \alpha _{42} & \alpha _{43} & \beta _{44}\end{bmatrix}$

求解过程：

$\begin{bmatrix} 1& 0 &0 & 0\\ \alpha _{21} & 1 & 0 &0 \\ \alpha _{31}& \alpha _{32} & 1 &0 \\ \alpha _{41}& \alpha _{42} & \alpha _{43} & 1 \end{bmatrix} \begin{bmatrix} \beta _{11} & \beta _{12} &\beta _{13} &\beta _{14} \\ 0 & \beta _{22}& \beta _{23} & \beta _{24}\\ 0 & 0& \beta _{33}&\beta _{34} \\ 0 & 0& 0& \beta _{44} \end{bmatrix}= \begin{bmatrix} a_{11}& a_{12} &a_{13} & a_{14}\\ a_{21}&a_{22} & a_{23}&a_{24} \\ a_{31}& a_{32} &a_{33} &a_{34} \\ a_{41} &a_{42} &a_{43} & a_{44} \end{bmatrix}$

首先：对于第一列未知数。

$a_{11}=\alpha _{11}*\beta _{11}$ （1）

$a_{21}=\alpha _{21}*\beta _{11}$ （2）

$a_{31}=\alpha _{31}*\beta _{11}$ （3）

$a_{41}=\alpha _{41}*\beta _{11}$ （4）

要想知道 $\alpha _{21},\alpha _{31},\alpha _{41}$ 需要先知道 $\beta _{11}$ ，但唯独式（1）可以求出 $\beta _{11}$ 。把 $\alpha _{11}=1$ 代入（1），可得：

$\beta _{11}=a_{11}$

注意：这一步分别用到了 $a_{11}$ 和 $\alpha _{11}$ (Input用红色方框表示)，求出 $\beta _{11}$ (output用绿色方框表示)，我们在图中标出。

然后再把 $\beta _{11}$ 代入（2）（3）（4），分别求出 $\alpha _{21},\alpha _{31},\alpha _{41}$ ，得到(同样，标出每一步所用到的数据和每一步的输出)：

$\alpha _{21}= a_{21}/\beta _{11}$

$\alpha _{31}=a_{31}/\beta _{11}$

$\alpha _{41}=a_{41}/\beta _{11}$

至此，第一列中的 $\beta$ 和 $\alpha$ 都计算完毕，且是先求出了 $\beta$ ，然后才能求 $\alpha$ 。

对于第二列未知数：

Tips:在求第二列之前，第一列的值都已经知道了。

$a_{12}=\alpha _{11}*\beta _{12}$ （5）

$a_{22}=\alpha _{21}*\beta _{12}+\alpha _{22}*\beta _{22}$ （6）

$a_{32}=\alpha _{31}*\beta _{12}+\alpha _{32}*\beta _{22}$ （7）

$a_{42}=\alpha _{41}*\beta _{12}+\alpha _{42}*\beta _{22}$ （8）

先观察（5）（6）（7）（8）式，除了第一列中已经算出来的未知数，求解（6）（7）（8）式都要用到 $\beta_{12}$ 和 $\beta_{22}$ 。

首先，已知 $\alpha _{11}=1$ ，代入（5）可求出 $\beta_{12}$ ：

$\beta _{12}=a_{12}/\alpha _{11}=a_{12}$

有了 $\beta_{12}$ 再加上 $\alpha _{11}=1$ 和第一列中已经算出的未知数 $\alpha_{21}$ ，代入（6）可求出 $\beta_{22}$ ：

$\beta _{22}=a_{22}-\alpha _{21}*\beta _{12}$

继而，能够分别求出 $\alpha_{32}$ 和 $\alpha_{42}$ ：

$\alpha _{32}=(a_{32}-\alpha _{31}*\beta _{12})/\beta _{22}$

$\alpha _{42}=(a_{42}-\alpha _{41}*\beta _{12})/\beta _{22}$

至此，第二列中的 $\beta$ 和 $\alpha$ 都计算完毕，同样也是先求出了 $\beta$ ，然后才求的 $\alpha$ 。

对于第三列未知数：

Tips:在求第三列之前，第一列和第二列的值都已经知道了。

$a_{13}=\alpha _{11}*\beta _{13}$ （9）

$a_{23}=\alpha _{21}*\beta _{13}+\alpha _{22}*\beta _{23}$ （10）

$a_{33}=\alpha _{31}*\beta _{13}+\alpha _{32}*\beta _{23}+\alpha _{33}*\beta _{33}$ （11）

$a_{43}=\alpha _{41}*\beta _{13}+\alpha _{42}*\beta _{23}+\alpha _{43}*\beta _{33}$ （12）

观察（9）（10）（11）（12）可知，欲求 $\alpha_{43}$ 要先知道 $\beta_{13}$ ， $\beta_{23}$ 和 $\beta_{33}$ 。而欲求 $\beta_{33}$ 又要先知道 $\beta_{13}$ 和 $\beta_{23}$ ，以此类推，可得我们最终的求解顺序应为 $\beta_{13}$ ， $\beta_{23}$ ， $\beta_{33}$ 和 $\alpha_{43}$ 。

已知 $\alpha _{11}=1$ ，代入（9）可求出 $\beta_{13}$ ：

$\beta _{13}=a_{13}/\alpha _{11}=a_{13}$

继而根据（10）（11）算出 $\beta_{23}$ ， $\beta_{33}$ ：

$\beta _{23}=a_{23}-\alpha _{21}*\beta _{13}$

$\beta _{33}=a_{33}-(\alpha _{31}*\beta _{13}+\alpha _{32}*\beta _{23})$

求得当前列所有所有 $\beta$ 后，最终才能求出 $\alpha_{43}$ ：

$\alpha _{43}=(a_{43}-(\alpha _{41}*\beta _{13}+\alpha _{42}*\beta _{23}))/\beta _{33}$

至此，第三列中的 $\beta$ 和 $\alpha$ 都计算完毕。

对于第四列未知数（也就是最后一列未知数）：

Tips:至此，前三列的值都已经求出来了。

$a_{14}=\alpha _{11}*\beta _{14}$ （13）

$a_{24}=\alpha _{21}*\beta _{14}+\alpha _{22}*\beta _{24}$ （14）

$a_{34}=\alpha _{31}*\beta _{14}+\alpha _{32}*\beta _{24}+\alpha _{33}*\beta _{34}$ （15）

$a_{44}=\alpha _{41}*\beta _{14}+\alpha _{42}*\beta _{24}+\alpha _{43}*\beta _{34}+\alpha _{44}*\beta _{44}$ （16）

这一列没有 $\alpha$ ，只需按序求出 $\beta$ 即可。再度把 $\alpha _{11}=1$ ，代入（13）可求出 $\beta_{14}$ ：

$\beta _{14}=a_{14}/\alpha _{11}=a_{14}$

继而才能根据（14）式求出 $\beta_{24}$ ( $\alpha _{22}=1$ )：

$\beta _{24}=a_{24}-\alpha _{21}*\beta _{14}$

依序，根据（15）式求出 $\beta_{34}$ ( $\alpha _{33}=1$ )：

$\beta _{34}=a_{34}-(\alpha _{31}*\beta _{14}+\alpha _{32}*\beta _{24})$

最后，根据（16）式求出 $\beta_{44}$ ( $\alpha _{44}=1$ )：

$\beta _{44}=a_{44}-(\alpha _{41}*\beta _{14}+\alpha _{42}*\beta _{24}+\alpha _{43}*\beta _{34})$

至此，第四列中的 $\beta$ 都计算完毕。

规律：

1，在上述计算过程中，矩阵A中的所有元素 $a_{ij}$ 都只调用了一次。与此同时，被计算出的相应位置的 $\alpha _{ij}$ 或 $\beta _{ij}$ 替换。

2，对每一列而言，在计算 $\alpha$ 之前需要先求出 $\alpha$ 上面的全部 $\beta$ 。

3，在计算第j列时，只会用到j列的左边的值(1~j列)，而不会用到j列右边的值(j+1~N列)。因此，必须要保证在第j列前，1~j列的值都已经计算完毕了。事实上，如果按照我们上面给出的，实际上也是自然而然得到的计算顺序，当你要计算第j列时，j列之前的1~j-1列的值都已经算好了。,

4，对于 $\beta _{ij}$ 的计算而言，可发现如下规律：他需要用到的元素有第j列中 $\beta _{ij}$ 之上的所有 $\beta$ ，和第i行上位于 $\beta _{ij}$ 左边，从i=1开始到次对角线上的全部 $\alpha _{ij}$ 。如下图所示，要想求出"x"处的 $\beta$ ，需要 $\beta$ 上面红色长方形中的 $\beta$ 和 $\beta$ 左边红色长方形中的 $\alpha$ ， $\alpha$ 的选取止于次对角线。(注意： $\beta$ 都在矩阵的上三角中，包括主对角线)

例如：

5，同样，对于 $\alpha _{ij}$ 的计算而言，需要用到的元素有：第i行 $\alpha _{ij}$ 左边的所有 $\alpha$ 元素，和第j列 $\alpha _{ij}$ 上面主对角线之上的所有 $\beta$ 。如下，要想求出"x"处的 $\alpha$ ，需要 $\alpha$ 上面次对角线上方红色长方形中的 $\beta$ (包括对角线上的 $\beta$ )和 $\alpha$ 左边红色长方形中的 $\alpha$ 。(注意： $\alpha$ 都在矩阵的下三角中，不包括主对角线)

例如：

6，对 $\beta$ 而言，所有的下标都满足i<=j。

7，对 $\alpha$ 而言，所有的下标都满足i>j。

结合上面的一些规律，在编程时，主要考虑以下几大点：

1，对每一列而言，按照从上到下的顺序求解每一列，这就保证了对于每一列而言，先求出了 $\beta$ ，然后再求的 $\alpha$ 。

2，对整个矩阵而言，按照从左到右的顺序进行逐列的求解。确保在计算第j列时，1~j-1中的未知数都已经得到了。

3，当i=j时， $\alpha_{ij}$ =1。

4，矩阵A中所有元素 $a_{ij}$ 都只能使用一次。

也就是按照下图中a,b,c,d,e。。。的顺序去设计程序，其中a,c,e,g,i表示 $\beta$ ,b,d,f,h,j表示 $\alpha$ 。

具体实现(A为NXN方阵)：

从第一列(j=1)开始以此类推，对每一列(j=1~n)有：

(1), 按照i=1~j的顺序，依次求出当前列中主对角线上(包括主对角线)的 $\beta _{ij}$ 。

$\beta _{ij}=a _{ij}-\sum_{k=1}^{i-1}\alpha_{ik}\beta _{kj}$

(2), 随着i的增加，对于i=j+1~N的部分，依次求出当前列中对角线下的 $\alpha_{ij}$ 。

$\alpha _{ij}=(a _{ij}-\sum_{k=1}^{j-1}\alpha_{ik}\beta _{kj})/\beta _{jj}$

(3), 对下一列执行上述两步操作。

事实上，在实际实现的时候，有两点和上面不同，也是这个算法的一个亮点，需要单独说明。当i=j时，也就是(1)中的最后一次计算 $\beta _{jj}$ ，实际上是放在(2)中进行的。而(2)又是分两步进行的，先是计算 $(a _{ij}-\sum_{k=1}^{j-1}\alpha_{ik}\beta _{kj})$ ，然后再除 $\beta _{jj}$ 。因为当i=j时，计算 $\alpha_{ij}$ 公式的前半部分和计算 $\beta _{ij}$ 的部分是一样的。把 $\alpha_{ij}$ 公式前半部分的中间计算结果通通放在A中，并替换A。这样一来，就可以顺理成章的把 $\beta _{jj}$ 的计算完成了，同时，也完成了对 $a_{jj}$ 的替换。随后，对i=j+1~n的部分除以 $\beta _{jj}$ ，得到 $\alpha_{ij}$ 。