| 矩阵分解 | LU分解 |
|---|---|
| 分解形式 | L（下三角矩阵）、U（上三角矩阵） |
| 目的 | 提高计算效率 |
| 前提 | （1）矩阵A为方阵；（2）矩阵可逆（满秩矩阵）；（3）消元过程中没有0主元出现，也就是消元过程中不能出现行交换的初等变换 |
LU分解其实就是将线性方程组：
A
x
=
b
Ax = b
Ax=b
分解为：
L
U
x
=
b
LUx = b
LUx=b
这样一来就会有：
{
L
y
=
b
U
x
=
y
\begin{cases}Ly = b \\ Ux = y \end{cases}
{Ly=bUx=y
先由下三角矩阵求解出y，再通过上三角矩阵完成对x的求解。
高斯迭代求解其实就是我们在小学时候学习的二元一次方程组的解法。将一个方程中的未知数，用含有另一未知数的代数式来表示。比如：
{
x
1
−
2
x
2
+
x
3
=
0
2
x
2
−
8
x
3
=
8
5
x
1
−
5
x
3
=
10
\left\{\begin{array}{l}x_{1}-2 x_{2}+x_{3}=0 \\ 2 x_{2}-8 x_{3}=8 \\ 5 x_{1}-5 x_{3}=10\end{array}\right.
⎩⎨⎧x1−2x2+x3=02x2−8x3=85x1−5x3=10
写成增广形式则有：
[
1
−
2
1
0
0
2
−
8
8
5
0
−
5
10
]
\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]
⎣⎡105−2201−8−50810⎦⎤
使用高斯迭代求解，先使用式1消去式3中的
x
1
x_1
x1，则有：
[
1
−
2
1
0
0
2
−
8
8
0
10
−
10
10
]
\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10\end{array}\right]
⎣⎡100−22101−8−100810⎦⎤
再使用式2消去式3中的
x
2
x_2
x2，则有：
[
1
−
2
1
0
0
2
−
8
8
0
0
30
−
30
]
\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]
⎣⎡100−2201−83008−30⎦⎤
最后通过这个上三角矩阵，可以很轻松的求解出x：
{
x
1
−
2
x
2
+
x
3
=
0
x
2
−
4
x
3
=
4
x
3
=
−
1
\left\{\begin{array}{l}x_{1}-2 x_{2}+x_{3}=0 \\ x_{2}-4 x_{3}=4 \\ x_{3}=-1\end{array}\right.
⎩⎨⎧x1−2x2+x3=0x2−4x3=4x3=−1
简单进行回代即可。
void Guass(float *A,float *x,float *b,int n)
{
float *At = (float *)malloc(sizeof(float)*n*n);
memcpy(At,A,sizeof(float)*n*n);
float *bt = (float *)malloc(sizeof(float)*n);
memcpy(bt,b,sizeof(float)*n);
for (int i=0;i<n;i++)
{
float pivot = At[i*n+i]; //对角线元素
for (int j=i+1;j<n;j++)
{
float scale=-At[j*n+i]/pivot; //计算消元时所乘的系数
for (int k=0;k<n;k++)
{
At[j*n+k]+=scale*At[i*n+k]; //消去元素，一行每一个元素都要相加
}
bt[j]+=scale*bt[i]; //增广形式，向量也需要操作
}
}
printf("At=\n");
printMat(At,n);
printf("bt=\n");
printVec(bt,n);
//回代解方程
for (int i=n-1;i>=0;i--)
{
float sum=0;
for (int j=i+1;j<n;j++)
{
sum+=At[i*n+j]*x[j]; //此处从最底层开始计算时不会进入此循环
}
x[i]=(bt[i]-sum)/At[i*n+i];
}
printf("x=\n");
printVec(x,n);
free(At);
free(bt);
free(x);
}
执行结果如下图：
可能大家已经看出来了最后消除完得到了一个上三角矩阵，其实高斯迭代最后所得的即为： U x = y Ux=y Ux=y。
但是有个疑问，下三角矩阵L又去哪了？
这个地方就涉及到线代中的矩阵初等变换了，如图所示：
我们在高斯迭代时对于矩阵做了两次初等行变换，也就相当于在矩阵左边乘了两个初等矩阵：
[
1
0
0
0
1
0
0
−
5
1
]
[
1
0
0
0
1
0
−
5
0
1
]
[
1
−
2
1
0
0
2
−
8
8
5
0
−
5
10
]
=
[
1
−
2
1
0
0
2
−
8
8
0
0
30
−
30
]
\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1\end{array}\right]\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]
⎣⎡10001−5001⎦⎤⎣⎡10−5010001⎦⎤⎣⎡105−2201−8−50810⎦⎤=⎣⎡100−2201−83008−30⎦⎤
将消元矩阵合并，用E表示则为：
E
A
=
[
1
0
0
0
1
0
−
5
−
5
1
]
[
1
−
2
1
0
0
2
−
8
8
5
0
−
5
10
]
=
[
1
−
2
1
0
0
2
−
8
8
0
0
30
−
30
]
=
U
E A=\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & -5 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]=U
EA=⎣⎡10−501−5001⎦⎤⎣⎡105−2201−8−50810⎦⎤=⎣⎡100−2201−83008−30⎦⎤=U
两边在左乘
E
−
1
E^{-1}
E−1则有：
A
=
[
1
−
2
1
0
0
2
−
8
8
5
0
−
5
10
]
=
E
−
1
U
=
[
1
0
0
0
1
0
5
5
1
]
[
1
−
2
1
0
0
2
−
8
8
0
0
30
−
30
]
=
L
U
A=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=E^{-1} U=\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]=L U
A=⎣⎡105−2201−8−50810⎦⎤=E−1U=⎣⎡105015001⎦⎤⎣⎡100−2201−83008−30⎦⎤=LU
这样就成功将矩阵A分解为了一个下三角矩阵（L）和一个上三角矩阵（U）的乘积。但是需要注意的是并不是所有矩阵都可以分解，要求参见表1的前提！
void Lu_doolittle(float *A,float *x,float *b,int n)
{
float *At = (float *)malloc(sizeof(float)*n*n);
memcpy(At,A,sizeof(float)*n*n);
float *bt = (float *)malloc(sizeof(float)*n);
memcpy(bt,b,sizeof(float)*n);
//为三角型矩阵L和U分配内存
float *L = (float *)malloc(sizeof(float)*n*n);
float *U = (float *)malloc(sizeof(float)*n*n);
memset(L,0,sizeof(float)*n*n);
memset(U,0,sizeof(float)*n*n);
//一边高斯消元一边记录相应的值，构成L和U
for (int i=0;i<n;i++)
{
float pivot =At[i*n+i];
for (int j=i+1;j<n;j++)
{
float scale=-At[j*n+i]/pivot;
At[j*n+i]=-scale;
for (int k=i+1;k<n;k++)
{
At[j*n+k]+=scale*At[i*n+k];
}
}
}
//给L矩阵赋值
for (int i=0;i<n;i++)
{
L[i*n+i]=1.00; //对角线置为1
for (int j=0;j<i;j++)
L[i*n+j]=At[i*n+j];//给下方小三角赋值即可
}
//给U矩阵赋值
for (int i=0;i<n;i++)
for (int j=i;j<n;j++)
U[i*n+j]=At[i*n+j];
printf("L=\n");
printMat(L,n);
printf("U=\n");
printMat(U,n);
float *y=(float *)malloc(sizeof(float)*n);
//从上向下解出向量y
for (int i=0;i<n;i++)
{
float sum=0;
for (int j=0;j<i;j++)
sum+=L[i*n+j]*y[j];
y[i]=(bt[i]-sum)/L[i*n+i];
}
//从下向上解出x
for (int i=n-1;i>=0;i--)
{
float sum=0;
for (int j=i+1;j<n;j++)
sum+=U[i*n+j]*x[j];
x[i]=(y[i]-sum)/U[i*n+i];
}
printf("x=\n");
printVec(x,n);
}
在进行LU分解时，涉及到除法（计算pilot）， 如果分母太大则会导致不稳定的情况 ，如下所示：
假设三位十进制浮点数下计算：
(
0.001
1.00
1.00
2.00
)
(
x
1
x
2
)
=
(
1.00
3.00
)
\left(\begin{array}{ll}0.001 & 1.00 \\ 1.00 & 2.00\end{array}\right)\left(\begin{array}{l}x_{1} \\ x_{2}\end{array}\right)=\left(\begin{array}{l}1.00 \\ 3.00\end{array}\right)
(0.0011.001.002.00)(x1x2)=(1.003.00)
则有：
l
21
=
a
21
a
11
=
1.00
0.001
=
1000
l_{21}=\frac{a_{21}}{a_{11}}=\frac{1.00}{0.001}=1000
l21=a11a21=0.0011.00=1000
L
A
=
(
1
0
−
1000
1
)
(
0.001
1.00
1.00
2.00
)
=
(
0.001
1.00
0
−
1000
)
=
U
L A=\left(\begin{array}{cc}1 & 0 \\ -1000 & 1\end{array}\right)\left(\begin{array}{cc}0.001 & 1.00 \\ 1.00 & 2.00\end{array}\right)=\left(\begin{array}{cc}0.001 & 1.00 \\ 0 & -1000\end{array}\right)=U
LA=(1−100001)(0.0011.001.002.00)=(0.00101.00−1000)=U
此处的-1000在计算时不够精确，
−
1000
×
1.00
+
1
×
2.00
=
−
0.100
×
1
0
4
+
0.200
×
10
=
−
0.100
×
1
0
4
+
0.0002
×
1
0
4
=
−
0.100
×
1
0
4
=
−
1000
-1000 \times 1.00+1 \times 2.00=-0.100 \times 10^{4}+0.200 \times 10\quad=-0.100 \times 10^{4}+0.0002 \times 10^{4}=-0.100 \times 10^{4}=-1000
−1000×1.00+1×2.00=−0.100×104+0.200×10=−0.100×104+0.0002×104=−0.100×104=−1000
此处的0.0002被舍去了，最后结果为：
x
=
(
0
,
1
)
T
x=(\mathbf{0}, \mathbf{1})^{\mathrm{T}}
x=(0,1)T，与精确值差距很大。
选主元其实就是为了避免分母很小影响精度的情况出现，因此选取消元后的低阶矩阵内绝对值最大的元素作为主元。此处仅接受列主元的LU分解，列主元只做行变换，相比全主元的方式，可以减少选择主元时的计算量，也可以避免记录交换信息。
列选主元的LU分解算法其实就是在进行高斯消元的过程中需要先选主元在消元。在矩阵变换中相当于先左乘了一个初等矩阵进行行交换，然后又乘了一个矩阵进行消元。
如下：
A
=
(
−
3
2
6
10
−
7
0
5
−
1
5
)
A=\left(\begin{array}{ccc}-3 & 2 & 6 \\ 10& -7 & 0 \\ 5 & -1 & 5\end{array}\right)
A=⎝⎛−31052−7−1605⎠⎞
选主元10，交换第一行和第二行（等同于左乘初等矩阵P），然后进行高斯消元再左乘矩阵F，则有：
P
1
=
(
0
1
0
1
0
0
0
0
1
)
,
F
1
=
(
1
3
10
1
−
5
10
1
)
P_{1}=\left(\begin{array}{lll}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1\end{array}\right), \quad F_{1}=\left(\begin{array}{ccc}1 & \\ \frac{3}{10} & 1 & \\ -\frac{5}{10} & & 1\end{array}\right)
P1=⎝⎛010100001⎠⎞,F1=⎝⎛1103−10511⎠⎞
此时：
A
(
1
)
=
F
1
P
1
A
=
(
10
−
7
0
0
−
1
10
6
0
5
2
5
)
A^{(1)}=F_{1} P_{1} A=\left(\begin{array}{ccc}10 & -7 & 0 \\ 0 & -\frac{1}{10} & 6 \\ 0 & \frac{5}{2} & 5\end{array}\right)
A(1)=F1P1A=⎝⎛1000−7−10125065⎠⎞
再次进行选主元，选择
5
2
\frac{5}{2}
25，交换二三行，再进行高斯消元得：
A
(
2
)
=
F
2
P
2
A
(
1
)
=
(
10
−
7
0
5
2
5
31
5
)
\boldsymbol{A}^{(2)}=\boldsymbol{F}_{2} \boldsymbol{P}_{2} \boldsymbol{A}^{(1)}=\left(\begin{array}{rrr}10 & -7 & 0 \\ & \frac{5}{2} & 5 \\ & & \frac{31}{5}\end{array}\right)
A(2)=F2P2A(1)=⎝⎛10−72505531⎠⎞
记
A
(
2
)
A^{(2)}
A(2)为U，则有
F
2
P
2
F
1
P
1
A
=
U
F_{2} P_{2} F_{1} P_{1} A=U
F2P2F1P1A=U，令
L
−
1
L^{-1}
L−1为
F
2
P
2
F
1
P
1
F_{2} P_{2} F_{1} P_{1}
F2P2F1P1可得：
A
=
L
U
=
(
−
3
10
−
1
25
1
1
0
0
1
2
1
0
)
(
10
−
7
0
5
2
5
31
5
)
A=LU=\left(\begin{array}{ccc}-\frac{3}{10} & -\frac{1}{25} & 1 \\ 1 & 0 & 0 \\ \frac{1}{2} & 1 & 0\end{array}\right)\left(\begin{array}{rrr}10 & -7 & 0 \\ & \frac{5}{2} & 5 \\ & \frac{31}{5}\end{array}\right)
A=LU=⎝⎛−103121−25101100⎠⎞⎝⎛10−72553105⎠⎞
伪代码如下所示：
伪代码如下所示：
每个节点处理一列，确定主元后，计算消元所需的值，将其和主元位置进行广播。其他节点接收后先进行主元的交换，然后进行消元。
