Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter 6 Linear systems
? Triangular Problems
? Banded Problems
? Full Problems
? Analysis
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
2
1
3
2
1
333231
2221
11
0
00
b
b
b
x
x
x
lll
ll
l
2323213133
2212122
1111
/)(
/)(
/
lxlxlbx
lxlbx
lbx
???
??
?
Triangular Problems
ii
i
j
jijii lxlbx /
1
1
???
?
???
?
?? ?
?
?
Forward substitution
for i=1:n
x(i)=b(i);
for j=1:i-1
x(i)=x(i)-L(i,j)*x(j);
end
x(i)=x(i)/L(i,i);
end
x(1)=b(1)/L(1,1);
for i=2:n
x(i)=(b(i)-L(i,1:i-1)*x(1:i-1))/L(i,i);
end
Row oriented
floops
Column-oriented
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5
2
6
897
051
002
3
2
1
x
x
x
?
?
?
?
?
?
?
?
??
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
16
1
7
1
3
5
2
89
05
3
2
x
x
31 ?x
5/12 ??x
)5/1(9168 3 ????x 40/713 ??x
)1,:2():2(
0
00
1
11
3113
2112
3
2
32
3332
22
nLxnb
lxb
lxb
lxb
x
x
x
lll
ll
l
nnnnnnn
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
????
?
?
1111 / lbx ?
jjjj lbx /?
),:1():1():1( jnjLxnjbnjb j ?????
saxpy
function x = LTriSol(L,b)
n = length(b);
x = zeros(n,1);
for j=1:n-1
x(j) = b(j)/L(j,j);
b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j);
end
x(n) = b(n)/L(n,n);
The version involves n2 flops,just like the row-oriented,dot
product version developed eariler.
Backward Substitution
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
3
2
1
3
2
1
33
2322
131211
00
0
b
b
b
x
x
x
u
uu
uuu
1131321213
2232322
3333
/)(
/)(
/
uxuxubx
uxubx
ubx
???
??
?
x(n)=b(n)/U(n,n);
for i=n-1:-1:1
x(i)=(b(i)-U(i,i+1:n)*x(i+1:n))/U(i,i)
end
function x = UTriSol(U,b)
n = length(b);
x = zeros(n,1);
for j=n:-1:2
x(j) = b(j)/U(j,j);
b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);
Column-oriented
saxpy
Backward Substitution
Multiple Right-Hand Sides
BLX ? )(,,)(,,kBkLX ?
X=zeros(n,r);
for k=1:r
X(:,k)=LTriSol(L,B(:,k));
end
x = zeros(n,1);
For k=1:r
for j=1:n-1
X(j,k) = B(j,k)/L(j,j);
B(j+1:n,k) = B(j+1:n,k) - L(j+1:n,j)*X(j,k);
end
X(n,k) =B(n,k)/L(n,n);
End
x = zeros(n,1);
for j=1:n-1
for k=1:r
X(j,k) = B(j,k)/L(j,j);
B(j+1:n,k) = B(j+1:n,k) - L(j+1:n,j)*X(j,k);
end
End
for k=1:r
X(n,k) =B(n,k)/L(n,n);
end
Exchange
index j,k
function X = LTriSolM(L,B)
[n,r] = size(B);
X = zeros(n,r);
for j=1:n-1
X(j,1:r) = B(j,1:r)/L(j,j);
B(j+1:n,1:r) = B(j+1:n,1:r) - L(j+1:n,j)*X(j,1:r);
end
X(n,1:r) = B(n,1:r)/L(n,n);
Vectorize,on k”
In high-performance computing environments,meaneuvers
like this are often the key to efficient matrix computation.
% Script File,ShowTri
% Inverse of the 5-by-5 Forsythe Matrix.
n = 5;
L = eye(n,n) - tril(ones(n,n),-1)
X = LTriSolM(L,eye(n,n))
L =
1 0 0 0 0
-1 1 0 0 0
-1 -1 1 0 0
-1 -1 -1 1 0
-1 -1 -1 -1 1
X =
1 0 0 0 0
1 1 0 0 0
2 1 1 0 0
4 2 1 1 0
8 4 2 1 1
Banded Problems
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
???
???
???
???
??
?
0000
000
000
000
000
0000
A
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
???
????
?????
??????
??????
?
0000
000
00
0
A
Tridiagonal matrix Upper Hessenberg matrix
LUA ?
bLyUxLxLUAx
yUx
bLy
?????
?
?
?
?
?
)()(
y=Ltrisol(L,b);
x=Utrisol(U,y);
Tridiagonal systems
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5
4
3
2
1
5
4
3
2
1
55
444
333
222
11
b
b
b
b
b
x
x
x
x
x
de
fde
fde
fde
fd
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5
44
33
22
11
5
4
3
2
55
444
333
222
11
1
1
1
1
1
u
fu
fu
fu
fu
l
l
l
l
de
fde
fde
fde
fd
11
11
:),(
/:)1,(
??
??
?????
????
iiiiiiii
iiiiii
flduufldii
ueluleii
n = length(d);
l = zeros(n,1);
u = zeros(n,1);
u(1) = d(1);
for i=2:n
l(i) = e(i)/u(i-1);
u(i) = d(i) - l(i)*f(i-1);
end
function x = LBiDiSol(l,b)
n = length(b);
x = zeros(n,1);
x(1) = b(1);
for i=2:n
x(i) = b(i) - l(i)*x(i-1);
end
Ly=b
function x = UBiDiSol(u,f,b)
n = length(b);
x = zeros(n,1);
x(n) = b(n)/u(n);
for i=n-1:-1:1
x(i) = (b(i) - f(i)*x(i+1))/u(i);
end
Ux=y
[l,u]=TriDiLU(d,e,f);
y=Lbidisol(l,b);
x=Ubidisol(u,f,y);
Operation Procedure Flops
A=LU TriDiLU 3n
Ly=b LBiDisol 2n
Ux=y UBiDiSol 3n
Hessenberg Systems
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
???
????
?????
??????
??????
0000
000
00
0
1121 / aa
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
???
????
?????
?????
??????
0000
000
00
0
0
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
???
????
????
?????
??????
0000
000
00
00
0
2232 / aa
function [v,U] = HessLU(A)
[n,n] = size(A);
v = zeros(n,1);
for k=1:n-1
v(k+1) = A(k+1,k)/A(k,k);
A(k+1,k:n) = A(k+1,k:n) - v(k+1)*A(k,k:n);
end
U = triu(A);
requires n2 flops
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
1
1
2
1
v
M
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
1
1
6
6
v
M
UAMM n ?? 11 ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
1
1
1
1
1
1
1
1
2
1
2
1
1
vv
M
1
1
1
1
?
?
??
?
nMML
LUA
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
1
1
1
1
1
1
6
5
4
3
2
1
5
1
4
1
3
1
2
1
1
v
v
v
v
v
MMMMM
[v,U]=HessLU(A);
y=Lbidisol(v,b);
x=UTriSol(U,y);
Operation Procedure Flops
A=LU HessLU n2
Ly=b LBiDisol 2n
Ux=y UTriDiSol n2
Full Problems
3716136
28564
1332
321
321
321
???
?????
???
xxx
xxx
xxx
2716
24
1332
32
32
321
???
???
???
xx
xx
xxx
63
24
1332
3
32
321
?
???
???
x
xx
xxx
32/)33(
14/)2(
23/6
231
32
3
????
?????
??
xxx
xx
x
Gaussian elimination
LUA ?
LUA ?
?
?
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
3
14
312
143
12
1
16136
564
312
Unit lower
triangular
,
37
28
13
143
12
1
3
2
1
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? y
y
y
y
Matrix factorizations
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??????
??????
??????
??????
??????
??????
111 / aai
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
?????
?????
?????
?????
??????
0
0
0
0
0
222 / aai
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
????
????
????
????
??????
00
00
00
00
00
for k=1:n-1
computer the multipliers required to eliminate x(k)
form equations k+1 through n and store in v(k+1:n).
Update equations k+1 through n.
end
for i=k+1:n
v(i)=A(i,k)/A(k,k);
end
v(k+1:n)=A(k+1:n,k)/A(k,k)
for k=1:n-1
v(k+1:n)=A(k+1:n,k)/A(k,k);
for i=k+1:n
A(i,k+1:n)=A(i,k+1:n)-v(i)*A(k,k+1:n);
end
end
U=triu(A)
n=6,k=3
A(4,4:6) ?A(4,4:6) - v(4)A(3,4:6)
A(5,4:6) ?A(5,4:6) - v(5)A(3,4:6)
A(6,4:6) ?A(6,4:6) - v(6)A(3,4:6)
)6:4,3(
)6:4,5(
)6:4,5(
)6:4,4(
)6:4,5(
)6:4,5(
)6:4,4(
)6:4,6:4(
6
5
4
A
v
v
v
A
A
A
A
A
A
A
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
for k=1:n-1
v(k+1:n)=A(k+1:n,k)/A(k,k)
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-v(k+1:n)*A(k,k+1:n);
end
U=triu(A);
Computer L
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
1
1
6
5
4
3
v
v
v
M
UAMM n ?? 11 ?
UMMMA n )( 1 11211 ????? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
1
1
6
5
4
1
3
v
v
v
M
1
1
1
1
?
?
??
?
nMML
LUA
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?????
1
1
1
1
1
1
)5(
6
)4(
6
)3(
6
)2(
6
)1(
6
)4(
5
)3(
5
)2(
5
)1(
5
)3(
4
)2(
4
)1(
4
)2(
3
)1(
3
)1(
2
1
5
1
4
1
3
1
2
1
1
vvvvv
vvvv
vvv
vv
v
MMMMM
function [L,U] = GE(A)
[n,n] = size(A);
for k=1:n-1
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
[L,U]= GE(A);
y = LTriSol(L,b);
x = UtriSol(U,y);
ShowGe.m
Stability
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
??
?
??
?
22
1211
212
1
01
01
11
10
u
uu
lx
x
??
?
??
??
1
1x
??
?
??
??
?
?
?
?
?
?
??
?
??
?
2
1
11
10
2
1
x
x
LUA ?
?
?
?
?
?
?
?
?
???
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
? 1
10
1
1
1
01
11
1
b
x
x
Ax ??
?
?
??
? ??
?
?
?
?
?
?
??
?
??
??
2
1
11
1
2
1 ??
?x
Delta x(1) x(2)
-----------------------------------------------------
1e-002 1.000000000000001 1.000000000000000
1e-004 0.999999999999890 1.000000000000000
1e-006 1.000000000028756 1.000000000000000
1e-008 0.999999993922529 1.000000000000000
1e-010 1.000000082740371 1.000000000000000
1e-012 0.999866855977416 1.000000000000000
1e-014 0.999200722162641 1.000000000000000
1e-016 2.220446049250313 1.000000000000000
1e-018 0.000000000000000 1.000000000000000
NoPivot.m
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
????
??????
63
53
43
33
00
00
00
00
00
a
a
a
a
Pivoting
Wouldn’t it be nice if |a33| was the largest entry in A(3:6,3)?
That would ensure that all the multipliers are less than or
equal to 1.
function [L,U,piv] = GEpiv(A)
[n,n] = size(A); piv = 1:n;
for k=1:n-1
[maxv,r] = max(abs(A(k:n,k))); q = r+k-1;
piv([k q]) = piv([q k]); A([k q],:) = A([q k],:);
if A(k,k) ~= 0
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end; end
L = eye(n,n) + tril(A,-1);
U = triu(A);
The integer vector piv
is a representation of a
permutation matrix p
piv=[ 3 1 5 4 2]
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
00010
01000
10000
00001
00100
P
)(
2
4
5
1
3
5
4
3
2
1
p i vb
b
b
b
b
b
b
b
b
b
b
P ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
:),( p i vAPA ?
GEpiv
LUPA ?
bAx ? PbPAx ?
LUPA ? PbLy ? yUx ?
[L,U,piv] = GEpiv(A);
y = LTriSol(L,b(piv));
x = UtriSol(U,y);
Operation Procedure Flops
PA=LU GEpiv 2n3/3
Ly=b LTriSol n2
Ux=y UTriSol n2
The pivoting amounts to an O(n2) overhead and so the stabilization
purchased by the row swapping does not seriously affect efficiency.
ShowGEpiv.m
The LU Mentality
dAc T 1???
[L,U,piv] = GEpiv(A);
y = LTrisol(L,d(piv));
x = UTriSol(U,y);
alpha = c’*x;
This is more efficient by a factor of 2 because the explicit
formation of A-1 requires about 4n3/3 flops
)1()( ?? jj vAv
b = v0;
V=zeros(n,k);
for j = 1:k
[L,U,piv] = GEpiv(A);
y = LTrisol(L,b(piv));
V(:,k) = UTriSol(U,y);
b = V(:,k);
end
b = v0;
V=zeros(n,k);
[L,U,piv] = GEpiv(A);
for j = 1:k
y = LTrisol(L,b(piv));
V(:,k) = UTriSol(U,y);
b = V(:,k);
end
k(2n3/3) flops 2n3/3+O(kn2 ) flops
The Matlab Linear System Tools
PA=LU [L,U,P] = LU(A)
piv = P*(1:n)’
[L,U,piv] = GEpiv(A);
y = LTriSol(L,b(piv));
x = UtriSol(U,y);
x = A\b
Analysis Residual Versus Error
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
254.
217.
659.913.
563.780.
2
1
x
x
?
?
?
?
?
?
?
?
0 8 7.
3 4 1.)1(
x ?
?
?
?
?
?
?
?
00.1
9 9 9.)2(
x
?
?
?
?
?
?
??
0
0 0 0 0 0 1.)1(
Axb ?
?
?
?
?
?
??
0 0 0 9 1 3.
0 0 0 7 8 0.)2(
Axb
)( ?e x a c tx
?
?
?
?
?
?
?
1
1)( e x a c t
x
Reasoning the face of such a dichotomy requires an
understanding about the context in which the linear system
arises,We may be in a situation where how well Ax predicts
b is paramount,In this case small residuals are important,In
other settings,accuracy is critical and the focus is on
nearness to the true solution.
It is clear from the discussion that (1) the notion of a,good”
solution can be ambiguous and (2) the intelligent appraisal of
algorithms requires a sharper understanding of the mathematics
behind the Ax=b problem
Problem Sensitivity and Nearness
?
?
?
?
?
?
?
659.913.
.,,5 6 3 0 0 1 0 9 5.780.~
A
?
?
?
?
?
?
??
00
...0 0 0 0 0 1 0 9 5.0~
AA
fbxEA ??? ?)(
11 Ae p sE ? 11 be p sf ?
How can we guarantee that A+E is nonsingular?
How close is to the true solution x?x?
Condition number
1
1
11 )(
?? AAA?
1)(1 ?A?
)()( 11 AA ??? ?
s i n g u l ar t ocl o s e is if l ar g e is )(1 AA?
?
?
?
?
?
?
?
?
???
?
?
?
?
?
? ??
7 8 0.9 1 3.
5 6 3.6 5 9.
10
6 5 9.9 1 3.
5 6 3.7 8 0,61
AA
61 10)( ?A?
?
1)(e p s
?
a n dr n o n s i n g u l a
is ? s y s t e ml i n e a r s t o r e d t h e n t h e
1)( e p s
If, W h e r et h a t
a n dr n o n s i n g u l a is s u p p o s e 6 T h e o r e m
`1
`1
1
`1
`1
1
?
?
?
?
?
?
?
?
??
?
?
?
??
?
?
x
x
A
x
xx
f l ( b )xf l ( A )
A
Rx,bb,Ax
RA
n
nn
?
?
.b y b o u n d e dy e s s e n t ia ll is ?in e r r o r r e la t iv e n o r m-1 hes a y t h a t t to
t h u sis t h e o r e m t h eofon c o n t r ib u t im a in T h e
b o u n d,u p p e r t h e t oli t t le sc o n t r ib u t e ?q u o t ie n t t h a t t h es h o w to
p o s s ib le isit t h e n s in g u la r,y n u m e r ic a lln o t isA If
s in g u la r,y n u m e r ic a ll asA t h in k to
ea p p r o p r ia t isit t h e n 1,t oc lo s e isq u a n t it y t h eIf
p la y, t or o lek e y a h a s e p sf a c t o r T h e
11
1
?
??
x
x/x
( A )??
The 2-norm condition number can be computed,at some
cost,using the function cond,A Cheap estimate of 1-norm
condition number can be obtained with condest
Backward Stability
bxEA ?? ?)(GEpiv AE e p s?
xAbxA ?e p s? ?? )(e p s
?
1 Ax
xx
??
?
The first essentially says that the algorithm produces small
residuals compared to the size of A and, The second heuristic
says that if the unit roundoff and condition satisfy
and,then(roughly) ahs approximately t-p correct
digits.
x?
t?? 10ep s
pA 10)(1 ??
Ge2.m
Pascal matrices
Pascal(n)
N Condition of P
4 6.9e+02
8 2.0e+07
12 8.7e+11
16 4.2e+16det(pascal(n))=1
Det(A)=0 A singular
Det(A) ≈0 A approximately singular
bxEA ?? ?)(