LDLT分解法

LDL T 分解法求解对称线性方程组
朱松盛 041002045 南京师范大学
一、 LDL T 分解法原理
利用矩阵 A 的T LDL 分解来解Ax = b ，的方法称为T LDL 分解法。

当求解方程组的系数矩阵是对称矩阵时，则用T LDL 分解法可以简化程序设计并减少计算量。

当矩阵 A 的各阶顺序主子式不为零时，A 有唯一的 Doolittle 分解 A = LU ，其中
⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=-1111121323121nn n n l l l l l l L ，⎪⎪⎪
⎪⎪
⎪⎭⎫
⎝⎛=nn n n n u u u u u u u u u u U 333223*********。

此时，当然有01
≠=
∏=n
i ii
u
A det ，所以矩阵U 的对角线元素0≠ii u ，
（ , , 2 , 1n i =），将矩阵U 每行依次提出ii u ，则有U D U ~
=，其中
⎪⎪⎪⎪⎪⎭⎫
⎝
⎛=nn u u u D
2211
，⎪⎪⎪⎪
⎪⎪⎪
⎪⎪⎪
⎪
⎪⎪⎭⎫
⎝
⎛=---11
11
11122
2222311111131112n n n n n n u u u u u u u u u u u u U
~
因而U LD A ~
=，显然这种A 的分解也是唯一的，当A 为对称矩阵时，由
A A T =，得()
U LD DL U L D U U LD T T T T T T ~
~~~===，由分解的唯一性，有
L U T =~，即T L U =~。

由此可得，若对称矩阵A 的各阶顺序主子式不为零时，A 可唯一分解为
T LDL A =，其中
⎪⎪⎪⎪⎪
⎪⎭⎫ ⎝⎛=-1111121323121nn n n l l l l l l L ，⎪⎪⎪⎪⎪⎭
⎫
⎝⎛=n d d d D 21，T L 为L 的转置矩阵。

当A 有T LDL 分解时，利用矩阵运算法则及相等原理易得计算ik l 及k d 的公式为
⎪⎪⎩
⎪⎪⎨⎧-=-=∑∑-=-=k
k m m km im ik ik k m m km kk k d
d l l a l d l a d /)(1
11
1
2 ， , , 2 , 1 , , 2 , 1n k k i n k ++== 为减少乘法次数，引入辅助量k ik ik d l u =，则上面公式可写成
⎪⎪⎪
⎩
⎪
⎪
⎪⎨⎧
=-=-=∑∑-=-=k ik ik k m km im ik ik k m km
km kk k d u l l u a u l u a d / 1111
，
n k k i n k , , 2 , 1 , , 2 , 1 ++== （1） 当设计程序时，为了减少存储空间，又由于矩阵A 在经过计算后就不需要再使
用了，因此可以采用原位存储的方法。

也就是可以将k d 存于矩阵A 的对角线元素即kk a 中，ik u 存于矩阵A 的上三角部分即ki a 中，而ik l 则可以存于矩阵A 的下三角部分即ik a 中。

当需要求解对称线性方程组b x A = 时，由T LDL A =分解，有()b x LDL T = ，这可以分解成三个简单的方程组b Ly =，y Dz =和z x L T =由此易得求解公式为
⎪⎪⎪⎩
⎪
⎪⎪⎨⎧
-=-=-===-=∑∑+=-=11 11
2 1
11
1
,,n n,i x l z x ,,n n,i d y z ,n ,,i y l b y n
i k k ki i i i i i i k k ik i i ,,/, （2） 同样，向量b 在经过计算后也不再需要，因此也可以采用相同的方法处理。

可
以在计算i y 时将i y 存于i b 中，然后计算i z 时也将i z 存于i b 中，最后计算i x 时还
是存于i b 中，由i b 输出计算结果即可。

公式（1）和（2）就是T LDL 分解的求解公式。

用T LDL 分解法求解对称方程组的计算量约为341
n ，这比用Gauss 消元法
求解的计算量要少。

对一般的对称方程组，会出现因某个0 i d 而不能进行T LDL 分解的情况，不过当矩阵A 的各阶顺序主子式都不为零时，则T LDL 分解总能进行到底，从而可用T LDL 分解法来求解。

二、 LDL T 分解法程序主要部分（C 语言）
程序中所有矩阵的下标都是从0～n －1，与公式中的1～n 不同，这只是习惯上的不同，其余都是一致的。

程序主要部分有两个模块LDLTDCMP 和LDLTBKSB 。

第一个是T LDL 分解，第二个是解方程组，A 经LDLTDCMP 计算后输出给LDLTBKSB 使用。

两个模块分开，可用于解不同b ，相同A 的多个对称线性方程组。

void LDLTDCMP (unsigned int n, double * *a) {
int k; int m; int i;
for (k = 0; k < n; k++) {
/* Calculate d[k], and saved in a[k][k] */ for (m = 0; m < k; m++)
a[k][k] = a[k][k] - a[m][k] * a[k][m]; if (a[k][k] == 0) {
printf ("\n\nERROR: LDL\' decompose failed !!\n");
printf (" Main submatrix of every rank of A cannot be zero.\n\n");
return; }
for (i = k + 1; i < n; i++) {
/* temporary varible u[i][k] is saved in a[k][i]*/ for (m = 0; m < k; m++)
a[k][i] = a[k][i] - a[m][i] * a[k][m]; /* Calculate l[i][k], and saved in a[i][k]*/ a[i][k] = a[k][i] / a[k][k];
} } }
void LDLTBKSB (unsigned int n, double * *a, double *b) {
int i; int k;
for (i = 0; i < n; i++) {
/* Calculate y[i], and saved in b[i] */ for (k = 0; k < i; k++)
b[i] = b[i] - a[i][k] * b[k]; }
for (i = n - 1; i >= 0; i--) {
/* Calculate z[i], and saved in b[i] */ b[i] = b[i] / a[i][i];
/* Calculate x[i], and saved in b[i] */ for (k = i + 1; k < n; k++)
b[i] = b[i] - a[k][i] * b[k]; } }
最终的测试程序还加入了参数的输入输出（都是文件读写），和一些信息提示。

三、 LDL T 分解法程序的检验
用⎪⎪⎪⎪⎪⎪⎭
⎫
⎝⎛=121
2
1211
21
212
11A ，()T b 321-=检验程序，得结果为()T
x 551-=，与实际计算结果一致。

四、 LDL T 分解法程序的使用说明
程序目录结构：
|-c_source C 源代码文件，能编译测试相应的算法函数 |-data 输入数据文本文件
|-func C 的函数代码（可被引用于其他程序中）
|-output 输出文件，包括执行文件、目标文件和输出结果的文本文件 |-spec 说明文挡
1．输入数据文件格式
输入数据文件名为LDLTPARA.DAT（在./data/下），为文本文件，文件内容格式如下：
# 以'#'开头的行将被忽略掉, 是注释行. 不要有空行.
# 第一个有用行必须以'n ='开头, 接着是系数矩阵的维数.
n = 3
# 接着输入增广矩阵（也就是系数矩阵A和b. A在左边, b在右边）.
# 每一行应该有n+1 个数, n个是矩阵A的, 一个b的.
# A b
1 0.5 0.5 1
0.5 1 0.5 -2
0.5 0.5 1 3
# 按维数输完A和b后, 接下来的所有行都会被忽略.
2．输出数据文件
输出数据文件名为LDLT_KEY.DAT（在./output/下），为文本文件，文件内容格式如下：
Original parameters:
n = 3
1 0.5 0.5 1
0.5 1 0.5 -2
0.5 0.5 1 3
After calculated:
1 0.5 0.5 1
0.5 0.75 0.25 -5
0.5 0.333 0.667 5
Answers x[]:
1 -5 5
其中，Original parameters 输出的是方程组原来的参数（维数和增广矩阵），After calculated 输出的是经过计算后矩阵A和b的内容，最后输出计算结果即方程组的解。