|
文章目录
前言一、矩阵的LU分解(三角分解)定义二、LU分解实现1.matlab代码2.数据验证
总结
前言
矩阵
A
A
A的
L
U
LU
LU分解是将矩阵A表示成一个下三角矩阵
L
L
L和上三角矩阵
U
U
U的乘积，即
A
=
L
U
A=LU
A=LU，这样的结构便于科学计算。
一、矩阵的LU分解(三角分解)定义
设矩阵
A
A
A是
m
×
n
m \times n
m×n阶矩阵；
L
L
L是
m
×
m
m \times m
m×m阶下三角矩阵，主对角元素为1；
U
U
U是
m
×
n
m \times n
m×n阶上三角矩阵。称
L
U
LU
LU为矩阵
A
A
A的一个三角分解，有：
A
m
×
n
=
[
a
11
a
12
a
13
⋯
a
1
n
a
21
a
22
a
23
⋯
a
2
n
a
31
a
32
a
33
⋯
a
3
n
⋮
⋮
⋮
⋱
⋮
a
m
1
a
m
2
a
m
3
⋯
a
m
n
]
(1-1)
\tag{1-1} A_{m\times n}=\begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots&\vdots&\vdots&\ddots&\vdots\\ a_{m1} & a_{m2} & a_{m3} &\cdots& a_{mn}\\ \end{bmatrix}
Am×n=⎣⎢⎢⎢⎢⎢⎡a11a21a31⋮am1a12a22a32⋮am2a13a23a33⋮am3⋯⋯⋯⋱⋯a1na2na3n⋮amn⎦⎥⎥⎥⎥⎥⎤(1-1)
L
m
×
m
=
[
1
0
0
⋯
0
l
21
1
0
⋯
0
l
31
l
32
1
⋯
0
⋮
⋮
⋮
⋱
⋮
l
m
1
l
m
2
l
m
3
⋯
1
]
(1-2)
\tag{1-2} L_{m\times m}= \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 &0& \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ l_{m1} & l_{m2} & l_{m3} &\cdots&1\\ \end{bmatrix}
Lm×m=⎣⎢⎢⎢⎢⎢⎡1l21l31⋮lm101l32⋮lm2001⋮lm3⋯⋯⋯⋱⋯000⋮1⎦⎥⎥⎥⎥⎥⎤(1-2)
U
m
×
n
=
[
u
11
u
12
u
13
⋯
u
1
n
0
u
22
u
23
⋯
u
2
n
0
0
u
33
⋯
u
3
n
⋮
⋮
⋮
⋱
⋮
0
0
0
⋯
u
m
n
]
(1-3)
\tag{1-3} U_{m\times n}=\begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0& u_{33} & \cdots & u_{3n} \\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0 & 0 & 0 &\cdots& u_{mn}\\ \end{bmatrix}
Um×n=⎣⎢⎢⎢⎢⎢⎡u1100⋮0u12u220⋮0u13u23u33⋮0⋯⋯⋯⋱⋯u1nu2nu3n⋮umn⎦⎥⎥⎥⎥⎥⎤(1-3)
A
m
×
n
=
L
m
×
m
U
m
×
n
⇒
[
a
11
a
12
a
13
⋯
a
1
n
a
21
a
22
a
23
⋯
a
2
n
a
31
a
32
a
33
⋯
a
3
n
⋮
⋮
⋮
⋱
⋮
a
m
1
a
m
2
a
m
3
⋯
a
m
n
]
=
[
u
11
u
12
u
13
⋯
u
1
n
l
21
u
11
l
21
u
12
+
u
22
l
21
u
13
+
u
23
⋯
l
21
u
1
n
+
u
2
n
l
31
u
11
l
31
u
12
+
l
32
u
22
l
31
u
13
+
l
32
u
23
+
u
33
⋯
l
31
u
1
n
+
l
32
u
2
n
+
u
3
n
⋮
⋮
⋮
⋱
⋮
l
m
1
u
11
l
m
1
u
12
+
l
m
2
u
22
l
m
1
u
13
+
l
m
2
u
23
+
l
m
3
u
33
⋯
l
m
1
u
1
n
+
l
m
2
u
2
n
+
⋯
+
l
m
,
m
−
1
u
m
−
1
,
n
+
u
m
n
]
A_{m\times n}=L_{m\times m}U_{m\times n} \begin{aligned} \Rightarrow \begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots&\vdots&\vdots&\ddots&\vdots\\ a_{m1} & a_{m2} & a_{m3} &\cdots& a_{mn}\\ \end{bmatrix} =\begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23} & \cdots & l_{21}u_{1n}+u_{2n} \\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33} & \cdots & l_{31}u_{1n}+l_{32}u_{2n}+u_{3n} \\ \vdots&\vdots&\vdots&\ddots&\vdots\\ l_{m1}u_{11} & l_{m1}u_{12}+l_{m2}u_{22} & l_{m1}u_{13}+l_{m2}u_{23}+l_{m3}u_{33} &\cdots&l_{m1}u_{1n}+l_{m2}u_{2n}+\cdots+l_{m,m-1}u_{m-1,n}+u_{mn}\\ \end{bmatrix} \end{aligned}
Am×n=Lm×mUm×n⇒⎣⎢⎢⎢⎢⎢⎡a11a21a31⋮am1a12a22a32⋮am2a13a23a33⋮am3⋯⋯⋯⋱⋯a1na2na3n⋮amn⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡u11l21u11l31u11⋮lm1u11u12l21u12+u22l31u12+l32u22⋮lm1u12+lm2u22u13l21u13+u23l31u13+l32u23+u33⋮lm1u13+lm2u23+lm3u33⋯⋯⋯⋱⋯u1nl21u1n+u2nl31u1n+l32u2n+u3n⋮lm1u1n+lm2u2n+⋯+lm,m−1um−1,n+umn⎦⎥⎥⎥⎥⎥⎤ 由上式，比较矩阵
A
A
A和
L
U
LU
LU中的各项，可以得出通项：
a
i
j
=
{
u
i
j
j
=
1
,
i
=
1
,
2
,
⋯
,
m
∑
k
=
1
i
−
1
l
i
k
u
k
j
+
u
i
j
i
≤
j
,
i
=
2
,
3
,
⋯
m
,
j
=
i
,
i
+
1
,
⋯
,
n
∑
k
=
1
j
−
1
l
i
k
u
k
j
i
>
j
,
j
=
2
,
3
,
⋯
,
n
,
i
=
j
+
1
,
⋯
,
m
(1-5)
a_{ij}= \begin{cases} u_{ij} & j=1,i=1,2,\cdots,m\\ \sum_{k=1}^{i-1} l_{ik}u_{kj}+u_{ij}&i \leq j,i=2,3,\cdots m,j=i,i+1,\cdots,n\\ \sum_{k=1}^{j-1}l_{ik}u_{kj} & i>j,j=2,3,\cdots,n,i=j+1,\cdots,m \end{cases} \tag{1-5}
aij=⎩⎪⎨⎪⎧uij∑k=1i−1likukj+uij∑k=1j−1likukjj=1,i=1,2,⋯,mi≤j,i=2,3,⋯m,j=i,i+1,⋯,ni>j,j=2,3,⋯,n,i=j+1,⋯,m(1-5) 可以求得上三角矩阵
U
m
×
n
U_{m \times n}
Um×n的表达式为：
u
i
j
=
{
a
i
j
−
∑
k
=
1
i
−
1
l
i
k
u
k
j
i
≤
j
0
i
>
j
(1-6)
u_{ij}=\begin{cases} a_{ij}-\sum_{k=1}^{i-1} l_{ik}u_{kj} &i \leq j \\ 0 & i>j \end{cases} \tag{1-6}
uij={aij−∑k=1i−1likukj0i≤ji>j(1-6) 同理，可以求得下三角矩阵
L
m
×
m
L_{m \times m}
Lm×m的表达式为：
l
i
j
=
{
0
i
<
j
1
i
=
j
(
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
u
k
j
)
/
u
j
j
i
>
j
(1-7)
l_{ij}=\begin{cases} 0 &ij \end{cases} \tag{1-7}
lij=⎩⎪⎨⎪⎧01(aij−∑k=1j−1likukj)/ujjij(1-7) 从
L
L
L的表达式(1-7)可以看出，当
u
j
j
=
0
u_{jj}=0
ujj=0时，递推表达式产生了问题。返回表达式(1-4),这时候不难看出，当
u
j
j
=
0
u_{jj}=0
ujj=0时，
l
i
j
l_{ij}
lij可以为任意值,等式恒成立。为了简化计算，当
u
j
j
=
0
u_{jj}=0
ujj=0时，取
l
i
j
=
0
(
i
>
j
）
l_{ij}=0(i>j）
lij=0(i>j）。 当矩阵
A
A
A的各阶主子式
d
e
t
(
A
k
)
=
d
e
t
(
A
(
1
:
k
,
1
:
k
)
)
,
k
=
1
,
2
,
⋯
,
m
det(A_{k})=det(A(1:k,1:k)),k=1,2,\cdots ,m
det(Ak)=det(A(1:k,1:k)),k=1,2,⋯,m不为0时，
L
U
LU
LU分解式唯一；当矩阵
A
A
A的
k
k
k阶主子式为
0
0
0时，
L
U
LU
LU分解式有无数个。通常在实际应用中，我们只需要一组
L
U
LU
LU分解式，为了简化计算，本文中，当
u
j
j
=
0
u_{jj}=0
ujj=0时，取
l
i
j
=
0
(
i
>
j
)
l_{ij}=0(i>j)
lij=0(i>j)。 当然了，当
u
j
j
=
0
u_{jj}=0
ujj=0时，我们也可以通过对矩阵
A
A
A进行初等行交换，使其各阶主子式不为
0
0
0。上述过程相当于对矩阵左乘一个置换矩阵，来找到一组
L
U
LU
LU分解式，即
P
A
=
L
U
PA=LU
PA=LU，
P
P
P为一置换矩阵，此为MATLAB中的做法（函数lu）。
L
U
LU
LU分解的更多信息参考：https://en.wikipedia.org/wiki/LU_decomposition#Closed_formula
二、LU分解实现
1.matlab代码
function [L,U] = lu_decompose(A)
% lu decompose
% L:下三角矩阵
% U:上三角矩阵
% A:输入矩阵
[m,n]=size(A);
L=eye(m);
L(:,1)=A(:,1)/A(1,1);%L第一列赋值
U=zeros(m,n);
U(1,:)=A(1,:);%U第一行赋值
for i=2:m
for j=2:n
if i> A=[1,2,3;4,5,6;7,8,9];
>> [L,U,P]=lu(A)
L =
1 0 0
1/7 1 0
4/7 1/2 1
U =
7 8 9
0 6/7 12/7
0 0 0
P =
0 0 1
1 0 0
0 1 0
>> [L1,U1]=lu_decompose(A)
L1 =
1 0 0
4 1 0
7 2 1
U1 =
1 2 3
0 -3 -6
0 0 0
>> [L1,U1]=lu_decompose(P*A)
L1 =
1 0 0
1/7 1 0
4/7 1/2 1
U1 =
7 8 9
0 6/7 12/7
0 0 0
总结
从矩阵的乘积的性质出发的得到上三角矩阵和下三角矩阵的表达式，并考虑上三角矩阵主对角元素为0的情况下上三角矩阵和下三角矩阵的表达式。
|