Matlab Math
线性方程组
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
Ax = b 有 x = A-1b,但实际上并不显式求 A-1
线性方程组
例子：
7x = 21
x = 21/7=3
如果求逆
x = 7-1× 21 =,142857 × 21 = 2.99997
这就需要一次除和一次乘，且精度更低
Backslash运算符
AX = B X = A\B 左除
XA = B X = B/A 右除
3-by-3的例子
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
6
4
7
515
623
0710
3
2
1
x
x
x
65 5
4623
7710
321
321
21
???
????
??
xxx
xxx
xx
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5.2
1.6
7
55.20
61.00
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1.6
5.2
7
61.00
55.20
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
2.6
5.2
7
2.600
55.20
0710
3
2
1
x
x
x
算法矩阵表示
?
?
?
?
?
?
?
?
?
? ?
?
2.600
55.20
0710
U
?
?
?
?
?
?
?
?
?
?
??
?
104.03.0
015.0
001
L
?
?
?
?
?
?
?
?
?
?
?
010
100
001
P
PALU ?
单位下三角阵 上三角阵 置换阵
置换阵
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0010
0100
0001
1000
P
PA交换 A的行,AP交换 A的列
置换向量 p =[4 1 3 2] PA = A(p,:)
Px = b
x = PTb
置换方程求解
三角阵的求解
Ux = b
x = zeros(n,1);
for k = n:-1:1
x(k) = b(k)/U(k,k);
i = (1:k-1)’;
b(i) = b(i) – x(k)*U(i,k);
end
x = zeros(n,1);
for k = n:-1:1
j = k+1:n
x(k) = (b(k) – U(k,j)*x(j))/U(k,k);
end
求出 x(n),然后把 x(n)消掉 ?
?
?
n
ij
ijij bxu
LU分解
C.F.Gauss,高斯消去法 (GE)
1955~1977,pivot选主元，舍入误差的影响
高斯消去法分为两步：向前消去和向后替换
APMPMPMU nn 112211,,,???
APPPULLL nn 121121,....,?? ?
121
121
.,,
.,,
PPPP
LLLL
n
n
?
?
?
? PALU ?
LU分解例子
?
?
?
?
?
?
?
?
?
?
?
?
?
?
515
623
0710
A
?
?
?
?
?
?
?
?
?
?
?
100
010
001
1P
?
?
?
?
?
?
?
?
?
?
?
104.00
010
001
2M
?
?
?
?
?
?
?
?
?
?
?
010
100
001
2P
?
?
?
?
?
?
?
?
?
?
?
?
105.0
013.0
001
1M
?
?
?
?
?
?
?
?
?
?
?
?
103.0
015.0
001
1L
?
?
?
?
?
?
?
?
?
?
?
?
104.00
010
001
2L
线性方程组求解
Ax = b
LU = PA
Ly = Pb
Ux = y
x = A\b
舍去和舍入
...5 7 1 4 2 8 5 71 4 2 8 5 7 1 4 2 8.3722 ??p
六位有效数字：
舍去浮点表示（ chopped floating-pointing representation）
1103 1 4 2 8 5.0)( ??pfl ch o p
舍入浮点表示（ rounded floating-pointing representation）
1103 1 4 2 8 6.0)( ??pfl ch o p
为何要选主元？
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
6
901.3
7
515
6099.23
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5.2
001.6
7
55.20
6001.00
0710
3
2
1
x
x
x
5位有效数字 在消元时第二个方程乘上 2.5× 103
6.001 × 2.5 × 103 =1.50025 × 104 舍入为 1.5002 × 104,再加上 2.5,
再一次舍入，成为 1.5004 × 104,则第三个方程为
-0.001x2 +(6)(0.99993) = 6.001
x3= 1.5004 × 104 / (1.5005 × 104 ) = 0.99993
1.5005 × 104 x3 = 1.5004 × 104
x2=-1.5 x1= -0.35
真解,(0,-1,1)T
部分选主元
如果每次乘数在量上是小于或等于 1的，计算出来
的解就是精确的。
要使乘数小于或等于 1，可以通过部分选主元来实
现，选要消去列的最大元素到对角元，交换 A和 b
的行（ x不要交换）。
lutx.m, matlab内置 lu分解的一个版本
bslashtx.m,matlab内置 backslash运算
lugui.m,动态显示 LU分解
舍入误差的影响
设计算的解为 x*,真解为 x = A-1b,误差为,e = x – x*,
残量为 r = b – Ax*
???
?
???
??
???
?
???
?
???
?
???
?
254.0
217.0
659.0913.0
563.0780.0
2
1
x
x
???
?
???
?
?
????
???
?
???
? ??
0 0 0 5 4 1.0
0 0 0 4 6 0.0*,
000.1
443.0* Axbrx
,000.1 000.1 ??
?
?
???
?
??x
六位精度算法 (选主元）得到精确解
???
?
???
?
?????
?
???
?
???
?
???
?
0 0 0 0 0 1.0
2 5 4 0 0 0.0
0 0 0 0 0 1.00
6 5 9 0 0 0.09 1 3 0 0 0.0
2
1
x
x
三位精度算法 (选主元）,
分析
在三位精度的算法中,x2是由两个与舍入误差相同阶的数相
除得到的，因此 x2可以是任意的数。然后完全任意的 x2代入
第一个方程得到 x1,所以第一方程的残量是小的。我们也可
以期望第二个方程的残量也是小的，因为矩阵是接近奇异的。
部分主元的高斯消去法可以确保产生小的残量
“小的”舍入误差的阶与三个数量有关：原来矩阵悉数矩阵
单元的大小，消去过程中系数矩阵单元的大小和计算解的元
素的大小。如果其中一个是“大的”，残量在绝对意义上不
必是小的。
即使残量是小的，但不能保证解的误差是小的。它们之间
的关系可以部分由矩阵的 条件数 来刻画
范数
如果 Ax = b,如何测量 x关于 A和 b的敏感性？
????
?
??
?
?? ?
?
pxx
pn
i
p
ip 1,
/1
1
ii
n
i
i
n
i
i xxxxxx m a x,,
2/12
1
2
1
1 ??
?
?
?
???
?
?? ?
??
??
x =(1:4)/5
n1 = norm(x,1),n2 = norm(x)
ninf = norm(x,inf)
条件数
x
Axm
x
AxM m i n,m a x ??
条件数：
1
m in
m a x
)(
?
?? AA
x
Ax
x
Ax
a?
扰动分析
bAx ? ? ? bbxxA ?? ???
xMb ? xmb ?? ?
b
bA
x
x ??? )(?
条件数是相对误差放大因子，也是刻画矩阵奇异性的测度。
???
?
???
?
???
?
?
???
?
???
?
?
???
?
?
0
1
,
7.9
1.4
,
6.67.9
8.21.4
xbA
???
?
???
?
???
?
?
???
?
?
97.0
34.0~
,
70.9
11.4~
xb
,63.1,1
,01.0,8.13
??
??
xx
bb
?
?
,63.1,0 0 0 7 2 4 6.0 ??
x
x
b
b ??
4.2249
0 0 0 7 2 4 6.0
63.1)( ??A?
例子
norm1
例子
假设矩阵的第一个元素 a11=0.1,所有其他元素和 b是整数，且 A
的条件数是 1012。假设用 IEEE双精度 53位浮点算法，且我们
假设可以精确求解计算机中的系统，那么由唯一的误差是由
0.1的二进制表示引起的，我们可以预期
45312 102210 ?? ????
x
x?
也就是说，简单的矩阵系数的存储就可能引起四位有效
位的丧失。
矩阵范数
x
Ax
A m a x? ?? ?? ?
j
iji
i
ijj aAaA m ax,m ax1
J.H.Wilkinson 证明高斯消去法的计算结果 x*精确满足
(A+E)x* = b
???AE
*** xEExAxb ??? ???
?
*
*
xA
Axb
*)(* 1 AxbAxx ??? ? ** 1 xEAxx ???
??? )(* * 1 AEAx xx ??? ?
矩阵范数
cond(A) 或 cond(A,2):使用 svd(A)
cond(A,1)计算 k1(A),使用 inv(A),工作量比 cond(A,2)小
cond(A,inf),同 cond(A’,1)
condest(A) 估计 k1(A),用 lu(A)用 Higham和 Tisseur的算
法，适合大的稀疏矩阵
rcound(A)估计 1/ k1(A),用 lu(A)和一个在 LINPACK和
LAPACK工程中的一个老算法（历史原因）
线性方程组
Cleve Morler著
陈文斌（ wbchen@fudan.edu.cn)
复旦大学 2002
Ax = b 有 x = A-1b,但实际上并不显式求 A-1
线性方程组
例子：
7x = 21
x = 21/7=3
如果求逆
x = 7-1× 21 =,142857 × 21 = 2.99997
这就需要一次除和一次乘，且精度更低
Backslash运算符
AX = B X = A\B 左除
XA = B X = B/A 右除
3-by-3的例子
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
6
4
7
515
623
0710
3
2
1
x
x
x
65 5
4623
7710
321
321
21
???
????
??
xxx
xxx
xx
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5.2
1.6
7
55.20
61.00
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1.6
5.2
7
61.00
55.20
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
2.6
5.2
7
2.600
55.20
0710
3
2
1
x
x
x
算法矩阵表示
?
?
?
?
?
?
?
?
?
? ?
?
2.600
55.20
0710
U
?
?
?
?
?
?
?
?
?
?
??
?
104.03.0
015.0
001
L
?
?
?
?
?
?
?
?
?
?
?
010
100
001
P
PALU ?
单位下三角阵 上三角阵 置换阵
置换阵
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
0010
0100
0001
1000
P
PA交换 A的行,AP交换 A的列
置换向量 p =[4 1 3 2] PA = A(p,:)
Px = b
x = PTb
置换方程求解
三角阵的求解
Ux = b
x = zeros(n,1);
for k = n:-1:1
x(k) = b(k)/U(k,k);
i = (1:k-1)’;
b(i) = b(i) – x(k)*U(i,k);
end
x = zeros(n,1);
for k = n:-1:1
j = k+1:n
x(k) = (b(k) – U(k,j)*x(j))/U(k,k);
end
求出 x(n),然后把 x(n)消掉 ?
?
?
n
ij
ijij bxu
LU分解
C.F.Gauss,高斯消去法 (GE)
1955~1977,pivot选主元，舍入误差的影响
高斯消去法分为两步：向前消去和向后替换
APMPMPMU nn 112211,,,???
APPPULLL nn 121121,....,?? ?
121
121
.,,
.,,
PPPP
LLLL
n
n
?
?
?
? PALU ?
LU分解例子
?
?
?
?
?
?
?
?
?
?
?
?
?
?
515
623
0710
A
?
?
?
?
?
?
?
?
?
?
?
100
010
001
1P
?
?
?
?
?
?
?
?
?
?
?
104.00
010
001
2M
?
?
?
?
?
?
?
?
?
?
?
010
100
001
2P
?
?
?
?
?
?
?
?
?
?
?
?
105.0
013.0
001
1M
?
?
?
?
?
?
?
?
?
?
?
?
103.0
015.0
001
1L
?
?
?
?
?
?
?
?
?
?
?
?
104.00
010
001
2L
线性方程组求解
Ax = b
LU = PA
Ly = Pb
Ux = y
x = A\b
舍去和舍入
...5 7 1 4 2 8 5 71 4 2 8 5 7 1 4 2 8.3722 ??p
六位有效数字：
舍去浮点表示（ chopped floating-pointing representation）
1103 1 4 2 8 5.0)( ??pfl ch o p
舍入浮点表示（ rounded floating-pointing representation）
1103 1 4 2 8 6.0)( ??pfl ch o p
为何要选主元？
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
6
901.3
7
515
6099.23
0710
3
2
1
x
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5.2
001.6
7
55.20
6001.00
0710
3
2
1
x
x
x
5位有效数字 在消元时第二个方程乘上 2.5× 103
6.001 × 2.5 × 103 =1.50025 × 104 舍入为 1.5002 × 104,再加上 2.5,
再一次舍入，成为 1.5004 × 104,则第三个方程为
-0.001x2 +(6)(0.99993) = 6.001
x3= 1.5004 × 104 / (1.5005 × 104 ) = 0.99993
1.5005 × 104 x3 = 1.5004 × 104
x2=-1.5 x1= -0.35
真解,(0,-1,1)T
部分选主元
如果每次乘数在量上是小于或等于 1的，计算出来
的解就是精确的。
要使乘数小于或等于 1，可以通过部分选主元来实
现，选要消去列的最大元素到对角元，交换 A和 b
的行（ x不要交换）。
lutx.m, matlab内置 lu分解的一个版本
bslashtx.m,matlab内置 backslash运算
lugui.m,动态显示 LU分解
舍入误差的影响
设计算的解为 x*,真解为 x = A-1b,误差为,e = x – x*,
残量为 r = b – Ax*
???
?
???
??
???
?
???
?
???
?
???
?
254.0
217.0
659.0913.0
563.0780.0
2
1
x
x
???
?
???
?
?
????
???
?
???
? ??
0 0 0 5 4 1.0
0 0 0 4 6 0.0*,
000.1
443.0* Axbrx
,000.1 000.1 ??
?
?
???
?
??x
六位精度算法 (选主元）得到精确解
???
?
???
?
?????
?
???
?
???
?
???
?
0 0 0 0 0 1.0
2 5 4 0 0 0.0
0 0 0 0 0 1.00
6 5 9 0 0 0.09 1 3 0 0 0.0
2
1
x
x
三位精度算法 (选主元）,
分析
在三位精度的算法中,x2是由两个与舍入误差相同阶的数相
除得到的，因此 x2可以是任意的数。然后完全任意的 x2代入
第一个方程得到 x1,所以第一方程的残量是小的。我们也可
以期望第二个方程的残量也是小的，因为矩阵是接近奇异的。
部分主元的高斯消去法可以确保产生小的残量
“小的”舍入误差的阶与三个数量有关：原来矩阵悉数矩阵
单元的大小，消去过程中系数矩阵单元的大小和计算解的元
素的大小。如果其中一个是“大的”，残量在绝对意义上不
必是小的。
即使残量是小的，但不能保证解的误差是小的。它们之间
的关系可以部分由矩阵的 条件数 来刻画
范数
如果 Ax = b,如何测量 x关于 A和 b的敏感性？
????
?
??
?
?? ?
?
pxx
pn
i
p
ip 1,
/1
1
ii
n
i
i
n
i
i xxxxxx m a x,,
2/12
1
2
1
1 ??
?
?
?
???
?
?? ?
??
??
x =(1:4)/5
n1 = norm(x,1),n2 = norm(x)
ninf = norm(x,inf)
条件数
x
Axm
x
AxM m i n,m a x ??
条件数：
1
m in
m a x
)(
?
?? AA
x
Ax
x
Ax
a?
扰动分析
bAx ? ? ? bbxxA ?? ???
xMb ? xmb ?? ?
b
bA
x
x ??? )(?
条件数是相对误差放大因子，也是刻画矩阵奇异性的测度。
???
?
???
?
???
?
?
???
?
???
?
?
???
?
?
0
1
,
7.9
1.4
,
6.67.9
8.21.4
xbA
???
?
???
?
???
?
?
???
?
?
97.0
34.0~
,
70.9
11.4~
xb
,63.1,1
,01.0,8.13
??
??
xx
bb
?
?
,63.1,0 0 0 7 2 4 6.0 ??
x
x
b
b ??
4.2249
0 0 0 7 2 4 6.0
63.1)( ??A?
例子
norm1
例子
假设矩阵的第一个元素 a11=0.1,所有其他元素和 b是整数，且 A
的条件数是 1012。假设用 IEEE双精度 53位浮点算法，且我们
假设可以精确求解计算机中的系统，那么由唯一的误差是由
0.1的二进制表示引起的，我们可以预期
45312 102210 ?? ????
x
x?
也就是说，简单的矩阵系数的存储就可能引起四位有效
位的丧失。
矩阵范数
x
Ax
A m a x? ?? ?? ?
j
iji
i
ijj aAaA m ax,m ax1
J.H.Wilkinson 证明高斯消去法的计算结果 x*精确满足
(A+E)x* = b
???AE
*** xEExAxb ??? ???
?
*
*
xA
Axb
*)(* 1 AxbAxx ??? ? ** 1 xEAxx ???
??? )(* * 1 AEAx xx ??? ?
矩阵范数
cond(A) 或 cond(A,2):使用 svd(A)
cond(A,1)计算 k1(A),使用 inv(A),工作量比 cond(A,2)小
cond(A,inf),同 cond(A’,1)
condest(A) 估计 k1(A),用 lu(A)用 Higham和 Tisseur的算
法，适合大的稀疏矩阵
rcound(A)估计 1/ k1(A),用 lu(A)和一个在 LINPACK和
LAPACK工程中的一个老算法（历史原因）
