Description
I. The solution to Ax = b may be written b = A−1x. This can be a good way to analyze
algorithms involving linear systems. But we try to avoid forming A−1
explicitly in
computations because it is more than twice as expensive as solving the linear equations.
A good way to form B = A−1
is to solve the matrix equation AB = I. Gauss elimination
applied to A gives A = LU, where the entries of L are the pivots used in elimination.
2
(a) Show that about n
3/3 work reduces AB = I to UB = L
−1
, where the entries of
U and L
−1 are known.
(b) Show that computing the entries of B from UB = L
−1
takes about n
3/2 work.
Hint: It takes one floating point operation (flop) per element for each of the n
elements of the bottom row of B, then two flops per element of the n − 1 row of
B, and so on to the top. The total is n(1 + 2 + · · · + n).
(c) Use this to verify the claim that computing A−1
is more than twice as expensive
as solving Ax = b.
II. Write a program for estimating the condition number of a matrix A. You may use
either the 1-norm or the ∞-norm (or try both and compare the results). You will need
to compute kAk, which is easy, and estimate kA−1k, which is more challenging. From
the properties of the norms, we know that if z is the solution to Az = y, then
kzk = kA
−1
yk ≤ kA
−1
kkyk,
so that
kzk
kyk
≤ kA
−1
k,
and this bound is achieved for some optimally chosen y. Thus, if we choose a y such
that the ratio kzk/kyk is as large as possible, then we will have a reasonable estimate
for kA−1k. Try two different approaches to choosing y:
(a) Choose y as the solution to the system AT
y = c where c is a vector each of
whose components is ±1, with the sign of each component chosen by the following
heuristic. Using the factorization A = LU, the system AT
y = c is solved in
two stages, successively solving U
T
v = c and L
T
y = v. At each step of the
first triangular solution, choose the corresponding component of c to be 1 or −1
depending on which will make the resulting component of v larger in magnitude.
Then solve the second triangular system in the usual way for y.
(b) Choose a small number, say five, different vectors y randomly and use the one
producing the largest ratio kzk/kyk.
Use both of the approaches on each of the following matrices:
A1 =
10 −7 0
−3 2 6
5 −1 5
, A2 =
−73 78 24
92 66 25
−80 37 10
.
How do the results using these two methods compare? To check the quality of your
estimates, compute A−1
explicitly to to determine its true norm (this computation can
also make use of the LU factorization already computed). How does your condition
numbers (the that uses estimated kA−1k as well as the one with exact kA−1k) compares
with the condition number given by Matlab.
3
III. Consider the linear system
1 1 +
1 − 1
x
y
=
1 + (1 + )
1
where is a small parameter to be specified. The exact solution is obviously
x
y
=
1
for any value of .
Use Gaussian elimination to solve the linear system. For = 1.0×10−n
, n = 0, 2, 4, 6, 8,
compute an estimate of the condition number of the matrix (using one or both methods
given in Problem II) and the relative error in each component of the solution. How
accurately is each component computed? How does the accuracy attained for each
component compare with expectations based on the condition number of the matrix.
What conclusions can you draw from this experiment?
IV. An n × n Hilbert matrix H has entries hij = 1/(i + j − 1) so it has the form
1 1/2 1/3 · · ·
1/2 1/3 1/4 · · ·
1/3 1/4 1/5 · · ·
.
.
.
.
.
.
.
.
.
.
.
.
.
Fon n = 2, 3, . . ., generate the Hilbert matrix of order n and also generate the n-vector
b = Hx, where x is the n-vector with all its components equal to 1. Use Gaussian
elimination to solve the resulting linear system Hx = b, obtaining the approximate
solution ˆx. Compute the ∞-norm of the residual r = b−Hxˆ and of the error ∆x = ˆx−x,
where x is the vector of all ones.
(a) How large can you take n before the error is 100%.
(b) Use a condition number estimator to obtain cond(H) for each value of n.
(c) As n varies, how does the number of correct digits in the components of the
computed solution relate to the condition number of the matrix?
V. The determinant of a triangular matrix is equal to the product of its diagonal entries.
Use this fact to develop a program for computing the determinant of an arbitrary
n × n matrix by using its LU factorization. How can you determine the proper sign
for the determinant? To avoid the risk of underflow and overflow, you may wish to
consider computing the logarithm of the determinant instead of the actual value of the
determinant.
