背景:
求解一些列具有相同系数矩阵的线性方程,如: A x = b 1 Ax=b_1 Ax=b1, A x = b 2 Ax=b_2 Ax=b2,… A x = b p Ax=b_p Ax=bp等。当矩阵 A A A可逆时,可以先求出矩阵 A A A的逆 A − 1 A^{-1} A−1,再计算 A − 1 b 1 A^{-1}b_1 A−1b1, A − 1 b 2 A^{-1}b_2 A−1b2等。
但是,也可以用LU分解法来解这一系列方程:先使用初等行变换化简解出 A x = b 1 Ax=b_1 Ax=b1,并同时得到矩阵 A A A的LU分解,剩下的方程使用LU分解法求解即可。
方法:
设矩阵 A A A是 m × n m \times n m×n的矩阵,它可以使用行变换(不进行行的交换)化简为阶梯阵,则矩阵 A A A可以写成 A = L U A=LU A=LU的形式。 L L L是 m × m m\times m m×m的单位下三角矩阵(即主对角线元素全为1的下三角矩阵); U U U是 m × n m\times n m×n的上三角阶梯形矩阵。
用途
当
A
=
L
U
A=LU
A=LU时,方程
A
x
=
b
Ax=b
Ax=b可以写成
L
(
U
x
)
=
b
L(Ux)=b
L(Ux)=b,令
U
x
=
y
Ux=y
Ux=y,则有:
L
y
=
b
U
x
=
y
Ly=b \qquad Ux=y
Ly=bUx=y
即先求解 L y = b Ly=b Ly=b得到 y y y,再求解 U x = y Ux=y Ux=y得到 x x x,因为矩阵 L L L和 U U U都是三角矩阵,所以求解上述两个方程比直接求解 A x = b Ax=b Ax=b要简单。
例题
(1)对矩阵
A
=
[
−
5
3
4
10
−
8
−
9
15
1
2
]
A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}
A=⎣⎡−510153−814−92⎦⎤进行LU分解。
解:
A
=
[
−
5
3
4
10
−
8
−
9
15
1
2
]
A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}
A=⎣⎡−510153−814−92⎦⎤有3行,所以,
L
L
L矩阵应该为
3
×
3
3\times 3
3×3的矩阵。
先求 U U U矩阵:
先把矩阵 A A A一步步消成上三角矩阵就能得到 U U U矩阵:
A = [ − 5 3 4 10 − 8 − 9 15 1 2 ] → [ − 5 3 4 0 − 2 − 1 0 10 14 ] → [ − 5 3 4 0 − 2 − 1 0 0 9 ] = U A=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&10&14\end{bmatrix}\rightarrow \begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=U A=⎣⎡−510153−814−92⎦⎤→⎣⎡−5003−2104−114⎦⎤→⎣⎡−5003−204−19⎦⎤=U
再求 L L L矩阵:
把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵,然后每列除以此列的主元即得到 L L L矩阵:
[ − 5 0 0 10 − 2 0 15 10 9 ] → [ 1 0 0 − 2 1 0 − 3 − 5 1 ] = L \begin{bmatrix}-5&0&0\\10&-2&0\\15&10&9\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}=L ⎣⎡−510150−210009⎦⎤→⎣⎡1−2−301−5001⎦⎤=L
(-5, 10, 15这一列每个元素除以主元-5;-2,10这一列每个元素除以主元-2;9这一列每个元素除以主元9,即得到 L L L矩阵)
验算:
L
U
=
[
1
0
0
−
2
1
0
−
3
−
5
1
]
[
−
5
3
4
0
−
2
−
1
0
0
9
]
=
[
−
5
3
4
10
−
8
−
9
15
1
2
]
=
A
LU=\begin{bmatrix}1&0&0\\-2&1&0\\-3&-5&1\end{bmatrix}\begin{bmatrix}-5&3&4\\0&-2&-1\\0&0&9\end{bmatrix}=\begin{bmatrix}-5&3&4\\10&-8&-9\\15&1&2\end{bmatrix}=A
LU=⎣⎡1−2−301−5001⎦⎤⎣⎡−5003−204−19⎦⎤=⎣⎡−510153−814−92⎦⎤=A
所以上述分解是正确的。
(2)对矩阵 A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix} A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤进行LU分解。
解:
A
=
[
2
4
−
1
5
−
2
−
4
−
5
3
−
8
1
2
−
5
−
4
1
8
−
6
0
7
−
3
1
]
A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}
A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤有4行,所以,
L
L
L矩阵应该为
4
×
4
4\times 4
4×4的矩阵,
L
L
L矩阵的第一列应该是矩阵
A
A
A的第一列除以其主元元素(2)的结果:
L = [ 1 0 0 0 − 2 1 0 0 1 1 0 − 3 1 ] L=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&\space&1&0\\-3&\space&\space&1\end{bmatrix} L=⎣⎢⎢⎡1−21−301 001 0001⎦⎥⎥⎤
先求 U U U矩阵:
先把矩阵 A A A一步步消成上三角矩阵就能得到 U U U矩阵:
A = [ 2 4 − 1 5 − 2 − 4 − 5 3 − 8 1 2 − 5 − 4 1 8 − 6 0 7 − 3 1 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 − 9 − 3 − 4 10 0 12 4 12 − 5 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 4 7 ] → [ 2 4 − 1 5 − 2 0 3 1 2 − 3 0 0 0 2 1 0 0 0 0 5 ] = U A=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&-9&-3&-4&10\\0&12&4&12&-5\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&4&7\end{bmatrix}\rightarrow \begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=U A=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤→⎣⎢⎢⎡200043−912−11−3452−412−2−310−5⎦⎥⎥⎤→⎣⎢⎢⎡20004300−11005224−2−317⎦⎥⎥⎤→⎣⎢⎢⎡20004300−11005220−2−315⎦⎥⎥⎤=U
再求 L L L矩阵:
把上面消元过程中的主元列的主元及其下面的元素都放入一个矩阵,然后每列除以此列的主元就得到 L L L矩阵:
[ 2 0 0 0 − 4 3 0 0 2 − 9 2 0 − 6 12 4 5 ] → [ 1 0 0 0 − 2 1 0 0 1 − 3 1 0 − 3 4 2 1 ] = L \begin{bmatrix}2&0&0&0\\-4&3&0&0\\2&-9&2&0\\-6&12&4&5\end{bmatrix}\rightarrow \begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}=L ⎣⎢⎢⎡2−42−603−91200240005⎦⎥⎥⎤→⎣⎢⎢⎡1−21−301−3400120001⎦⎥⎥⎤=L
(2,-,4, 2, -6这一列每个元素除以主元2;3,-9,12这一列每个元素除以主元3;2,4这一列每个元素除以主元2,5这一列除以主元5,即得到 L L L矩阵)
验算:
L
U
=
[
1
0
0
0
−
2
1
0
0
1
−
3
1
0
−
3
4
2
1
]
[
2
4
−
1
5
−
2
0
3
1
2
−
3
0
0
0
2
1
0
0
0
0
5
]
=
[
2
4
−
1
5
−
2
−
4
−
5
3
−
8
1
2
−
5
−
4
1
8
−
6
0
7
−
3
1
]
=
A
LU=\begin{bmatrix}1&0&0&0\\-2&1&0&0\\1&-3&1&0\\-3&4&2&1\end{bmatrix}\begin{bmatrix}2&4&-1&5&-2\\0&3&1&2&-3\\0&0&0&2&1\\0&0&0&0&5\end{bmatrix}=\begin{bmatrix}2&4&-1&5&-2\\-4&-5&3&-8&1\\2&-5&-4&1&8\\-6&0&7&-3&1\end{bmatrix}=A
LU=⎣⎢⎢⎡1−21−301−3400120001⎦⎥⎥⎤⎣⎢⎢⎡20004300−11005220−2−315⎦⎥⎥⎤=⎣⎢⎢⎡2−42−64−5−50−13−475−81−3−2181⎦⎥⎥⎤=A
所以上述分解是正确的。
关于MATLAB中的LU分解函数
MATLAB中提供了LU分解的函数lu( ),有两种用法,一种是返回三个矩阵(L,U和置换矩阵P),一种是只返回L和U矩阵。例如:
F =
1 3 0
4 4 8
1 2 3
>> [L1 U1 P] = lu(F) //使用MATLAB内置的lu( )函数求矩阵F的LU分解,返回三个矩阵L,U和P的用法
L1 =
1.0000 0 0
0.2500 1.0000 0
0.2500 0.5000 1.0000
U1 =
4 4 8
0 2 -2
0 0 2
P =
0 1 0
1 0 0
0 0 1
>> [L2 U2] = lu(F) //使用MATLAB内置的lu( )函数求矩阵F的LU分解,返回2个矩阵L和U的用法
L2 =
0.2500 1.0000 0
1.0000 0 0
0.2500 0.5000 1.0000
U2 =
4 4 8
0 2 -2
0 0 2
对于上述矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡141342083⎦⎤,手动LU分解为:
L U = [ 1 0 0 4 1 0 1 1 8 1 ] [ 1 3 0 0 − 8 8 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F LU=\begin{bmatrix}1&0&0\\4&1&0\\1&\frac{1}{8}&1\end{bmatrix}\begin{bmatrix}1&3&0\\0&-8&8\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F LU=⎣⎡1410181001⎦⎤⎣⎡1003−80082⎦⎤=⎣⎡141342083⎦⎤=F
上面MATLAB第一种方法(L1,U1和P三个矩阵)算出来的结果:
L 1 U 1 = [ 1 0 0 1 4 1 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 4 4 8 1 3 0 1 2 3 ] ≠ F L_1U_1=\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix}\neq F L1U1=⎣⎡141410121001⎦⎤⎣⎡4004208−22⎦⎤=⎣⎡411432803⎦⎤=F
MATLAB第二种方法(L2和U2 两个矩阵)算出来的结果:
L 2 U 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] [ 4 4 8 0 2 − 2 0 0 2 ] = [ 1 3 0 4 4 8 1 2 3 ] = F L_2U_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}\begin{bmatrix}4&4&8\\0&2&-2\\0&0&2\end{bmatrix}=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=F L2U2=⎣⎡411411021001⎦⎤⎣⎡4004208−22⎦⎤=⎣⎡141342083⎦⎤=F
可见输出两个矩阵的方法得到LU矩阵是正确的,第一种方法得到的 L 1 U 1 ≠ F L_1U_1\neq F L1U1=F。两种方法得到的LU矩阵都和手动算的结果不一样。
原因:
第一种方法:输入[L1 U1 P] = lu(F)
命令的时候,MATLAB首先生成
P
F
=
[
4
4
8
1
3
0
1
2
3
]
PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix}
PF=⎣⎡411432803⎦⎤,其中
P
P
P为置换矩阵(Permutation Matrix),
P
=
[
0
1
0
1
0
0
0
0
1
]
P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}
P=⎣⎡010100001⎦⎤,所以,这种方法求得的分解是
L
1
U
1
=
P
F
L_1U_1=PF
L1U1=PF而不是
L
1
U
1
=
F
L_1U_1=F
L1U1=F。
第二种方法:输入[L2 U2] = lu(F)
命令的时候,算出来的
U
2
U_2
U2矩阵和第一种方法一样(但是和手算的不一样);算出来的
L
2
L_2
L2矩阵其实是
L
2
=
P
T
L
1
L_2=P^TL_1
L2=PTL1,即:
L 2 = [ 1 4 1 0 1 0 0 1 4 1 2 1 ] = P T L 1 = [ 0 1 0 1 0 0 0 0 1 ] [ 1 0 0 1 4 1 0 1 4 1 2 1 ] L_2=\begin{bmatrix}\frac{1}{4}&1&0\\1&0&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix}=P^TL_1=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&0&0\\\frac{1}{4}&1&0\\\frac{1}{4}&\frac{1}{2}&1\end{bmatrix} L2=⎣⎡411411021001⎦⎤=PTL1=⎣⎡010100001⎦⎤⎣⎡141410121001⎦⎤
第二种方法算出来的 L 2 U 2 = F L_2U_2=F L2U2=F是成立的。
关于置换矩阵(Permutation Matrix)
上MATLAB计算LU分解时,两种方法都用到了置换矩阵
P
P
P。
置换矩阵是(0,1)矩阵,即只有0和1两种元素的矩阵,其作用是让矩阵的行进交换的矩阵。置换矩阵
P
P
P是改变相应阶数的单位矩阵得到的,例如对单位矩阵
I
I
I交换第一第二行得到
P
21
=
[
0
1
0
1
0
0
0
0
1
]
P_{21}=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}
P21=⎣⎡010100001⎦⎤,那么得到的置换矩阵
P
21
P_{21}
P21的功能就是,使得和它相乘的矩阵也交换第一第二行。
例如:
P
21
=
[
0
1
0
1
0
0
0
0
1
]
P_{21}=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}
P21=⎣⎡010100001⎦⎤对于上述矩阵
F
=
[
1
3
0
4
4
8
1
2
3
]
F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}
F=⎣⎡141342083⎦⎤,让
P
21
P_{21}
P21乘以
F
F
F就是把
F
F
F的第二行和第一行交换:
P
21
F
=
[
0
1
0
1
0
0
0
0
1
]
[
1
3
0
4
4
8
1
2
3
]
=
[
4
4
8
1
3
0
1
2
3
]
P_{21}F=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix}=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix}
P21F=⎣⎡010100001⎦⎤⎣⎡141342083⎦⎤=⎣⎡411432803⎦⎤
类似的,还有:
P 32 = [ 1 0 0 0 0 1 0 1 0 ] P_{32}=\begin{bmatrix}1&0&0\\0&0&1\\0&1&0\end{bmatrix} P32=⎣⎡100001010⎦⎤
P 31 = [ 0 0 1 0 1 0 1 0 0 ] P_{31}=\begin{bmatrix}0&0&1\\0&1&0\\1&0&0\end{bmatrix} P31=⎣⎡001010100⎦⎤
关于部分主元消去法(Partial Pivoting)
部分主元消去法在求解线性系统时常被用到(如LU分解和使用A\b
求
A
x
=
b
Ax=b
Ax=b等),在计算机的浮点数计算中能提高精确度。
例如,上述第一种方法求LU分解时,对于矩阵 F = [ 1 3 0 4 4 8 1 2 3 ] F=\begin{bmatrix}1&3&0\\4&4&8\\1&2&3\end{bmatrix} F=⎣⎡141342083⎦⎤,MATLAB首先生成 P F = [ 4 4 8 1 3 0 1 2 3 ] PF=\begin{bmatrix}4&4&8\\1&3&0\\1&2&3\end{bmatrix} PF=⎣⎡411432803⎦⎤,把 F F F的第一第二行进行交换,让第一列最大的元素4放在主元位置。
类似的,使用部分主元消去法时,每一行主元位置的元素与其下面的元素中最大的元素放在主元位置(使用行交换将最大的元素换到主元位置),后面使用高斯消元法,这就是部分主元消去法。所以,上面求解中MATLAB使用部分主元消去法时,先使用置换矩阵 P P P来调整主元使之最大,然后再求解。
部分主元消去法就是多了一个把最大的值调到主元位置的操作,却能提高计算精度。