§ 2 三角分解法 /* Matrix Factorization */
? 高斯消元法的矩阵形式 /* Matrix Form of G.E,*/,
Step 1,)0(/
111111 ?? aaam ii
记 L1 =
1
..
.
1
1
1
21
nm
m
?
?,则
?][ )1()1(1 bAL ?
)1(1)1(1)1(11,.,baa n
)2(A )2(b?
Step n ? 1,
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
)(
)2(
2
)1(
1
)(
)2(
2
)2(
22
)1(
1
)1(
12
)1(
11
121
..
.
..
....
...
...
...
n
n
n
nn
n
n
nn
b
b
b
a
aa
aaa
bALLL
?其中
Lk =
1
..
.
1
1
,
,1
kn
kk
m
m
?
?
?
§ 2 Matrix Factorization – Matrix Form of G.E,
??1kL
1
..
.
1
1
,
,1
kn
kk
m
m
?
?? ??? 1 11211,.,nLLL
1
1
1
jim,
记为 L
单位下三角阵
/* unitary lower-triangular matrix */
记 U =
)(
)2(
2
)2(
22
)1(
1
)1(
12
)1(
11
..
..,,
.,,
.,,
n
nn
n
n
a
aa
aaa
LUA ??
A 的 LU 分解
/* LU factorization */
Hey hasn’t GE given me
enough headache? Why
do I have to know its
matrix form!
When you have to
solve the system for
different with a fixed
A,
b?
Could you be more
specific,please?
Factorize A first,then
for every you nly have to
solve two simple triangular
systems and
,
b?
byL ?? ?
yxU ?? ?
§ 2 Matrix Factorization – Matrix Form of G.E,
定理 若 A的所有顺序主子式 /* determinant of leading
principal submatrices */ 均不为 0，则 A 的 LU 分解唯一（其
中 L 为 单位 下三角阵）。
证明,由 § 1中定理可知,LU 分解存在。下面证明唯一性。
若不唯一，则可设 A = L1U1 = L2U2,推出
??121 UU 211122211 LLUULL ??? ?
Upper-triangular
Lower-triangular
With diagonal
entries 1
I? ?
注, L 为一般下三角阵而 U 为 单位 上三角阵的分解称为
Crout 分解 。
实际上只要考虑 A* 的 LU 分解，即,则
即是 A 的 Crout 分解。
ULA ~~* ?
*~*~ LUA ?
§ 2 Matrix Factorization – Doolittle
? 道立特分解法 /* Doolittle Factorization */,
—— LU 分解的紧凑格式 /* compact form */
反复计算,
很浪费哦 ……
通过比较法直接导出 L 和 U 的计算公式。 思路
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
n
nnnn
n
u
uu
l
l
aa
aa
..
.
..
....
...
1...
......
1
1
...
..
.
..
.
..
.
..
.
..,111
1
21
1
111
?
?
),m i n (
1
ji
k
jkki ul
?? jia
§ 2 Matrix Factorization – Doolittle
??? ),m i n( 1 jik jkkiji ula
固定 i,
对 j = i,i+1,…,n 有
ijkj
i
k
ikij uula ?? ?
?
?
1
1
lii = 1
kj
i
k
ikijij ulau ?
?
?
???
1
1
a
将 i, j 对换，对 j = i,i+1,…,n 有
iijiki
i
k
jkji ulula ?? ?
?
?
1
1
ii
i
k
kijkjiji uulal /)(
1
1
??
?
??? b
Algorithm,Doolittle Factorization
Step 1,u1j = a1j; lj1 = aj1 / u11; ( j = 1,…,n )
Step 2,compute and for i = 2,…,n?1;
Step 3,
a b
????? 11nk knnknnnn ulau
一般采用 列主元
法增强稳定性。但注意
也必须做相应的
行交换。
b?
§ 2 Matrix Factorization – Choleski
? 平方根法 /* Choleski’s Method */,
——对称 /* symmetric */ 正定 /* positive definite */
矩阵的分解法
定义 一个矩阵 A = ( aij )n?n 称为 对称阵,如果 aij = aji 。
定义 一个矩阵 A 称为 正定阵,如果 对任意非
零向量 都成立。
0?xAx T ??
x?
?回顾,对称正定阵的几个重要性质
? A?1 亦对称正定，且 aii > 0 ?
若不然，则 0?? ?xA 存在非零解，即
0?xAx T ?? 存在非零解。
?
1111 )()(,???? ???? AAIAAIAA TTT
?
对任意,存在,使得,
即 。
0???x 0???y xyA ???
xAy ?? 1?? 011 ??? ?? yAyyAAAyxAx TTT ??????
?
其中 0?? xAxa Tii ?? Tx )0...1...0(??
第 i 位
? A 的顺序主子阵 /* leading principal submatrices */ Ak 亦对
称正定
对称性显然。对任意 有
,其中 。
kk Rx ?? 0??
0?? xAxxAx TkkTk ???? nk Rxx ????
?
?
???
?? 0
0
????
? A 的特征值 /* eigen value */ ?i > 0
设对应特征值 ? 的非零特征向量
为,则 。 20 xxxxAx TT ????? ?? ???x?
? A 的全部顺序主子式 det ( Ak ) > 0
因为 ?
?
? n
i
iA
1
)de t ( ?
§ 2 Matrix Factorization – Choleski
将 对称 正定阵 A 做 LU 分解
U = uij =
u11 u
ij / uii
1 1
1
u22
unn
记为
UD~
A 对称 TUL ~? 即 TL D LA ?
记 D1/2 =
11u
22u
nnu
Why is uii > 0? Since det(Ak) >
0
2/1~ LDL ?则 仍是下三角阵
TLLA ~~?
nnRL ??
定理 设矩阵 A对称正定，则存在非奇异下三角阵
使得 。若限定 L 对角元为正，则分解唯一。 TLLA ?
注, 对于对称正定阵 A,从 可知对任意 k ? i
有 。 即 L 的元素不会增大，误差可控，不
需选主元。
? ?? ik ikii la 1 2
iiik al ?||
§ 2 Matrix Factorization – Choleski
Algorithm,Choleski’s Method
To factor the symmetric positive definite n?n matrix A into LLT,where L
is lower triangular,
Input,the dimension n; entries aij for 1? i,j ? n of A,
Output,the entries lij for 1? j ? i and 1? i ? n of L,
Step 1 Set ;
Step 2 For j = 2,…,n,set ;
Step 3 For i = 2,…,n?1,do steps 4 and 5
Step 4 Set ;
Step 5 For j = i+1,…,n,set ;
Step 6 Set ;
Step 7 Output ( lij for j = 1,…,i and i = 1,…,n );
STOP,
1111 al ?
1111 / lal jj ?
? ???? 11 2ik ikiiii lal
? ? iiik ikjkjiji lllal ? ???? 11
? ???? 11 2nk nknnnn lal
因为 A 对称，所以只需存半个 A，即
其中
? ? ? ?nnn aaaaannA,..,,...,,,,2/)1( 1222111??
? ?jiiAa ij ???? 2/)1(
运算量为 O( 3/6),比普通 LU
分解少一半，但有 n 次开方。用 A = LDLT
分解，可省开方时间 (p.50-51)。
HW,p.54
#2,#5,#6
§ 2 Matrix Factorization – Tridiagonal System
? 追赶法解 三对角 方程组
/* Crout Reduction for Tridiagonal Linear System */
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
nnnn
nnn
f
f
f
x
x
x
ba
cba
cba
cb
2
1
2
1
111
222
11
Step 1,对 A 作 Crout 分解
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
2
1
n
nn
A
?
?
??
?
?
直接比较等式两边
的元素，可得到计
算公式 （ p.52） 。
Step 2,追 ——即解, fyL ?? ?,
1
11 ?fy ? ),.,,,2()( 1 niyrfy
i
iiii ??? ??
Step 3,赶 ——即解, yxU ?? ? )1,.,,,1(,1 ????? ? nixyxyx iiiinn ?
与 G.E.类似，一旦
i = 0 则算法中断，故并非任何 三对角阵都可以用此方法分解。
§ 2 Matrix Factorization – Tridiagonal System
定理 若 A 为 对角占优 /* diagonally dominant */ 的三对角
阵，且满足,则追赶
法可解以 A 为系数矩阵的方程组。
0,0,0||||,0|||| 11 ?????? iinn caabcb
Hey,what does diagonally
dominant mean It means that the diagonal entries of the matrix are very LARGE,
Well,how large is
LARGE?
They satisfy the following inequality,
?
?
?
ij ijii
aa ||||
注,
? 如果 A 是 严格 对角占优阵，则不要求三对角线上的
所有元素非零。
? 根据不等式
可知：分解过程中，矩阵元素不会过分增大，算法
保证稳定。
? 运算量为 O(6n)。
||||||||||,1|| 1 iiiiiiii abbab ?????? ????
§ 2 Matrix Factorization – Tridiagonal System
Lab 06,Crout Reduction for Tridiagonal Linear Systems
Apply Crout Reduction to solve a given n× n tridiagonal linear
system
Input
There are several sets of inputs,For each set,
The 1st line contains an integer 100 ? n ? 0 which is the size of a
matrix,n = ?1 signals the end of file,
The 2nd line contains n?1 real numbers,
The 3rd line contains n real numbers,
The 4th line contains n?1 real numbers,
The 5th line contains n real numbers,
The numbers are separated by spaces,
naaa,..32
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
nnnn
nnn
f
f
f
x
x
x
ba
cba
cba
cb
2
1
2
1
111
222
11
nbbb,..21
121,.,?nccc
nfff,..21
§ 2 Matrix Factorization – Tridiagonal System
Output
Each entry of the solution is to be printed as in the C fprintf,
fprintf(outfile,"%16.8e\n",x );
If the method fails to give a solution,print the message
, The?Crout?method?failed.\n”, /* here ? represents a space*/
The outputs of two test cases must be seperated by a blank line,
Sample Input
(? represents a space)
5
–1?–1?–1?–1
4?4?4?4?4
–1?–1?–1?–1
100?200?200?200?100
3
–1?–1
0?3?3
–0.5?–1
8?8?7
–1
Sample Output
(? represents a space)
?4.61538462e+001
?8.46153846e+001
?9.23076923e+001
?8.46153846e+001
?4.61538462e+001
The Crout method failed,
? 高斯消元法的矩阵形式 /* Matrix Form of G.E,*/,
Step 1,)0(/
111111 ?? aaam ii
记 L1 =
1
..
.
1
1
1
21
nm
m
?
?,则
?][ )1()1(1 bAL ?
)1(1)1(1)1(11,.,baa n
)2(A )2(b?
Step n ? 1,
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
)(
)2(
2
)1(
1
)(
)2(
2
)2(
22
)1(
1
)1(
12
)1(
11
121
..
.
..
....
...
...
...
n
n
n
nn
n
n
nn
b
b
b
a
aa
aaa
bALLL
?其中
Lk =
1
..
.
1
1
,
,1
kn
kk
m
m
?
?
?
§ 2 Matrix Factorization – Matrix Form of G.E,
??1kL
1
..
.
1
1
,
,1
kn
kk
m
m
?
?? ??? 1 11211,.,nLLL
1
1
1
jim,
记为 L
单位下三角阵
/* unitary lower-triangular matrix */
记 U =
)(
)2(
2
)2(
22
)1(
1
)1(
12
)1(
11
..
..,,
.,,
.,,
n
nn
n
n
a
aa
aaa
LUA ??
A 的 LU 分解
/* LU factorization */
Hey hasn’t GE given me
enough headache? Why
do I have to know its
matrix form!
When you have to
solve the system for
different with a fixed
A,
b?
Could you be more
specific,please?
Factorize A first,then
for every you nly have to
solve two simple triangular
systems and
,
b?
byL ?? ?
yxU ?? ?
§ 2 Matrix Factorization – Matrix Form of G.E,
定理 若 A的所有顺序主子式 /* determinant of leading
principal submatrices */ 均不为 0，则 A 的 LU 分解唯一（其
中 L 为 单位 下三角阵）。
证明,由 § 1中定理可知,LU 分解存在。下面证明唯一性。
若不唯一，则可设 A = L1U1 = L2U2,推出
??121 UU 211122211 LLUULL ??? ?
Upper-triangular
Lower-triangular
With diagonal
entries 1
I? ?
注, L 为一般下三角阵而 U 为 单位 上三角阵的分解称为
Crout 分解 。
实际上只要考虑 A* 的 LU 分解，即,则
即是 A 的 Crout 分解。
ULA ~~* ?
*~*~ LUA ?
§ 2 Matrix Factorization – Doolittle
? 道立特分解法 /* Doolittle Factorization */,
—— LU 分解的紧凑格式 /* compact form */
反复计算,
很浪费哦 ……
通过比较法直接导出 L 和 U 的计算公式。 思路
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
n
nnnn
n
u
uu
l
l
aa
aa
..
.
..
....
...
1...
......
1
1
...
..
.
..
.
..
.
..
.
..,111
1
21
1
111
?
?
),m i n (
1
ji
k
jkki ul
?? jia
§ 2 Matrix Factorization – Doolittle
??? ),m i n( 1 jik jkkiji ula
固定 i,
对 j = i,i+1,…,n 有
ijkj
i
k
ikij uula ?? ?
?
?
1
1
lii = 1
kj
i
k
ikijij ulau ?
?
?
???
1
1
a
将 i, j 对换，对 j = i,i+1,…,n 有
iijiki
i
k
jkji ulula ?? ?
?
?
1
1
ii
i
k
kijkjiji uulal /)(
1
1
??
?
??? b
Algorithm,Doolittle Factorization
Step 1,u1j = a1j; lj1 = aj1 / u11; ( j = 1,…,n )
Step 2,compute and for i = 2,…,n?1;
Step 3,
a b
????? 11nk knnknnnn ulau
一般采用 列主元
法增强稳定性。但注意
也必须做相应的
行交换。
b?
§ 2 Matrix Factorization – Choleski
? 平方根法 /* Choleski’s Method */,
——对称 /* symmetric */ 正定 /* positive definite */
矩阵的分解法
定义 一个矩阵 A = ( aij )n?n 称为 对称阵,如果 aij = aji 。
定义 一个矩阵 A 称为 正定阵,如果 对任意非
零向量 都成立。
0?xAx T ??
x?
?回顾,对称正定阵的几个重要性质
? A?1 亦对称正定，且 aii > 0 ?
若不然，则 0?? ?xA 存在非零解，即
0?xAx T ?? 存在非零解。
?
1111 )()(,???? ???? AAIAAIAA TTT
?
对任意,存在,使得,
即 。
0???x 0???y xyA ???
xAy ?? 1?? 011 ??? ?? yAyyAAAyxAx TTT ??????
?
其中 0?? xAxa Tii ?? Tx )0...1...0(??
第 i 位
? A 的顺序主子阵 /* leading principal submatrices */ Ak 亦对
称正定
对称性显然。对任意 有
,其中 。
kk Rx ?? 0??
0?? xAxxAx TkkTk ???? nk Rxx ????
?
?
???
?? 0
0
????
? A 的特征值 /* eigen value */ ?i > 0
设对应特征值 ? 的非零特征向量
为,则 。 20 xxxxAx TT ????? ?? ???x?
? A 的全部顺序主子式 det ( Ak ) > 0
因为 ?
?
? n
i
iA
1
)de t ( ?
§ 2 Matrix Factorization – Choleski
将 对称 正定阵 A 做 LU 分解
U = uij =
u11 u
ij / uii
1 1
1
u22
unn
记为
UD~
A 对称 TUL ~? 即 TL D LA ?
记 D1/2 =
11u
22u
nnu
Why is uii > 0? Since det(Ak) >
0
2/1~ LDL ?则 仍是下三角阵
TLLA ~~?
nnRL ??
定理 设矩阵 A对称正定，则存在非奇异下三角阵
使得 。若限定 L 对角元为正，则分解唯一。 TLLA ?
注, 对于对称正定阵 A,从 可知对任意 k ? i
有 。 即 L 的元素不会增大，误差可控，不
需选主元。
? ?? ik ikii la 1 2
iiik al ?||
§ 2 Matrix Factorization – Choleski
Algorithm,Choleski’s Method
To factor the symmetric positive definite n?n matrix A into LLT,where L
is lower triangular,
Input,the dimension n; entries aij for 1? i,j ? n of A,
Output,the entries lij for 1? j ? i and 1? i ? n of L,
Step 1 Set ;
Step 2 For j = 2,…,n,set ;
Step 3 For i = 2,…,n?1,do steps 4 and 5
Step 4 Set ;
Step 5 For j = i+1,…,n,set ;
Step 6 Set ;
Step 7 Output ( lij for j = 1,…,i and i = 1,…,n );
STOP,
1111 al ?
1111 / lal jj ?
? ???? 11 2ik ikiiii lal
? ? iiik ikjkjiji lllal ? ???? 11
? ???? 11 2nk nknnnn lal
因为 A 对称，所以只需存半个 A，即
其中
? ? ? ?nnn aaaaannA,..,,...,,,,2/)1( 1222111??
? ?jiiAa ij ???? 2/)1(
运算量为 O( 3/6),比普通 LU
分解少一半，但有 n 次开方。用 A = LDLT
分解，可省开方时间 (p.50-51)。
HW,p.54
#2,#5,#6
§ 2 Matrix Factorization – Tridiagonal System
? 追赶法解 三对角 方程组
/* Crout Reduction for Tridiagonal Linear System */
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
nnnn
nnn
f
f
f
x
x
x
ba
cba
cba
cb
2
1
2
1
111
222
11
Step 1,对 A 作 Crout 分解
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
2
1
n
nn
A
?
?
??
?
?
直接比较等式两边
的元素，可得到计
算公式 （ p.52） 。
Step 2,追 ——即解, fyL ?? ?,
1
11 ?fy ? ),.,,,2()( 1 niyrfy
i
iiii ??? ??
Step 3,赶 ——即解, yxU ?? ? )1,.,,,1(,1 ????? ? nixyxyx iiiinn ?
与 G.E.类似，一旦
i = 0 则算法中断，故并非任何 三对角阵都可以用此方法分解。
§ 2 Matrix Factorization – Tridiagonal System
定理 若 A 为 对角占优 /* diagonally dominant */ 的三对角
阵，且满足,则追赶
法可解以 A 为系数矩阵的方程组。
0,0,0||||,0|||| 11 ?????? iinn caabcb
Hey,what does diagonally
dominant mean It means that the diagonal entries of the matrix are very LARGE,
Well,how large is
LARGE?
They satisfy the following inequality,
?
?
?
ij ijii
aa ||||
注,
? 如果 A 是 严格 对角占优阵，则不要求三对角线上的
所有元素非零。
? 根据不等式
可知：分解过程中，矩阵元素不会过分增大，算法
保证稳定。
? 运算量为 O(6n)。
||||||||||,1|| 1 iiiiiiii abbab ?????? ????
§ 2 Matrix Factorization – Tridiagonal System
Lab 06,Crout Reduction for Tridiagonal Linear Systems
Apply Crout Reduction to solve a given n× n tridiagonal linear
system
Input
There are several sets of inputs,For each set,
The 1st line contains an integer 100 ? n ? 0 which is the size of a
matrix,n = ?1 signals the end of file,
The 2nd line contains n?1 real numbers,
The 3rd line contains n real numbers,
The 4th line contains n?1 real numbers,
The 5th line contains n real numbers,
The numbers are separated by spaces,
naaa,..32
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
nnnn
nnn
f
f
f
x
x
x
ba
cba
cba
cb
2
1
2
1
111
222
11
nbbb,..21
121,.,?nccc
nfff,..21
§ 2 Matrix Factorization – Tridiagonal System
Output
Each entry of the solution is to be printed as in the C fprintf,
fprintf(outfile,"%16.8e\n",x );
If the method fails to give a solution,print the message
, The?Crout?method?failed.\n”, /* here ? represents a space*/
The outputs of two test cases must be seperated by a blank line,
Sample Input
(? represents a space)
5
–1?–1?–1?–1
4?4?4?4?4
–1?–1?–1?–1
100?200?200?200?100
3
–1?–1
0?3?3
–0.5?–1
8?8?7
–1
Sample Output
(? represents a space)
?4.61538462e+001
?8.46153846e+001
?9.23076923e+001
?8.46153846e+001
?4.61538462e+001
The Crout method failed,
