This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Outline
Scientiﬁc Computing: An Introductory Survey
1 Existence, Uniqueness, and Conditioning 2 Solving Linear Systems Department of Computer Science
University of Illinois at UrbanaChampaign 3 Special Types of Linear Systems Copyright c 2002. Reproduction permitted
for noncommercial, educational use only. 4 Software for Linear Systems Chapter 2 – Systems of Linear Equations Prof. Michael T. Heath Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Michael T. Heath 1 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Systems of Linear Equations Scientiﬁc Computing 2 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Singularity and Nonsingularity Given m × n matrix A and mvector b, ﬁnd unknown
nvector x satisfying Ax = b n × n matrix A is nonsingular if it has any of following
equivalent properties System of equations asks “Can b be expressed as linear
combination of columns of A?” 1 2 rank(A) = n 4 Solution may or may not exist, and may or may not be
unique det(A) = 0 3 If so, coefﬁcients of linear combination are given by
components of solution vector x Inverse of A, denoted by A−1 , exists For any vector z = 0, Az = 0 For now, we consider only square case, m = n Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems 3 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Existence and Uniqueness Scientiﬁc Computing 4 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Geometric Interpretation
In two dimensions, each equation determines straight line
in plane Existence and uniqueness of solution to Ax = b depend
on whether A is singular or nonsingular Solution is intersection point of two lines
Can also depend on b, but only in singular case
If two straight lines are not parallel (nonsingular), then
intersection point is unique If b ∈ span(A), system is consistent
A
nonsingular
singular
singular b
arbitrary
b ∈ span(A)
b ∈ span(A)
/ Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems If two straight lines are parallel (singular), then lines either
do not intersect (no solution) or else coincide (any point
along line is solution) # solutions
one (unique)
inﬁnitely many
none In higher dimensions, each equation determines
hyperplane; if matrix is nonsingular, intersection of
hyperplanes is unique solution
5 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example: Nonsingularity Scientiﬁc Computing 6 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Example: Singularity 2 × 2 system
2 × 2 system
2x1 + 3x2 = b1
5x1 + 4x2 = b2 Ax = or in matrixvector notation
Ax = x1
b
= 1 =b
x2
b2 With b = 4 7 T , then x = 1 2 Scientiﬁc Computing T , there is no solution T is nonsingular regardless of value of b Michael T. Heath x1
b
= 1 =b
x2
b2 is singular regardless of value of b 23
54 For example, if b = 8 13
solution 23
46 T T With b = 4 8 , x = γ (4 − 2γ )/3 is solution for any
real number γ , so there are inﬁnitely many solutions is unique 7 / 88 Michael T. Heath Scientiﬁc Computing 8 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Vector Norms Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Example: Vector Norms
Drawing shows unit sphere in two dimensions for each
norm Magnitude, modulus, or absolute value for scalars
generalizes to norm for vectors
We will use only pnorms, deﬁned by
1/p n x xi p = p i=1 for integer p > 0 and nvector x
Important special cases
n
i=1 xi 
n
2 1/2
i=1 xi  1norm: x 1 = 2norm: x 2 = ∞norm: x ∞ Norms have following values for vector shown = maxi xi  Michael T. Heath x 9 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds 2 Michael T. Heath = 2.0 x ∞ = 1.6 Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Equivalence of Norms 10 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Properties of Vector Norms In general, for any vector x in Rn , x
However, we also have
√
x 1 ≤ n x 2,
x x < interactive example > Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems = 2.8 1 2 ≤ √ nx 1 ≥x ∞, x For any vector norm 2 ≥x 1 ≤n x ∞ x > 0 if x = 0
γ x = γ  · x for any scalar γ
x+y ≤ x + y In more general treatment, these properties taken as
deﬁnition of vector norm Thus, for given n, norms differ by at most a constant, and
hence are equivalent: if one is small, they must all be
proportionally small. Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems (triangle inequality) ∞ Scientiﬁc Computing Useful variation on triangle inequality
 x − y ≤ x−y 11 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Matrix Norms 12 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Matrix Norms
Matrix norm corresponding to vector 1norm is maximum
absolute column sum Matrix norm corresponding to given vector norm is deﬁned
by
Ax
A = maxx=0
x n A 1 aij  = max
j i=1 Matrix norm corresponding to vector ∞norm is maximum
absolute row sum Norm of matrix measures maximum stretching matrix does
to any vector in given vector norm n A ∞ aij  = max
i j =1 Handy way to remember these is that matrix norms agree
with corresponding vector norms for n × 1 matrix
Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing 13 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Properties of Matrix Norms Scientiﬁc Computing 14 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Condition Number
Condition number of square nonsingular matrix A is
deﬁned by
cond(A) = A · A−1 Any matrix norm satisﬁes
A > 0 if A = 0 By convention, cond(A) = ∞ if A is singular γ A = γ  · A for any scalar γ Since A+B ≤ A + B Matrix norms we have deﬁned also satisfy A · A−1 = AB ≤ A · B
Ax ≤ A · x for any vector x max
x=0 Ax
x · min
x=0 Ax
x −1 condition number measures ratio of maximum stretching to
maximum shrinking matrix does to any nonzero vectors
Large cond(A) means A is nearly singular Michael T. Heath Scientiﬁc Computing 15 / 88 Michael T. Heath Scientiﬁc Computing 16 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Properties of Condition Number Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Computing Condition Number
Deﬁnition of condition number involves matrix inverse, so it
is nontrivial to compute For any matrix A, cond(A) ≥ 1 Computing condition number from deﬁnition would require
much more work than computing solution whose accuracy
is to be assessed For identity matrix, cond(I ) = 1
For any matrix A and scalar γ , cond(γ A) = cond(A)
For any diagonal matrix D = diag(di ), cond(D ) = max di 
min di  In practice, condition number is estimated inexpensively as
byproduct of solution process
Matrix norm A is easily computed as maximum absolute
column sum (or row sum, depending on norm used) < interactive example > Estimating A−1 at low cost is more challenging
Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing 17 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Computing Condition Number, continued Scientiﬁc Computing 18 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Error Bounds
Condition number yields error bound for computed solution
to linear system From properties of norms, if Az = y , then ˆ
Let x be solution to Ax = b, and let x be solution to
ˆ
Ax = b + ∆b z
≤ A−1
y ˆ
If ∆x = x − x, then and bound is achieved for optimally chosen y ˆ
b + ∆b = A(x) = A(x + ∆x) = Ax + A∆x Efﬁcient condition estimators heuristically pick y with large
ratio z / y , yielding good estimate for A−1 which leads to bound
∆x
∆b
≤ cond(A)
x
b Good software packages for linear systems provide
efﬁcient and reliable condition estimator for possible relative change in solution x due to relative
change in righthand side b
< interactive example >
Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing 19 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Error Bounds, continued Scientiﬁc Computing 20 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Error Bounds – Illustration Similar result holds for relative change in matrix: if
ˆ
(A + E )x = b, then In two dimensions, uncertainty in intersection point of two
lines depends on whether lines are nearly parallel ∆x
E
≤ cond(A)
ˆ
x
A
If input data are accurate to machine precision, then bound
for relative error in solution x becomes
ˆ
x−x
≤ cond(A)
x mach Computed solution loses about log10 (cond(A)) decimal
digits of accuracy relative to accuracy of input
Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing < interactive example >
21 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Error Bounds – Caveats Scientiﬁc Computing 22 / 88 Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Residual
ˆ
Residual vector of approximate solution x to linear system
Ax = b is deﬁned by Normwise analysis bounds relative error in largest
components of solution; relative error in smaller
components can be much larger ˆ
r = b − Ax Componentwise error bounds can be obtained, but
somewhat more complicated ˆ
In theory, if A is nonsingular, then x − x = 0 if, and only
if, r = 0, but they are not necessarily small
simultaneously Conditioning of system is affected by relative scaling of
rows or columns Since Illconditioning can result from poor scaling as well as near
singularity
Rescaling can help the former, but not the latter Michael T. Heath Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing ∆x
r
≤ cond(A)
ˆ
ˆ
x
A·x small relative residual implies small relative error in
approximate solution only if A is wellconditioned
23 / 88 Michael T. Heath Scientiﬁc Computing 24 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Singularity and Nonsingularity
Norms
Condition Number
Error Bounds Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Residual, continued Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Solving Linear Systems ˆ
If computed solution x exactly satisﬁes
To solve linear system, transform it into one whose solution
is same but easier to compute ˆ
(A + E )x = b
then r
ˆ
x A ≤ What type of transformation of linear system leaves
solution unchanged? E
A so large relative residual implies large backward error in
matrix, and algorithm used to compute solution is unstable We can premultiply (from left) both sides of linear system
Ax = b by any nonsingular matrix M without affecting
solution Stable algorithm yields small relative residual regardless of
conditioning of nonsingular system Solution to M Ax = M b is given by
x = (M A)−1 M b = A−1 M −1 M b = A−1 b Small residual is easy to obtain, but does not necessarily
imply computed solution is accurate
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 25 / 88 Example: Permutations Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 26 / 88 Example: Diagonal Scaling Permutation matrix P has one 1 in each row and column
and zeros elsewhere, i.e., identity matrix with rows or
columns permuted Row scaling: premultiplying both sides of system by
nonsingular diagonal matrix D , DAx = Db, multiplies
each row of matrix and righthand side by corresponding
diagonal entry of D , but solution x is unchanged Note that P −1 = P T
Premultiplying both sides of system by permutation matrix,
P Ax = P b, reorders rows, but solution x is unchanged Column scaling: postmultiplying A by D , ADx = b,
multiplies each column of matrix by corresponding
diagonal entry of D , which rescales original solution Postmultiplying A by permutation matrix, AP x = b,
reorders columns, which permutes components of original
solution x = (AD )−1 b = D −1 A−1 b x = (AP )−1 b = P −1 A−1 b = P T (A−1 b)
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 27 / 88 Triangular Linear Systems Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 28 / 88 Triangular Matrices What type of linear system is easy to solve? Two speciﬁc triangular forms are of particular interest If one equation in system involves only one component of
solution (i.e., only one entry in that row of matrix is
nonzero), then that component can be computed by
division lower triangular : all entries above main diagonal are zero,
aij = 0 for i < j
upper triangular : all entries below main diagonal are zero,
aij = 0 for i > j If another equation in system involves only one additional
solution component, then by substituting one known
component into it, we can solve for other component Successive substitution process described earlier is
especially easy to formulate for lower or upper triangular
systems If this pattern continues, with only one new solution
component per equation, then all components of solution
can be computed in succession. Any triangular matrix can be permuted into upper or lower
triangular form by suitable row permutation System with this property is called triangular
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 29 / 88 ForwardSubstitution Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy xi = bi − ij xj / ii , n i = 2, . . . , n xn = bn /unn , xi = bi − j =1 for j = 1 to n
if j j = 0 then stop
xj = bj / jj
for i = j + 1 to n
bi = bi − ij xj
end
end
Michael T. Heath 30 / 88 Backsubstitution for upper triangular system U x = b i− 1 11 , Scientiﬁc Computing BackSubstitution Forwardsubstitution for lower triangular system Lx = b x1 = b1 / Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems uij xj / uii , i = n − 1, . . . , 1 j =i+1 { loop over columns }
{ stop if matrix is singular }
{ compute solution component }
{ update righthand side } Scientiﬁc Computing 31 / 88 for j = n to 1
if ujj = 0 then stop
xj = bj /ujj
for i = 1 to j − 1
bi = bi − uij xj
end
end
Michael T. Heath { loop backwards over columns }
{ stop if matrix is singular }
{ compute solution component }
{ update righthand side } Scientiﬁc Computing 32 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example: Triangular Linear System Elimination
To transform general linear system into triangular form, we
need to replace selected nonzero entries of matrix by
zeros 2 4 −2 x1
2
0 1
1 x2 = 4
00
4
x3
8 This can be accomplished by taking linear combinations of
rows Using backsubstitution for this upper triangular system,
last equation, 4x3 = 8, is solved directly to obtain x3 = 2 Consider 2vector a =
Next, x3 is substituted into second equation to obtain
x2 = 2 Michael T. Heath Scientiﬁc Computing a1
a2 If a1 = 0, then Finally, both x3 and x2 are substituted into ﬁrst equation to
obtain x1 = −1 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 1
0
−a2 /a1 1 Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 33 / 88 a1
a
=1
a2
0 Elementary Elimination Matrices Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 34 / 88 Elementary Elimination Matrices, continued More generally, we can annihilate all entries below k th
position in nvector a by transformation 1 ···
0
0 ··· 0
a1
a1
. .. . . .
. ..
.
.
.
.
.
. . . . .
.
.
. . . 0 · · ·
1
0 · · · 0 ak ak Mk a = 0 · · · −mk+1 1 · · · 0 ak+1 = 0 . ..
.
. ..
. . . .
.
.
.
. . . . .
.
.
.
.
.
0 ···
−mn 0 · · · 1
an
0 Matrix Mk , called elementary elimination matrix, adds
multiple of row k to each subsequent row, with multipliers
mi chosen so that result is zero
Mk is unit lower triangular and nonsingular
Mk = I − mk eT , where mk = [0, . . . , 0, mk+1 , . . . , mn ]T
k
and ek is k th column of identity matrix
−
−
Mk 1 = I + mk eT , which means Mk 1 = Lk is same as
k
Mk except signs of multipliers are reversed where mi = ai /ak , i = k + 1, . . . , n
Divisor ak , called pivot, must be nonzero
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 35 / 88 Elementary Elimination Matrices, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 36 / 88 Example: Elementary Elimination Matrices 2
For a = 4,
−2 If Mj , j > k , is another elementary elimination matrix, with
vector of multipliers mj , then
Mk Mj 100
2
2
M1 a = −2 1 0 4 = 0
1 0 1 −2
0 = I − mk eT − mj eT + mk eT mj eT
j
j
k
k
= I − mk eT − mj eT
j
k which means product is essentially “union,” and similarly
for product of inverses, Lk Lj Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems and Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 37 / 88 L1 = Michael T. Heath Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 38 / 88 To reduce general linear system Ax = b to upper
triangular form, ﬁrst choose M1 , with a11 as pivot, to
annihilate ﬁrst column of A below ﬁrst row L2 = −
M2 1 1
0
0
0
1
0
=
0 −1/2 1 System becomes M1 Ax = M1 b, but solution is unchanged Next choose M2 , using a22 as pivot, to annihilate second
column of M1 A below second row and 100
M1 M2 = −2 1 0 ,
1 1 /2 1 Scientiﬁc Computing Gaussian Elimination Note that 100 2 1 0 ,
=
−1 0 1 Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example, continued −
M1 1 100
2
2
M2 a = 0 1 0 4 = 4
0 1 /2 1
−2
0 System becomes M2 M1 Ax = M2 M1 b, but solution is still
unchanged 1
0
0
1
0
L1 L 2 = 2
−1 −1/2 1 Scientiﬁc Computing Process continues for each successive column until all
subdiagonal entries have been zeroed 39 / 88 Michael T. Heath Scientiﬁc Computing 40 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Gaussian Elimination, continued Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy LU Factorization
Product Lk Lj is unit lower triangular if k < j , so
−
−1
L = M −1 = M1 1 · · · Mn−1 = L1 · · · Ln−1 Resulting upper triangular linear system
Mn−1 · · · M1 Ax = Mn−1 · · · M1 b is unit lower triangular M Ax = M b By design, U = M A is upper triangular can be solved by backsubstitution to obtain solution to
original linear system Ax = b So we have Process just described is called Gaussian elimination with L unit lower triangular and U upper triangular A = LU Thus, Gaussian elimination produces LU factorization of
matrix into triangular factors
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 41 / 88 LU Factorization, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 42 / 88 Example: Gaussian Elimination
Use Gaussian elimination to solve linear system 2
4 −2 x1
2
9 −3 x2 = 8 = b
Ax = 4
−2 −3
7 x3
10 Having obtained LU factorization, Ax = b becomes
LU x = b, and can be solved by forwardsubstitution in
lower triangular system Ly = b, followed by
backsubstitution in upper triangular system U x = y To annihilate subdiagonal entries of ﬁrst column of A, 100
2
4 −2
2 4 −2
9 −3 = 0 1
1 ,
M1 A = −2 1 0 4
101
−2 −3
7
01
5 100
2
2
M1 b = −2 1 0 8 = 4
101
10
12 Note that y = M b is same as transformed righthand side
in Gaussian elimination
Gaussian elimination and LU factorization are two ways of
expressing same solution process Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 43 / 88 Example, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 44 / 88 Example, continued
We have reduced original system to equivalent upper
triangular system 2 4 −2 x1
2
0 1
1 x2 = 4 = M b
Ux =
00
4
x3
8 To annihilate subdiagonal entry of second column of M1 A, 1
00
2 4 −2
2 4 −2
1 0 0 1
1 = 0 1
1 = U ,
M2 M1 A = 0
0 −1 1
01
5
00
4 which can now be solved by backsubstitution to obtain −1
x = 2
2 1
00
2
2
1 0 4 = 4 = M b
M2 M1 b = 0
0 −1 1
12
8 Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 45 / 88 Example, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 46 / 88 Row Interchanges
Gaussian elimination breaks down if leading diagonal entry
of remaining unreduced matrix is zero at any stage To write out LU factorization explicitly, 100 100
100 2 1 0 0 1 0 = 2 1 0 = L
L1 L2 =
−1 0 1 0 1 1
−1 1 1 Easy ﬁx: if diagonal entry in column k is zero, then
interchange row k with some subsequent row having
nonzero entry in column k and then proceed as usual
If there is no nonzero on or below diagonal in column k ,
then there is nothing to do at this stage, so skip to next
column so that 2
4 −2
100
2 4 −2
4 = 2 1 0 0 1
9 −3
1 = LU
A=
−2 −3
7
−1 1 1
00
4 Zero on diagonal causes resulting upper triangular matrix
U to be singular, but LU factorization can still be completed
Subsequent backsubstitution will fail, however, as it should
for singular matrix Michael T. Heath Scientiﬁc Computing 47 / 88 Michael T. Heath Scientiﬁc Computing 48 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Partial Pivoting Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy LU Factorization with Partial Pivoting
With partial pivoting, each Mk is preceded by permutation
Pk to interchange rows to bring entry of largest magnitude
into diagonal pivot position
Still obtain M A = U , with U upper triangular, but now In principle, any nonzero value will do as pivot, but in
practice pivot should be chosen to minimize error
propagation
To avoid amplifying previous rounding errors when
multiplying remaining portion of matrix by elementary
elimination matrix, multipliers should not exceed 1 in
magnitude M = Mn−1 Pn−1 · · · M1 P1
M −1 L=
is still triangular in general sense, but not
necessarily lower triangular
Alternatively, we can write This can be accomplished by choosing entry of largest
magnitude on or below diagonal as pivot at each stage PA = LU Such partial pivoting is essential in practice for numerically
stable implementation of Gaussian elimination for general
linear systems
< interactive example >
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems where P = Pn−1 · · · P1 permutes rows of A into order
determined by partial pivoting, and now L is lower
triangular Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 49 / 88 Complete Pivoting Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 50 / 88 Example: Pivoting Complete pivoting is more exhaustive strategy in which
largest entry in entire remaining unreduced submatrix is
permuted into diagonal pivot position
Requires interchanging columns as well as rows, leading
to factorization
P AQ = L U Need for pivoting has nothing to do with whether matrix is
singular or nearly singular
For example, with L unit lower triangular, U upper triangular, and P and
Q permutations
Numerical stability of complete pivoting is theoretically
superior, but pivot search is more expensive than for partial
pivoting
Numerical stability of partial pivoting is more than
adequate in practice, so it is almost always used in solving
linear systems by Gaussian elimination
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems is nonsingular yet has no LU factorization unless rows are
interchanged, whereas is singular yet has LU factorization
51 / 88 Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example: Small Pivots Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 52 / 88 Example, continued To illustrate effect of small pivots, consider Using small pivot, and correspondingly large multiplier, has
caused loss of information in transformed matrix
If rows interchanged, then pivot is 1 and multiplier is − , so 1
11 A= where is positive number smaller than mach
If rows are not interchanged, then pivot is and multiplier is
−1/ , so
1
0
10
M=
, L=
,
−1/ 1
1/ 1
1
0 1 − 1/ = 0
1 1
0 −1/ U= 1
0 −1/ Michael T. Heath Scientiﬁc Computing 0
,
1 1
1
0 1− in ﬂoatingpoint arithmetic
Thus,
10
LU =
1 1
=
=A
10 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems 1
− M= in ﬂoatingpoint arithmetic, but then
1
LU =
1/ 11
11 A= Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy U= 01
10 A= L=
= 10
,
1
11
01 11
11
=
01
1 which is correct after permutation Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 53 / 88 Pivoting, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 54 / 88 Residual
ˆ
ˆ
Residual r = b − Ax for solution x computed using
Gaussian elimination satisﬁes Although pivoting is generally required for stability of
Gaussian elimination, pivoting is not required for some
important classes of matrices r
A Diagonally dominant ≤ E
≤ ρ n2
A mach where E is backward error in matrix A and growth factor ρ
is ratio of largest entry of U to largest entry of A n aij  < ajj , ˆ
x j = 1, . . . , n i=1, i=j Without pivoting, ρ can be arbitrarily large, so Gaussian
elimination without pivoting is unstable Symmetric positive deﬁnite
A = AT and xT Ax > 0 for all x = 0 Michael T. Heath Scientiﬁc Computing With partial pivoting, ρ can still be as large as 2n−1 , but
such behavior is extremely rare
55 / 88 Michael T. Heath Scientiﬁc Computing 56 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Residual, continued Example: Small Residual
Use 3digit decimal arithmetic to solve There is little or no growth in practice, so
r
A Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy ˆ
x ≤ E
≈n
A 0.641 0.242
0.321 0.121 mach Gaussian elimination with partial pivoting yields triangular
system which means Gaussian elimination with partial pivoting
yields small relative residual regardless of conditioning of
system 0.641
0.242
0
0.000242 Thus, small relative residual does not necessarily imply
computed solution is close to “true” solution unless system
is wellconditioned Michael T. Heath Scientiﬁc Computing x1
0.883
=
x2
−0.000383 Backsubstitution then gives solution
ˆ
x = 0.782 1.58 T Exact residual for this solution is
−0.000622
ˆ
r = b − Ax =
−0.000202 Complete pivoting yields even smaller growth factor, but
additional margin of stability usually is not worth extra cost Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems x1
0.883
=
x2
0.442 Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 57 / 88 Example, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 58 / 88 Implementation of Gaussian Elimination
Gaussian elimination has general form of triplenested loop Residual is as small as we can expect using 3digit
arithmetic, but exact solution is
x = 1.00 1.00 for
for T for
aij = aij − (aik /akk )akj
end
end
end so error is almost as large as solution
Cause of this phenomenon is that matrix is nearly singular
(cond(A) > 4000)
Division that determines x2 is between two quantities that
are both on order of rounding error, and hence result is
essentially arbitrary
When arbitrary value for x2 is substituted into ﬁrst
equation, value for x1 is computed so that ﬁrst equation is
satisﬁed, yielding small residual, but poor solution
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Indices i, j , and k of for loops can be taken in any order,
for total of 3! = 6 different arrangements
These variations have different memory access patterns,
which may cause their performance to vary widely on
different computers Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 59 / 88 Uniqueness of LU Factorization Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 60 / 88 Elementary elimination matrices Mk , their inverses Lk ,
and permutation matrices Pk used in formal description of
LU factorization process are not formed explicitly in actual
implementation Provided row pivot sequence is same, if we have two LU
ˆˆ
ˆ
ˆ
factorizations P A = LU = LU , then L−1 L = U U −1 = D
is both lower and upper triangular, hence diagonal U overwrites upper triangle of A, multipliers in L overwrite
strict lower triangle of A, and unit diagonal of L need not
be stored ˆ
If both L and L are unit lower triangular, then D must be
ˆ
ˆ
identity matrix, so L = L and U = U
Uniqueness is made explicit in LDU factorization
P A = LDU , with L unit lower triangular, U unit upper
triangular, and D diagonal
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Row interchanges usually are not done explicitly; auxiliary
integer vector keeps track of row order in original locations Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 61 / 88 Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 62 / 88 Even with many righthand sides b, inversion never
overcomes higher initial cost, since each matrixvector
multiplication A−1 b requires n2 operations, similar to cost
of forward and backsubstitution Forward and backsubstitution for single righthandside
vector together require about n2 multiplications and similar
number of additions Inversion gives less accurate answer; for example, solving
3x = 18 by division gives x = 18/3 = 6, but inversion gives
x = 3−1 × 18 = 0.333 × 18 = 5.99 using 3digit arithmetic Can also solve linear system by matrix inversion:
x = A−1 b Matrix inverses often occur as convenient notation in
formulas, but explicit inverse is rarely required to
implement such formulas Computing A−1 is tantamount to solving n linear systems,
requiring LU factorization of A followed by n forward and
backsubstitutions, one for each column of identity matrix Scientiﬁc Computing Scientiﬁc Computing Inversion vs. Factorization LU factorization requires about n3 /3 ﬂoatingpoint
multiplications and similar number of additions n3 , Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Complexity of Solving Linear Systems Michael T. Heath Scientiﬁc Computing Storage Management Despite variations in computing it, LU factorization is
unique up to diagonal scaling of factors Operation count for inversion is about
expensive as LU factorization Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems For example, product A−1 B should be computed by LU
factorization of A, followed by forward and
backsubstitutions using each column of B three times as
63 / 88 Michael T. Heath Scientiﬁc Computing 64 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems GaussJordan Elimination GaussJordan Elimination, continued
GaussJordan elimination requires about n3 /2
multiplications and similar number of additions, 50% more
expensive than LU factorization In GaussJordan elimination, matrix is reduced to diagonal
rather than triangular form
Row combinations are used to annihilate entries above as
well as below diagonal
Elimination matrix used for given column vector a is of form 1
.
.
. 0 0 0 .
.
.
0 ···
..
.
···
···
···
..
.
··· 0
.
.
. −m 1
.
.
. 0
.
.
. 1
0
0
.
.
.
0 −mk−1
1
−mk+1
.
.
.
−mn 0
0
1
.
.
.
0 ···
..
.
···
···
···
..
.
··· During elimination phase, same row operations are also
applied to righthandside vector (or vectors) of system of
linear equations 0
a1
0
. . . .
. .
. . . 0 ak−1 0 0 ak = ak ak+1 0 0 . . .
. . . .
.
.
0
an
1 Once matrix is in diagonal form, components of solution
are computed by dividing each entry of transformed
righthand side by corresponding diagonal entry of matrix
Latter requires only n divisions, but this is not enough
cheaper to offset more costly elimination phase
< interactive example > where mi = ai /ak , i = 1, . . . , n
Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 65 / 88 Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Solving Modiﬁed Problems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 66 / 88 ShermanMorrison Formula
Sometimes refactorization can be avoided even when
matrix does change If righthand side of linear system changes but matrix does
not, then LU factorization need not be repeated to solve
new system ShermanMorrison formula gives inverse of matrix
resulting from rankone change to matrix whose inverse is
already known Only forward and backsubstitution need be repeated for
new righthand side (A − uv T )−1 = A−1 + A−1 u(1 − v T A−1 u)−1 v T A−1 This is substantial savings in work, since additional
triangular solutions cost only O(n2 ) work, in contrast to
O(n3 ) cost of factorization Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems where u and v are nvectors
Evaluation of formula requires O(n2 ) work (for
matrixvector multiplications) rather than O(n3 ) work
required for inversion Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 67 / 88 = A−1 b + A−1 u(1 − v T A−1 u)−1 v T A−1 b 68 / 88 (with 3, 2 entry changed) of system whose LU factorization
was computed in earlier example
One way to choose update vectors is 0
0 0 and v = 1
u=
−2
0 which can be implemented by following steps
−1 Solve Az = u for z , so z = A u
Solve Ay = b for y , so y = A−1 b
Compute x = y + ((v T y )/(1 − v T z ))z If A is already factored, procedure requires only triangular
solutions and inner products, so only O(n2 ) work and no
explicit inverses
Scientiﬁc Computing Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Consider rankone modiﬁcation 2
4 −2 x1
2
4
9 −3 x2 = 8
−2 −1
7 x3
10 x = (A − uv T )−1 b Michael T. Heath Scientiﬁc Computing Example: RankOne Updating of Solution To solve linear system (A − uv T )x = b with new matrix,
use ShermanMorrison formula to obtain Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems RankOne Updating of Solution so matrix of modiﬁed system is A − uv T Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 69 / 88 Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example, continued Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 70 / 88 Scaling Linear Systems
In principle, solution to linear system is unaffected by
diagonal scaling of matrix and righthandside vector Using LU factorization of A to solve Az = u and Ay = b, −3/2
−1 1/2 and y = 2
z=
−1/2
2 In practice, scaling affects both conditioning of matrix and
selection of pivots in Gaussian elimination, which in turn
affect numerical accuracy in ﬁniteprecision arithmetic Final step computes updated solution −1
−3/2
−7
2
vT y 2 + 1/2 = 4
x=y+
z=
1 − vT z
1 − 1/2
2
−1/2
0 It is usually best if all entries (or uncertainties in entries) of
matrix have about same size
Sometimes it may be obvious how to accomplish this by
choice of measurement units for variables, but there is no
foolproof method for doing so in general We have thus computed solution to modiﬁed system
without factoring modiﬁed matrix
Michael T. Heath Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Scientiﬁc Computing Scaling can introduce rounding errors if not done carefully
71 / 88 Michael T. Heath Scientiﬁc Computing 72 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Example: Scaling Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy Iterative Reﬁnement
Given approximate solution x0 to linear system Ax = b,
compute residual
r0 = b − Ax0 Linear system
10
0 x1
1
=
x2 Now solve linear system Az0 = r0 and take has condition number 1/ , so is illconditioned if is small x1 = x0 + z0 If second row is multiplied by 1/ , then system becomes
perfectly wellconditioned as new and “better” approximate solution, since
Ax1 = A(x0 + z0 ) = Ax0 + Az0 Apparent illconditioning was due purely to poor scaling = ( b − r0 ) + r0 = b In general, it is usually much less obvious how to correct
poor scaling Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Process can be repeated to reﬁne solution successively
until convergence, potentially producing solution accurate
to full machine precision Triangular Systems
Gaussian Elimination
Updating Solutions
Improving Accuracy 73 / 88 Iterative Reﬁnement, continued Michael T. Heath Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems 74 / 88 Symmetric Systems
Banded Systems
Iterative Methods Special Types of Linear Systems Iterative reﬁnement requires double storage, since both
original matrix and its LU factorization are required Work and storage can often be saved in solving linear
system if matrix has special properties Due to cancellation, residual usually must be computed
with higher precision for iterative reﬁnement to produce
meaningful improvement Examples include
Symmetric : A = AT , aij = aji for all i, j For these reasons, iterative improvement is often
impractical to use routinely, but it can still be useful in some
circumstances Positive deﬁnite : xT Ax > 0 for all x = 0
Band : aij = 0 for all i − j  > β , where β is bandwidth of A
Sparse : most entries of A are zero For example, iterative reﬁnement can sometimes stabilize
otherwise unstable algorithm
Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Scientiﬁc Computing 75 / 88 Symmetric Systems
Banded Systems
Iterative Methods Symmetric Positive Deﬁnite Matrices l11 l21
0 l22 implies
a11 , l21 = a21 /l11 , Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems l22 = 2
a22 − l21 Scientiﬁc Computing 77 / 88 Symmetric Systems
Banded Systems
Iterative Methods Scientiﬁc Computing 78 / 88 Symmetric Systems
Banded Systems
Iterative Methods Symmetric Indeﬁnite Systems Features of Cholesky algorithm for symmetric positive
deﬁnite matrices For symmetric indeﬁnite A, Cholesky factorization is not
applicable, and some form of pivoting is generally required
for numerical stability All n square roots are of positive numbers, so algorithm is
well deﬁned
No pivoting is required to maintain numerical stability
Only lower triangle of A is accessed, and hence upper
triangular portion need not be stored
Only n3 /6 multiplications and similar number of additions
are required Factorization of form
P AP T = LDLT
with L unit lower triangular and D either tridiagonal or
block diagonal with 1 × 1 and 2 × 2 diagonal blocks, can be
computed stably using symmetric pivoting strategy Thus, Cholesky factorization requires only about half work
and half storage compared with LU factorization of general
matrix by Gaussian elimination, and also avoids need for
pivoting In either case, cost is comparable to that of Cholesky
factorization < interactive example >
Scientiﬁc Computing Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Cholesky Factorization, continued Michael T. Heath Symmetric Systems
Banded Systems
Iterative Methods for j = 1 to n
for k = 1 to j − 1
for i = j to n
aij = aij − aik · ajk
end
end
√
ajj = ajj
for k = j + 1 to n
akj = akj /ajj
end
end where L is lower triangular with positive diagonal entries
Algorithm for computing it can be derived by equating
corresponding entries of A and LLT
In 2 × 2 case, for example, √ 76 / 88 One way to write resulting general algorithm, in which
Cholesky factor L overwrites original matrix A, is A = L LT l11 = Scientiﬁc Computing Cholesky Factorization If A is symmetric and positive deﬁnite, then LU
factorization can be arranged so that U = LT , which gives
Cholesky factorization a11 a21
l
0
= 11
a21 a22
l21 l22 Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems 79 / 88 Michael T. Heath Scientiﬁc Computing 80 / 88 Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Symmetric Systems
Banded Systems
Iterative Methods Band Matrices Symmetric Systems
Banded Systems
Iterative Methods Tridiagonal Matrices
Consider tridiagonal matrix Gaussian elimination for band matrices differs little from
general case — only ranges of loops change b1 a2 A=0 .
.
. Typically matrix is stored in array by diagonals to avoid
storing zero entries
If pivoting is required for numerical stability, bandwidth can
grow (but no more than double) 0 For ﬁxed small bandwidth, band solver can be extremely
simple, especially if pivoting is not required for stability
Michael T. Heath Scientiﬁc Computing 0 b2
..
. c2
..
. ···
..
.
..
. ..
.
··· an−1
0 bn − 1
an 0
.
.
. 0 cn−1 bn Gaussian elimination without pivoting reduces to
d1 = b1
for i = 2 to n
mi = ai /di−1
di = bi − mi ci−1
end General purpose solver for arbitrary bandwidth is similar to
code for Gaussian elimination for general matrices Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems c1 81 / 88 Michael T. Heath Symmetric Systems
Banded Systems
Iterative Methods Tridiagonal Matrices, continued Scientiﬁc Computing Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems 82 / 88 Symmetric Systems
Banded Systems
Iterative Methods General Band Matrices LU factorization of A is then given by 1 m2 L= 0 .
.
.
0 0
1
..
.
.. .
··· ···
..
.
..
.
mn−1
0 ···
.. . 1
mn 0
.
.
. . ,
.
. 0
1 Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems d1 0 .
U =.
.
.
.
.
0 c1 0 d2
..
. c2
..
.
.. ··· .
··· ···
..
.
..
.
dn−1
0 In general, band system of bandwidth β requires O(βn)
storage, and its factorization requires O(β 2 n) work 0
.
.
. 0 cn−1 dn Scientiﬁc Computing Compared with full system, savings is substantial if β 83 / 88 Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems Symmetric Systems
Banded Systems
Iterative Methods Iterative Methods for Linear Systems LAPACK is more recent replacement for LINPACK featuring
higher performance on modern computer architectures,
including some parallel computers For some types of problems, iterative methods have
signiﬁcant advantages over direct methods Both LINPACK and LAPACK are available from Netlib We will study speciﬁc iterative methods later when we
consider solution of partial differential equations
85 / 88 Michael T. Heath
Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems LINPACK and LAPACK
BLAS Basic Linear Algebra Subprograms Level
1 Work
O(n) BLAS encapsulate basic operations on vectors and
matrices so they can be optimized for given computer
architecture while highlevel routines that call them remain
portable 2 O(n2 ) 3 O(n3 ) Higherlevel BLAS encapsulate matrixvector and
matrixmatrix operations for better utilization of memory
hierarchies such as cache and virtual memory with paging
Generic Fortran versions of BLAS are available from
Netlib, and many computer vendors provide custom
versions optimized for their particular systems
Scientiﬁc Computing Scientiﬁc Computing 86 / 88 LINPACK and LAPACK
BLAS Examples of BLAS Highlevel routines in LINPACK and LAPACK are based on
lowerlevel Basic Linear Algebra Subprograms (BLAS) Michael T. Heath LINPACK and LAPACK
BLAS Solving linear systems of such fundamental importance in
scientiﬁc computing that LINPACK has become standard
benchmark for comparing performance of computers In theory, it might take inﬁnite number of iterations to
converge to exact solution, but in practice iterations are
terminated when residual is as small as desired Scientiﬁc Computing 84 / 88 LINPACK is software package for solving wide variety of
systems of linear equations, both general dense systems
and special systems, such as symmetric or banded Iterative methods begin with initial guess for solution and
successively improve it until desired accuracy attained Michael T. Heath Scientiﬁc Computing LINPACK and LAPACK Gaussian elimination is direct method for solving linear
system, producing exact solution in ﬁnite number of steps
(in exact arithmetic) Existence, Uniqueness, and Conditioning
Solving Linear Systems
Special Types of Linear Systems
Software for Linear Systems n Examples
saxpy
sdot
snrm2
sgemv
strsv
sger
sgemm
strsm
ssyrk Function
Scalar × vector + vector
Inner product
Euclidean vector norm
Matrixvector product
Triangular solution
Rankone update
Matrixmatrix product
Multiple triang. solutions
Rankk update Level3 BLAS have more opportunity for data reuse, and
hence higher performance, because they perform more
operations per data item than lowerlevel BLAS
87 / 88 Michael T. Heath Scientiﬁc Computing 88 / 88 ...
View
Full
Document
This note was uploaded on 10/16/2011 for the course MECHANICAL 581 taught by Professor Wasfy during the Fall '11 term at IUPUI.
 Fall '11
 Wasfy
 Mechanical Engineering, Sula

Click to edit the document details