一、A = LU
线性代数很多关键的概念实际上就是矩阵的分解(factorization)。原始矩阵
A
A
A 变成两个或三个特殊矩阵的乘积。第一个分解,实际上也是最重要的分解,来自消元法。因子
L
L
L 和
U
U
U 都是三角形矩阵,分解
A
=
L
U
A=LU
A=LU 来自消元法。
矩阵
U
U
U 是上三角矩阵,它的主元都在对角线上,消元步骤将
A
A
A 变为
U
U
U。现在我们要反向这些步骤(将
U
U
U 变为
A
A
A ),通过一个下三角矩阵
L
L
L 就可以。
L
L
L 的元素正好就是乘数
l
i
j
l_{ij}
lij —— 从行
i
i
i 减去乘数乘主元行
j
j
j。
以
2
×
2
2×2
2×2 的矩阵为例,矩阵
A
A
A 有四个元素
2
,
1
,
6
,
8
2,1,6,8
2,1,6,8,要消去的元素是
6
6
6。从行
2
2
2 减去
3
3
3 乘行
1
1
1,这个正向步骤使用消元矩阵
E
21
E_{21}
E21,乘数
l
21
=
3
l_{21}=3
l21=3。从
U
U
U 到
A
A
A 的反向步骤使用
L
=
E
21
−
1
L=E^{-1}_{21}
L=E21−1(行
2
2
2 加上
3
3
3 乘行
1
1
1)。
正向从
A
到
U
:
E
21
A
=
[
1
0
−
3
1
]
[
2
1
6
8
]
=
[
2
1
0
5
]
=
U
反向从
U
到
A
:
E
21
−
1
U
=
[
1
0
3
1
]
[
2
1
0
5
]
=
[
2
1
6
8
]
=
A
正向从\,A\,到\,U:\kern 5ptE_{21}A=\begin{bmatrix}1&0\\-3 &1\end{bmatrix}\begin{bmatrix}2&1\\6&8\end{bmatrix}=\begin{bmatrix}2&1\\0&5\end{bmatrix}=U\\[1ex]反向从\,U\,到\,A:\kern 5ptE_{21}^{-1}U=\begin{bmatrix}1&0\\3&1\end{bmatrix}\begin{bmatrix}2&1\\0&5\end{bmatrix}=\begin{bmatrix}2&1\\6&8\end{bmatrix}=A
正向从A到U:E21A=[1−301][2618]=[2015]=U反向从U到A:E21−1U=[1301][2015]=[2618]=A上面第二行就是分解
L
U
=
A
LU=A
LU=A,将
E
21
−
1
E_{21}^{-1}
E21−1 用
L
L
L 代替。更大的矩阵会有很多
E
′
s
E's
E′s,
L
\pmb L
L 包含它们所有的逆矩阵。
A
A
A 到
U
U
U 的每一个步骤都要左乘一个矩阵
E
i
j
E_{ij}
Eij,将位置
(
i
,
j
)
(i,j)
(i,j) 的元素变为
0
0
0。为了清晰起见,假设没有行交换,如果
A
A
A 是
3
×
3
3×3
3×3 的矩阵,需要在左边依次乘
E
21
E_{21}
E21,
E
31
E_{31}
E31,
E
32
E_{32}
E32,乘数
l
i
j
l_{ij}
lij 将会使得位置
(
2
,
1
)
(2,1)
(2,1),
(
3
,
1
)
(3,1)
(3,1),
(
3
,
2
)
(3,2)
(3,2) 位置处的元素都变为
0
0
0,它们都在对角线下方,消元法在得到一个上三角矩阵后结束。
现在将
E
′
s
E's
E′s 移到另外一边,将它们的逆矩阵乘上
U
U
U:
( E 32 E 31 E 21 ) A = U 变为 A = ( E 21 − 1 E 31 − 1 E 32 − 1 ) U 就是 A = L U ( 2.6.1 ) (E_{32}E_{31}E_{21})A=U\kern 10pt变为\kern 10ptA=(E_{21}^{-1}E_{31}^{-1}E_{32}^{-1})U\kern 10pt就是\kern 10pt\pmb{A=LU}\kern 10pt(2.6.1) (E32E31E21)A=U变为A=(E21−1E31−1E32−1)U就是A=LU(2.6.1)
逆矩阵是反序相乘,三个逆矩阵的乘积就是 L L L。我们得到了 A = L U A=LU A=LU。
二、解释与例子
第一点:每一个逆矩阵
E
−
1
E^{-1}
E−1 都是下三角矩阵。它的非主对角线元素是
l
i
j
l_{ij}
lij,用来恢复
−
l
i
j
-l_{ij}
−lij 产生的减法。
E
E
E 和
E
−
1
E^{-1}
E−1 的主对角线都是
1
1
1。上面的例子中
l
21
=
3
l_{21}=3
l21=3,
E
=
[
1
0
−
3
1
]
E=\begin{bmatrix}1&0\\-3&1\end{bmatrix}
E=[1−301],
L
=
E
−
1
=
[
1
0
3
1
]
L=E^{-1}=\begin{bmatrix}1&0\\3&1\end{bmatrix}
L=E−1=[1301]。
第二点:式(2.6.1)展示了一个下三角矩阵(
E
i
j
E_{ij}
Eij 的乘积)乘
A
A
A,也展示了所有的
E
i
j
−
1
E_{ij}^{-1}
Eij−1 乘
U
U
U 会得到
A
A
A。
E
i
j
\pmb{E_{ij}}
Eij 的逆矩阵的乘积得到的下三角矩阵就是
L
\pmb L
L。
我们处理这些逆矩阵的一个原因是想要分解
A
A
A,而不是
U
U
U。它的 “反向形式” 得到了
A
=
L
U
A=LU
A=LU。另一个原因是我们会得到更多的信息,
L
L
L 是一个很好的矩阵,这也是第三点。
第三点:每个乘数
l
i
j
l_{ij}
lij 可以直接放入
L
L
L 的
(
i
,
j
)
(i,j)
(i,j) 位置,不需要改变。通常矩阵的乘法会将这些位置弄乱,但是在
L
L
L 里不会。因为逆矩阵的正确顺序,使得
l
l
l 没有发生变化。在式(2.6.3)中会给出原因。
第四点:每个
E
−
1
E^{-1}
E−1 的对角线都是
1
1
1,
L
L
L 也是如此。
A = L U 消元过程中没有行交换 . 上三角矩阵 U 的 主元在它的对角线上 . 下三角矩阵 L 的主 元都是 1 ,乘数 l i j 在 L 的对角线下方 . \begin{matrix}A=LU\end{matrix}\kern 15pt\begin{matrix}\pmb{消元过程中没有行交换}.上三角矩阵\,U\,的\\主元在它的对角线上.下三角矩阵\,L\,的主\\元都是\,1\,,乘数\,l_{ij}\,在\,L\,的对角线下方.\end{matrix} A=LU消元过程中没有行交换.上三角矩阵U的主元在它的对角线上.下三角矩阵L的主元都是1,乘数lij在L的对角线下方.
【例1】消元法从行 2 2 2 减去 1 2 \displaystyle\frac{1}{2} 21 乘行 1 1 1,最后一步从行 3 3 3 减去 2 3 \displaystyle\frac{2}{3} 32 乘行 2 2 2。下三角矩阵 L L L 有 l 21 = 1 2 l_{21}=\displaystyle\frac{1}{2} l21=21, l 32 = 2 3 l_{32}=\displaystyle\frac{2}{3} l32=32。 L U LU LU 的乘积得到 A A A: A = [ 2 1 0 1 2 1 0 1 2 ] = [ 1 0 0 1 2 1 0 0 2 3 1 ] [ 2 1 0 0 3 2 1 0 0 4 3 ] = L U A=\begin{bmatrix}2&1&0\\1&2&1\\0&1&2\end{bmatrix}=\begin{bmatrix}1&0&0\\\displaystyle\frac{1}{2}&1&0\\0&\displaystyle\frac{2}{3}&1\end{bmatrix}\begin{bmatrix}2&1&0\\0&\displaystyle\frac{3}{2}&1\\0&0&\displaystyle\frac{4}{3}\end{bmatrix}=LU A= 210121012 = 12100132001 20012300134 =LU因为 A A A 的 ( 3 , 1 ) (3,1) (3,1) 元素是 0 0 0,所以 ( 3 , 1 ) (3,1) (3,1) 位置的乘数是 0 0 0,即不需要进行操作。
【例2】将 A A A 左上角的元素 2 2 2 改为 1 1 1,变成 B B B。则所有的主元都是 1 1 1,所有的乘数也是 1 1 1。保持这种模式,当 B B B 是 4 × 4 4\times4 4×4 的矩阵时: B = [ 1 1 0 0 1 2 1 0 0 1 2 1 0 0 1 2 ] = [ 1 1 1 0 1 1 0 0 1 1 ] [ 1 1 0 0 1 1 0 1 1 1 ] B=\begin{bmatrix}1&1&0&0\\1&2&1&0\\0&1&2&1\\0&0&1&2\end{bmatrix}=\begin{bmatrix}1&&&\\1&1&&\\0&1&1&\\0&0&1&1\end{bmatrix}\begin{bmatrix}1&1&0&0\\&1&1&0\\&&1&1\\&&&1\end{bmatrix} B= 1100121001210012 = 1100110111 1110110011 假设没有行交换,如何可以得知 L L L 和 U U U 中哪些元素为 0 0 0 呢? 当 A 的某一行从 0 开始,则 L 的该行也从 0 开始 当 A 的某一列从 0 开始,则 U 的该列也从 0 开始 \pmb{当\,A\,的某一行从\,0\,开始,则\,L\,的该行也从\,0\,开始}\\\pmb{当\,A\,的某一列从\,0\,开始,则\,U\,的该列也从\,0\,开始} 当A的某一行从0开始,则L的该行也从0开始当A的某一列从0开始,则U的该列也从0开始如果某一行从 0 0 0 开始,那么就不需要消元, L L L 相应的位置就是 0 0 0,这将节省电脑的时间。同样的,如果某一列从 0 0 0 开始, U U U 相应的位置也为 0 0 0。但是,如果 0 0 0 在中间,因为消元法是前向消除,这些位置在 L L L 或 U U U 的对应的位置大概率就不再是 0 0 0。那么为什么 L L L 相应的位置是乘数 l i j l_{ij} lij,而不发生混乱呢?
关键原因是 A \pmb A A 为什么等于 L U \pmb{LU} LU:当主元行下方的行减去时,它们还是 A A A 的原始行吗?不是!因为在消元过程中它们被改变了。那么它们是 U U U 的行吗?是的!因为主元不再改变。当计算 U U U 的第三行时,我们会减去乘数乘 U U U 前面的行(不是 A A A 的行): R o w 3 o f U = ( R o w 3 o f A ) − l 31 ( R o w 1 o f U ) − l 32 ( R o w 2 o f U ) ( 2.6.2 ) Row\,\,3\,\,of\,\,U=(Row\,\,3\,\,of\,\,A)-l_{31}(Row\,\,1\,\,of\,\,U)-l_{32}(Row\,\,2\,\,of\,\,U)\kern 19pt(2.6.2) Row3ofU=(Row3ofA)−l31(Row1ofU)−l32(Row2ofU)(2.6.2)改写这个方程,看看行 [ l 31 l 32 1 ] \begin{bmatrix}l_{31}&l_{32}&1\end{bmatrix} [l31l321] 是如何与 U U U 相乘的: ( R o w 3 o f A ) = l 31 ( R o w 1 o f U ) + l 32 ( R o w 2 o f U ) + 1 ( R o w 3 o f U ) ( 2.6.3 ) (Row\,\,3\,\,of\,\,A)=l_{31}(Row\,\,1\,\,of\,\,U)+l_{32}(Row\,\,2\,\,of\,\,U)+1(Row\,\,3\,\,of\,\,U)\kern 10pt(2.6.3) (Row3ofA)=l31(Row1ofU)+l32(Row2ofU)+1(Row3ofU)(2.6.3)正好就是 ( A = L U ) (A=LU) (A=LU) 的第三行。 L L L 的行 3 3 3 的分量是 l 31 l_{31} l31, l 32 l_{32} l32, 1 1 1。无论 A A A 有多大,所有的行都是这样。如果没有行交换,则有 A = L U A=LU A=LU。
平衡形式 L D U \pmb{LDU} LDU: A = L U A=LU A=LU 是不对称的,因为 U U U 的对角线上是主元,而 L L L 的对角线都是 1 1 1。将 U U U 分成一个举着 D D D 和一个新的矩阵 U U U 相乘,矩阵 D D D 是对角线矩阵,它的对角线就是 U U U 的主元,而新的矩阵 U U U 的对角线就会变成 1 1 1,其它元素除以该行的主元:
新的上三角矩阵
U
U
U 的对角线都是
1
1
1。与正常的
L
U
LU
LU 不同的是,新形式中间有一个
D
D
D:下三角矩阵
L
\pmb L
L 乘对角线矩阵
D
\pmb D
D 乘上三角矩阵
U
\pmb U
U。
矩阵的三角分解可以写 A = L U 或 A = L D U 矩阵的三角分解可以写\kern 10pt\pmb{A=LU}\,或 \,\pmb{A=LDU} 矩阵的三角分解可以写A=LU或A=LDU
当看到 L D U LDU LDU 的形式时,我们就可以知道 U U U 的对角线都是 1 1 1,每一行都除以其第一个非零元素 —— 主元( L U LU LU中 U U U 的主元 ): [ 1 0 3 1 ] [ 2 8 0 5 ] 进一步分解为 [ 1 0 3 1 ] [ 2 5 ] [ 1 4 0 1 ] ( 2.6.4 ) \begin{bmatrix}1&0\\3&1\end{bmatrix}\begin{bmatrix}2&8\\0&5\end{bmatrix}\kern 10pt进一步分解为\kern 10pt\begin{bmatrix}1&0\\3&1\end{bmatrix}\begin{bmatrix}2&\\&5\end{bmatrix}\begin{bmatrix}1&4\\0&1\end{bmatrix}\kern 10pt(2.6.4) [1301][2085]进一步分解为[1301][25][1041](2.6.4)主元 2 2 2 和 5 5 5 进入了 D D D,行分别除以 2 2 2 和 5 5 5,得到新的 U U U,它的对角线都是 1 1 1。乘数 3 3 3 仍然在 L L L 内。
三、一个方形系统 = 两个三角形系统
矩阵
L
L
L 包含了高斯消元法的记忆,它保存了每次进行消元时的乘数。我们可以使用这些求解
A
x
=
b
A\boldsymbol x=\boldsymbol b
Ax=b。
当存在右侧的
b
\boldsymbol b
b 时,则需要
L
L
L。因子
L
L
L 和
U
U
U 完全取决于左侧(矩阵
A
A
A)。在
A
x
=
b
A\boldsymbol x=\boldsymbol b
Ax=b 的右侧,我们先使用
L
−
1
L^{-1}
L−1 再使用
U
−
1
U^{-1}
U−1。该求解步骤会处理两个三角形矩阵。
1、分解(通过对左侧的矩阵 A A A 进行消元,得到 L L L 和 U U U)
2、求解(使用 L L L 对 b \boldsymbol b b 进行前向消元,然后使用 U U U 进行回代求得 x \boldsymbol x x)
以前,我们使用增广矩阵
[
A
b
]
\begin{bmatrix}A&\boldsymbol b\end{bmatrix}
[Ab] 同时处理
A
A
A 和
b
\boldsymbol b
b。但是电脑大多数会将两侧分开,
L
L
L 和
U
U
U 都保存有消元的记忆,无论何时我们都可以处理
b
\boldsymbol b
b。因为这样求解一个单一的系统只需要一个子程序。
那么如何使用
b
\boldsymbol b
b 呢?首先对右侧使用前向消元法(乘数存储在
L
L
L 中,现在可以使用),它会将
b
\boldsymbol b
b 变成一个新的右侧
c
\boldsymbol c
c,我们现在求解的是
L
c
=
b
L\boldsymbol c=\boldsymbol b
Lc=b,然后使用回代求解
U
x
=
c
U\boldsymbol x=\boldsymbol c
Ux=c,原始系统
A
x
=
b
A\boldsymbol x=\boldsymbol b
Ax=b 就被分解成了两个三角系统:
前向与反向 求解 L c = b ,然后求解 U x = c ( 2.6.5 ) \pmb{前向与反向}\kern 10pt求解\kern 5ptL\boldsymbol c=\boldsymbol b,然后求解\kern 5ptU\boldsymbol x=\boldsymbol c\kern 10pt(2.6.5) 前向与反向求解Lc=b,然后求解Ux=c(2.6.5)
要验证
x
\boldsymbol x
x 就是要求的解,
U
x
=
c
U\boldsymbol x=\boldsymbol c
Ux=c 两侧同时左乘
L
L
L,得到
L
U
x
=
L
c
LU\boldsymbol x=L\boldsymbol c
LUx=Lc 就是
A
x
=
b
A\boldsymbol x=\boldsymbol b
Ax=b。
强调: 这些步骤并没有新的知识。我们使用前向消元求解三角系统
L
c
=
b
L\boldsymbol c=\boldsymbol b
Lc=b,然后回代求解
U
x
=
c
U\boldsymbol x=\boldsymbol c
Ux=c。
【例3】以 A x = b A\boldsymbol x=\boldsymbol b Ax=b 前向消元开始,在 U x = c U\boldsymbol x=\boldsymbol c Ux=c 结束: A x = b u + 2 v = 5 4 u + 9 v = 21 变为 u + 2 v = 5 v = 1 U x = c A\boldsymbol x=\boldsymbol b\kern 10pt\begin{matrix}u+2v=5\\4u+9v=21\end{matrix}\kern 10pt变为\kern 10pt\begin{matrix}u+2v=5\\\kern 24ptv=1\end{matrix}\kern 10ptU\boldsymbol x=\boldsymbol c Ax=bu+2v=54u+9v=21变为u+2v=5v=1Ux=c乘数 4 4 4 保存在 L L L 中,右侧使用 4 4 4 将 21 21 21 变成了 1 1 1: L c = b 下三角系统 [ 1 0 4 1 ] [ c ] = [ 5 21 ] 求得 c = [ 5 1 ] L\boldsymbol c=\boldsymbol b\kern 10pt下三角系统\kern 10pt\begin{bmatrix}1&0\\4&1\end{bmatrix}[\boldsymbol c]=\begin{bmatrix}5\\21\end{bmatrix}\kern 5pt求得\kern 5pt\boldsymbol c=\begin{bmatrix}5\\1\end{bmatrix} Lc=b下三角系统[1401][c]=[521]求得c=[51] U x = c 上三角系统 [ 1 2 0 1 ] [ x ] = [ 5 1 ] 求得 x = [ 3 1 ] U\boldsymbol x=\boldsymbol c\kern 10pt上三角系统\kern 10pt\begin{bmatrix}1&2\\0&1\end{bmatrix}[\boldsymbol x]=\begin{bmatrix}5\\1\end{bmatrix}\kern 5pt求得\kern 5pt\boldsymbol x=\begin{bmatrix}3\\1\end{bmatrix} Ux=c上三角系统[1021][x]=[51]求得x=[31] L L L 和 U U U 所使用的也就是以前 A A A 所使用的 n 2 n^2 n2 的存储空间。
四、消元法的成本
这里讨论消元法的成本 —— 即计算时间的问题。我们在计算机上解方程,就需要考虑计算成本,在科学计算时我们可能会遇到大型系统,三维空间的问题就很容易有一百万个未知数,如果计算成本太高的话,我们不可能让计算机计算成百上千年。
消元法的第一阶段是将列
1
1
1 的第一主元以下的元素全部变为
0
0
0,第一行以下的元素全部都需要改变,而改变一个元素需要一次乘法和一次减法,所以第一阶段大约需要
n
2
n^2
n2 次乘法和
n
2
n^2
n2 次减法,实际上会少一些,因为第一行不变,实际上需要
n
(
n
−
1
)
n(n-1)
n(n−1) 次乘法和加法,这里计算的是一个大致成本。
第二阶段我们需要将列
2
2
2 的第二主元下方的元素变为
0
0
0,此时我们要考虑的矩阵会小一些,是一个
(
n
−
1
)
×
(
n
−
1
)
(n-1)\times(n-1)
(n−1)×(n−1) 的矩阵,所以这一阶段大约是
(
n
−
1
)
2
(n-1)^2
(n−1)2 次乘法与减法。越往下进行所要考虑的矩阵越小,最终要得到矩阵
U
U
U 则粗略估计需要的次数为
n
2
+
(
n
−
1
)
2
+
⋯
+
2
2
+
1
2
n^2+(n-1)^2+\cdots+2^2+1^2
n2+(n−1)2+⋯+22+12。
上式平方和的公式为
1
3
n
(
n
+
1
2
)
(
n
+
1
)
\displaystyle\frac{1}{3}n(n+\frac{1}{2})(n+1)
31n(n+21)(n+1),当
n
n
n 很大时,就可以忽略里面的
1
2
\displaystyle\frac{1}{2}
21 和
1
1
1,总和大约就是
1
3
n
3
\displaystyle\frac{1}{3}n^3
31n3。
x
2
x^2
x2 从
0
0
0 到
n
n
n 的积分就是
1
3
n
3
\displaystyle\frac{1}{3}n^3
31n3,需要注意的是积分是连续的,而这里是离散的。
对矩阵 A 使用消元法大概需要 1 3 n 3 次乘法 和 1 3 n 3 次减法 对矩阵 A 使用消元法大概需要 \,\displaystyle\pmb{\frac{1}{3}n^3}\, \pmb{次乘法}和 \displaystyle\frac{1}{3}n^3 次减法 对矩阵A使用消元法大概需要31n3次乘法和31n3次减法
下面考虑右侧的
b
\boldsymbol b
b,我们要计算
L
c
=
b
L\boldsymbol c=\boldsymbol b
Lc=b,得到
c
\boldsymbol c
c。首先,我们从
b
2
,
⋯
,
b
n
b_2,\cdots,b_n
b2,⋯,bn 减去乘数乘
b
1
b_1
b1,这里需要
n
−
1
n-1
n−1 步,第二阶段就不需要考虑
b
1
b_1
b1,共需要
n
−
2
n-2
n−2 步,最后一阶段需要
1
1
1 步。
最后考虑回代,通过
U
x
=
c
U\boldsymbol x=\boldsymbol c
Ux=c 求解
x
\boldsymbol x
x。首先,计算
x
n
x_n
xn 需要
1
1
1 步,仅需要除以最后一个主元;然后计算
x
n
−
1
x_{n-1}
xn−1 需要
2
2
2 步,这里需要代入
x
n
x_n
xn,然后除以第
n
−
1
n-1
n−1 主元;最后计算
x
1
x_1
x1 时需要
n
n
n 步,要代入
(
n
−
1
)
(n-1)
(n−1) 个未知数,然后除以第一主元。所有计算右侧的
b
\boldsymbol b
b 正好需要需要
n
2
n^2
n2 步(从前到后再回代):
[
(
n
−
1
)
+
(
n
−
2
)
+
⋯
+
1
]
+
[
1
+
2
+
⋯
+
(
n
−
1
)
+
n
]
=
n
2
(
2.6.6
)
[(n-1)+(n-2)+\cdots+1]+[1+2+\cdots+(n-1)+n]=n^2\kern 10pt(2.6.6)
[(n−1)+(n−2)+⋯+1]+[1+2+⋯+(n−1)+n]=n2(2.6.6)可以看到,右侧的成本要比左侧小很多。
求解 右侧需要 n 2 次乘法 和 n 2 次减法 \pmb{求解}\kern 15pt右侧需要\,\pmb{n^2\,次乘法}和\,n^2\,次减法 求解右侧需要n2次乘法和n2次减法
一个带状矩阵 B B B 只在主对角线的上方和下方有 w w w 个非零对角线,带状外的其它元素在消元过程中都保持 0 0 0 不变( L L L 和 U U U 中)。第一列需要 w 2 w^2 w2 次乘法和减法(在主元下方产生 w w w 个 0 0 0,每个 0 0 0 需要使用长度为 w w w 的主元行),所以要执行完消元过程共需要 n w 2 nw^2 nw2 次乘法和减法,得到 U U U。这样会节省很多时间。
带状矩阵 A 到 U 1 3 n 3 减少到 n w 2 求解 n 2 减少到 2 n w \pmb{带状矩阵}\kern 14pt\pmb{A\,到\,U}\kern 8pt\frac{1}{3}n^3减少到nw^2\kern 10pt\pmb{求解}\kern 6ptn^2\,减少到\,2nw 带状矩阵A到U31n3减少到nw2求解n2减少到2nw
一个三对角矩阵( w = 1 w=1 w=1)可以计算的很快,不需要存储 0 0 0。
五、主要内容总结
- 高斯消元法(没有行交换)将 A A A 分解成 L L L 乘 U U U。
- 下三角矩阵 L L L 包含用来乘主元行的乘数 l i j l_{ij} lij,它们使得 A A A 变成 U U U。乘积 L U LU LU 将这些行反向加回去可以恢复成 A A A。
- 在右侧我们求解 L c = b L\boldsymbol c=\boldsymbol b Lc=b(前向),然后求解 U x = c U\boldsymbol x=\boldsymbol c Ux=c(反向)。
- 分解: 左侧共有 n 3 − n 3 \displaystyle\frac{n^3-n}{3} 3n3−n 次乘法和减法(这个结果没有取近似)。
- 求解: 右侧共有 n 2 n^2 n2 次乘法和减法。
- 对于带状矩阵,需要的步骤 1 3 n 3 \displaystyle\frac{1}{3}n^3 31n3 变为 n w 2 nw^2 nw2, n 2 n^2 n2 变为 2 n w 2nw 2nw。
六、例题
【例4】下三角帕斯卡矩阵
L
L
L 包含著名的 “帕斯卡三角形”,这里我们分解帕斯卡。
对称帕斯卡矩阵
P
\pmb P
P 是帕斯卡矩阵
L
\pmb L
L 和
U
\pmb U
U 的乘积。对称的
P
P
P 矩阵以帕斯卡三角命名,所以它的每个元素都是其上方和左侧元素之和。MATLAB 中,
n
×
n
n\times n
n×n 的对称
P
P
P 矩阵写成 pascal(n)。
问题:建立一个下 - 上三角分解的
P
=
L
U
P=LU
P=LU。
pascal(4)
=
[
1
1
1
1
1
2
3
4
1
3
6
10
1
4
10
20
]
=
[
1
0
0
0
1
1
0
0
1
2
1
0
1
3
3
1
]
[
1
1
1
1
0
1
2
3
0
0
1
3
0
0
0
1
]
=
L
U
\textrm{pascal(4)}=\begin{bmatrix}\pmb1&\pmb1&\pmb1&\pmb1\\\pmb1&\pmb2&\pmb3&4\\\pmb1&\pmb3&6&10\\\pmb1&4&10&20\\\end{bmatrix}=\begin{bmatrix}\pmb1&0&0&0\\\pmb1&\pmb1&0&0\\\pmb1&\pmb2&\pmb1&0\\\pmb1&\pmb3&\pmb3&\pmb1\end{bmatrix}\begin{bmatrix}\pmb1&\pmb1&\pmb1&\pmb1\\0&\pmb1&\pmb2&\pmb3\\0&0&\pmb1&\pmb3\\0&0&0&\pmb1\end{bmatrix}=LU
pascal(4)=
1111123413610141020
=
1111012300130001
1000110012101331
=LU预测并检验
5
×
5
5\times5
5×5 的帕斯卡矩阵的下一行和下一列。
解: 计算
L
U
LU
LU 可以得到
P
P
P。下面从对称
P
P
P 矩阵出发,利用消元法得到上三角矩阵
U
U
U:
P
=
[
1
1
1
1
1
2
3
4
1
3
6
10
1
4
10
20
]
→
[
1
1
1
1
0
1
2
3
0
2
5
9
0
3
9
19
]
→
[
1
1
1
1
0
1
2
3
0
0
1
3
0
0
3
10
]
→
[
1
1
1
1
0
1
2
3
0
0
1
3
0
0
0
1
]
=
U
P=\begin{bmatrix}1&1&1&1\\1&2&3&4\\1&3&6&10\\1&4&10&20\end{bmatrix}\rightarrow\begin{bmatrix}1&1&1&1\\0&1&2&3\\0&2&5&9\\0&3&9&19\end{bmatrix}\rightarrow\begin{bmatrix}1&1&1&1\\0&1&2&3\\0&0&1&3\\0&0&3&10\end{bmatrix}\rightarrow\begin{bmatrix}1&1&1&1\\0&1&2&3\\0&0&1&3\\0&0&0&1\end{bmatrix}=U
P=
1111123413610141020
→
10001123125913919
→
10001100121313310
→
1000110012101331
=U上面步骤所用到的乘数
l
i
j
l_{ij}
lij 会进入到下三角矩阵
L
L
L,
P
=
L
U
P=LU
P=LU 是一个特别整洁有序的例子。注意到
U
U
U 的在对角线上的主元都是
1
1
1。
若使用 MATLAB 来计算,指令 lu(pascal(4)) 无法生成上述的
U
U
U,这是因为 lu 的子程序会在每一列选取最大的主元,这样第二主元就变成了
3
3
3,而不是
1
1
1,但是使用乔里斯基(Cholesky)分解不会发生行交换,可以产生上述结果:U = chol(pascal(4))
【例5】问题:求解
P
x
=
b
=
(
1
,
0
,
0
,
0
)
P\boldsymbol x=\boldsymbol b=(1,0,0,0)
Px=b=(1,0,0,0)。方程的右侧等于
I
I
I 的第一列,这表明
x
\boldsymbol x
x 会是
P
−
1
P^{-1}
P−1 的第一列。这就是高斯 - 若尔当消元法,会匹配
P
P
−
1
=
I
PP^{-1}=I
PP−1=I 的列。我们已知帕斯卡矩阵
L
L
L 和
U
U
U 是
P
P
P 的两个因子:
两个三角系统
L
c
=
b
(
前向
)
U
x
=
c
(
后向
)
\pmb{两个三角系统}\kern 20ptL\boldsymbol c=\boldsymbol b\,(前向)\kern 10ptU\boldsymbol x=\boldsymbol c\,(后向)
两个三角系统Lc=b(前向)Ux=c(后向)解: 下三角系统
L
c
=
b
L\boldsymbol c=\boldsymbol b
Lc=b 由上到下求解
c
1
=
1
c
1
+
c
2
=
0
c
1
+
2
c
2
+
c
3
=
0
c
1
+
3
c
2
+
3
c
3
+
c
4
=
0
解得
c
1
=
+
1
c
2
=
−
1
c
3
=
+
1
c
4
=
−
1
\begin{matrix}c_1\kern 74pt=1\\c_1+c_2\kern 53pt=0\\c_1+2c_2+c_3\kern 27pt=0\\c_1+3c_2+3c_3+c_4=0\end{matrix}\kern 10pt解得\kern 10pt\begin{matrix}c_1=+1\\c_2=-1\\c_3=+1\\c_4=-1\end{matrix}
c1=1c1+c2=0c1+2c2+c3=0c1+3c2+3c3+c4=0解得c1=+1c2=−1c3=+1c4=−1利用
L
−
1
L^{-1}
L−1 执行前向消元,得到上三角形系统
U
x
=
c
U\boldsymbol x=\boldsymbol c
Ux=c,使用回代求解
x
\boldsymbol x
x,上三角系统由下到上求解:
x
1
+
x
2
+
x
3
+
x
4
=
1
x
2
+
2
x
3
+
3
x
4
=
−
1
x
3
+
3
x
4
=
1
x
4
=
−
1
解得
x
1
=
+
4
x
2
=
−
6
x
3
=
+
4
x
4
=
−
1
\begin{matrix}x_1+x_2+x_3+x_4=1\\\kern 22ptx_2+2x_3+3x_4=-1\\\kern 44ptx_3+3x_4=1\\\kern 81ptx_4=-1\end{matrix}\kern 10pt解得\kern 10pt\begin{matrix}x_1=+4\\x_2=-6\\x_3=+4\\x_4=-1\end{matrix}
x1+x2+x3+x4=1x2+2x3+3x4=−1x3+3x4=1x4=−1解得x1=+4x2=−6x3=+4x4=−1使用 inv(pascal(4)) 指令求得
P
P
P 的逆矩阵,可以看到解就是
P
−
1
P^{-1}
P−1 的第一列。
【例6】对于
6
×
6
6\times6
6×6 的第二差分常对角(constant-diagonal)矩阵
K
K
K,将主元和乘数放入
K
=
L
U
K=LU
K=LU。(
L
L
L 和
U
U
U 有两个非零对角线,因为
K
K
K 有三个)。找到
L
−
1
L^{-1}
L−1 的
i
,
j
i,j
i,j 元素的公式。
**解:
format rat; % 切换成分数显示
K = toeplitz([2 -1 0 0 0 0]) % 生成常对角矩阵,特普利兹 toeplitz 矩阵,这里参数是其第一行
[L U] = lu(K); % 对 K 进行 LU 分解
inv(L) % 求 L 的逆矩阵
从 MATLAB 运行的结果可以得出其
L
−
1
L^{-1}
L−1 的
i
,
j
i,j
i,j 位置元素的公式是:
j
i
,
i
≥
j
\frac{j}{i},\kern 10pti\geq j
ij,i≥j