上一篇我们对下三角矩阵的求解给出了一个方便的求解,利用消元代入可以在
Θ
(
N
2
)
\Theta(N^2)
Θ(N2) 的时间内完成,对于上三角矩阵,我们仍然可以利用类似的方法在相同的时间内求解。
对于一个非三角矩阵系统何对其进行求解,我们将在接下来几篇博文中进行讨论,而在这篇博里我们会对求解做一个浅显的分析和简易尝试。
【1】下三角矩阵的线性方程求解
【2】矩阵的LU分解初步:一个对角线上元素非零的方阵◀
对于如下的线性方程(其中矩阵对角线上的元素不为零,矩阵非奇异)
[
a
11
a
12
.
.
.
a
1
n
a
21
a
22
.
.
.
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
a
m
1
a
m
2
a
m
3
.
.
a
m
n
]
[
x
1
x
2
.
.
.
x
n
]
=
[
b
1
b
2
.
.
.
b
n
]
\begin{bmatrix} a_{11} & a_{12}&. &. &. &a_{1n}\\ a_{21} &a_{22} &. &. &. & a_{2n}\\ .& .& .& & & .\\ .& .& & .& & .\\ .& .& & & .&. \\ a_{m1}&a_{m2} &a_{m3} &. &. &a_{mn} \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2}\\.\\.\\.\\ x_{n}\end{bmatrix}= \begin{bmatrix}b_{1}\\ b_{2}\\.\\ .\\ .\\ b_{n}\end{bmatrix}
⎣⎢⎢⎢⎢⎢⎢⎡a11a21...am1a12a22...am2...am3........a1na2n...amn⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎡x1x2...xn⎦⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡b1b2...bn⎦⎥⎥⎥⎥⎥⎥⎤
我们记
A
=
(
a
i
j
)
A=(a_{ij})
A=(aij),
x
=
(
x
i
)
x=(x_{i})
x=(xi),
b
=
(
b
i
)
b=(b_{i})
b=(bi),于是我们可以简单地将方程记作
A
x
=
b
Ax=b
Ax=b
为了求解这个方程,我们定义两个三角矩阵
L
L
L、
U
U
U其中,
L
L
L为一个单位下三角矩阵,
U
U
U为一个上三角矩阵。如果有
L
U
=
A
LU=A
LU=A成立,那么我们将这一对
L
L
L、
U
U
U称作
A
A
A的一个
L
U
LU
LU分解,这样线性方程又可以记作
L
U
x
=
b
LUx=b
LUx=b在这里,我们令
U
x
=
y
Ux=y
Ux=y,这样上式变成
L
y
=
b
Ly=b
Ly=b如此一来,好处显而易见:我们将一个非三角线性系统转化为两个三角线性系统,这样我们就可以采用上一篇博文中介绍的方法求解。
以上我们简要介绍了
L
U
LU
LU分解的思路。为了对
A
A
A完成
L
U
LU
LU分解,我们可以尝试高斯消元(Gaussian elimination)的方式,实际上
L
U
LU
LU分解确实是一种特殊的高斯消元。
在《算法导论》(第三版,480页)我们可以看到其利用“主元”的方式(俺们小学时候解二元一次方程组的那种方法,同乘上一个系数,一减就完事了,这里写的好高大上啊…)参考其推导如下:
KaTeX parse error: Undefined control sequence: \ at position 281: …1 \end{bmatrix}\̲ ̲
而我们逆用矩阵乘法公式可以继续分解为两个矩阵相乘(可以试着验真)
[
a
11
ω
T
v
A
1
]
=
[
1
0
v
a
11
1
]
[
a
11
ω
T
0
A
1
−
v
ω
T
a
11
]
\begin{bmatrix} a_{11} & \omega^T\\ v& A^ 1 \end{bmatrix}=\begin{bmatrix} 1 & 0\\ \frac{v}{a_{11}}&1 \end{bmatrix}\begin{bmatrix} a_{11} & \omega^T\\ 0& A^ 1- \frac{v \omega^T}{a_{11}} \end{bmatrix}
[a11vωTA1]=[1a11v01][a110ωTA1−a11vωT]
这里我们采用
a
11
a_{11}
a11为主元,而我们又令
A
1
−
v
ω
T
a
11
=
L
1
U
1
A^ 1- \frac{v \omega^T}{a_{11}}=L^1U^1
A1−a11vωT=L1U1,从而又可以得到
A
=
[
1
0
v
a
11
L
1
]
[
a
11
ω
T
0
U
1
]
A=\begin{bmatrix} 1 & 0\\ \frac{v}{a_{11}}&L^1 \end{bmatrix}\begin{bmatrix} a_{11} & \omega^T\\ 0& U^1 \end{bmatrix}
A=[1a11v0L1][a110ωTU1]
这个过程我们可以使用递归完成看出其递归性,我们能够使用递归算法来完成分解。递归是个老生常谈的问题,形式简洁但是使用起来可能会存在一些效率问题,我们可以将其重写为循环的形式,以下给出核心代码
for (i = 0; i < Row - 1; i++)
{
for (j = i; j < Row - 1; j++)
MatA[j + 1][i] /= MatA[i][i];
for (j = i; j < Row - 1; j++)
for (k = i; k < Row - 1; k++)
MatA[j + 1][k + 1] -= MatA[i][k+1] * MatA[j+1][i];
}
这段代码中我们得到了一个包含 L U L U LU的新矩阵,实际上在这个矩阵中 L L L和 U U U以对角线分隔
for (i = 0; i<Row; i++)
{
L[i][i] = 1;//L是一个单位下三角矩阵,因此对角线时1
for (j = 0; j <=i; j++)
L[i][j] = MatA[i][j];
}
for (i = 0; i<Row; i++)
{
for (j = 0; j <= i; j++)
U[j][i] = MatA[j][i];
}