LU矩阵分解
LU分解可以将一个矩阵
A
\mathbf A
A 分解为一个单位下三角矩阵
L
\mathbf L
L 和一个上三角矩阵
U
\mathbf U
U 的乘积。单位下三角单位矩阵主对角线上的元素为 1 。有如下关系:
A
=
L
U
(1)
\mathbf A = \mathbf L \mathbf U \tag1
A=LU(1)
在程序实现上,(1)式中的矩阵 依次对应到下方三个矩阵 以及 两个三角矩阵相乘后元素之间的关系:
核心代码
算法如下:
for (int i = 0; i < N - 1; i++)
{
double u = U[i * N + i];
for (int j = i + 1; j < N; j++)
{
U[j * N + i] /= u;
for (int k = i + 1; k < N; k++)
{
U[j * N + k] -= U[i * N + k] * U[j * N + i];
}//Of for k
}// Of for j
}
数组 U 初始存储矩阵 A \mathbf A A, 经过计算后存储上三角矩阵 和 除主对角线以外的下三角矩阵 的元素。
在代码的实现上和 图1 的过程有一些区别。代码实现更像是一点一点地把矩阵 A \mathbf A A 给稀释的感觉。
第一行
u
11
:
u
14
u_{11} : u_{14}
u11:u14 可以直接解出。例如在第一轮外循环后,下三角矩阵第一列被解出,上三角第二行被解出。数组 U 中 存储的矩阵元素为:
[
u
11
u
12
u
13
u
14
l
21
u
22
u
23
u
24
l
31
l
32
u
22
l
32
u
23
+
u
33
l
32
u
24
+
u
34
l
41
l
42
u
22
l
42
u
23
+
l
43
u
33
l
42
u
24
+
l
43
u
34
+
u
44
]
\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14}\\ l_{21} & u_{22} & u_{23} & u_{24}\\ l_{31} & l_{32}u_{22} & l_{32}u_{23}+u_{33} &l_{32}u_{24}+u_{34}\\ l_{41} & l_{42}u_{22} &l_{42}u_{23}+l_{43}u_{33} & l_{42}u_{24}+l_{43}u_{34}+u_{44} \end{bmatrix}
⎣⎢⎢⎡u11l21l31l41u12u22l32u22l42u22u13u23l32u23+u33l42u23+l43u33u14u24l32u24+u34l42u24+l43u34+u44⎦⎥⎥⎤
第二轮外循环后,下三角矩阵第二列被解出,上三角矩阵第三行被解出。数组 U 中 存储的矩阵元素为:
[
u
11
u
12
u
13
u
14
l
21
u
22
u
23
u
24
l
31
l
32
u
33
u
34
l
41
l
42
l
43
u
33
l
43
u
34
+
u
44
]
\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14}\\ l_{21} & u_{22} & u_{23} & u_{24}\\ l_{31} & l_{32} & u_{33} & u_{34}\\ l_{41} & l_{42} &l_{43}u_{33} & l_{43}u_{34}+u_{44} \end{bmatrix}
⎣⎢⎢⎡u11l21l31l41u12u22l32l42u13u23u33l43u33u14u24u34l43u34+u44⎦⎥⎥⎤
第三轮外循环…
直到解出 L , U \mathbf L, \mathbf U L,U 矩阵。
LUP矩阵分解
LUP分解是在LU分解的基础上增加主元的选取。
相对于LU分解,LUP分解增加了置换矩阵
P
\mathbf P
P , 有如下关系:
P
A
=
L
U
(2)
\mathbf P \mathbf A = \mathbf L \mathbf U \tag{2}
PA=LU(2)
算法如下:
如 图2 所示,3~6行增加了对矩阵 L , U , P \mathbf L, \mathbf U, \mathbf P L,U,P 的换行,7~12行算法与LU分解相同。
为什么增加了置换矩阵 P \mathbf P P?其实我也是糊涂的(嘿嘿…
大概有以下几个原因:
1.被分解的矩阵
A
\mathbf A
A 其实被要求是非奇异矩阵。再有,如果主元直接为0,LU算法的在计算过程中会出现分母为0的情况。
2.如果主元不为 0 但是接近 0 呢?为了避免除数过小引起精度问题,所以会尽量选择大的主元。
经过行变换后的矩阵再经过LU分解就能尽量避免以上问题。
利用LUP分解对矩阵求逆
由(2)式,有:
P − 1 A − 1 = L − 1 U − 1 \mathbf P^{-1} \mathbf A^{-1} = \mathbf L^{-1} \mathbf U^{-1} P−1A−1=L−1U−1
L , U \mathbf L, \mathbf U L,U 为三角矩阵, P \mathbf P P 为置换矩阵,到这里其实就已经容易求出逆矩阵了。