Mathematical Sciences Institute
MATH3512/MATH6112— Matrix Computation
Assignment 2
Thursday 21 September (Week 7), 11.55pm )
Submit online to wattle site. Scanned copies of hand written maths is fine. Computational questions can be uploaded seperately as a PDF file. This can be prepared with software such as Word or LaTeX. (You might consider online resources such as ShareLaTeX). The total point score for this assignment is 45.
Question 1 (Stability) 10pts

{Similar to Question 15.1, Page 112, Trefethen and Bau.}
Higham gives the following model for rounding, and for the floating point operations , corresponding to the real arithmetic operations +, -, ×, ÷ respectively.
(Higham Theorem 2.2) For x ? R , if x lies in the range of the floating point number system F , then
fl(x) = x(1 + d), |d | u.
(Higham (2.4)) For floating point numbers y,z ? F and the floating point operation ~, corresponding to the real arithmetic operation *,
y ~ z = (y * z)(1 + d), |d | 6 u.
For each question below, state whether the algorithm is backward stable, forward stable but not backward stable, or unstable, according to the definition given in the Notes, and prove it or at least give a reasonably convincing argument. Be sure to follow the floating point model given above.
(a) Data: x ? R. 2 pts
Algorithm: Approximate 2x by computing fl(x) ? fl(x).
(b) Data: x ? R. 4 pts
Algorithm: Approximate 2x2 by computing fl(x) ? fl(x) ? fl(x).
(c) Data: x ? R{0}. 4 pts
Algorithm: Approximate 2 by computing fl(x) ? fl(x) fl(x).
Question 2 (Experimental mathematics with random matrices) 15pts
{Based on Question 12.2, page 96, Trefethen and Bau.}
The goal of this problem is to explore some of the properties of random matrices. Your experiments may lead to conjectures and more refined experiments. Do not try to prove anything. Do produce well-designed plots. Define a random real matrix to be an m × m matrix whose entries are independent samples from the normal distribution with mean zero and variance 1v/m. [Using ipython -pylab: A=randn(m,m)/sqrt(m).] The scale factor 1/ m is introduced to make the limiting behaviour clean as m ? 8.
For each of the following attributes of a random matrix, experiment as described below. For the plots, vary m in powers of 2, e.g. m = 2, 4, 8, 16, 32, 64 . . ., and use about 1000 trials per value of m. Your write up should be no more than three pages long.
(a) What do the eigenvalues look like? Try plotting both the complex eigenvalues and histograms for the real and imaginary parts. 2 pts
(b) How does the spectral radius behave as m ? 8? 2 pts
(c) How many real eigenvalues are there? Try plotting histograms of the number of (approximately) real eigenvalues. How does the histogram behave as m increases?
2 pts
(d) What about norms? How does the 2-norm of a random matrix behave as m ? 8. Of course e must have ?(A) = kAk. Does this inequality appear to approach an equality as m ? 8. 2 pts
(e) What is the distribution of condition numbers? Try plotting histograms of the 2-norm condition number. How does the histogram behave as m increases? 2 pts
(f) How do the answers to a) to d) change if we consider random triangular instead of full matrices? 5 pts Question 3 (Determinant computation) 10pts
Consider the calculation of the determinant of a 2 by 2 matrix
-
a b#
x = det= ad - bc.
c d
In floating point arithmetic this calculation can lead to catastrophic cancellation and hence loss of relative accuracy if ad is close to bc .
Let us suppose we have a machine that can perform a fused multiply-add operation with just a single rounding error:
fl(x * y ± z) = (x * y ± z)(1 + d), |d | u.
On this machine we can use this property to create a determinant calculation which has high relative accuracy.
Assuming that the variables a,b,c,d are exactly represented within the floating point system, consider the algorithm
w = b * c, e = w - b * c, x = (a * d - w) + e.
(a) Assuming a machine with a fused multiply-add operation, show that this algorithm satisfies
x = [x(1 + d3) - bcd1(d3 - d2)] (1 + d4)
where x is the computed value of x and the di ’s are the relative floating point errors
obtained when computingD D w , e and x . 3 pts
(b) Show that if |d | = u, then 3 pts
1 u
= 1 +?, where |? | = .
1 + d 1 - u
(c) Using parts a) and b), show that 3 pts
u 2|bc|.
|x - x| = (|x| + |x|) + 2u
1 - u D
(d) Explain why the result in part c) justifies the statement that the algorithm has high relative accuracy. 1 pt Question 4 (Stability of Gaussian Elimination) 10pts
Higham (9.13 and 9.14) gives a historical perspective and notes on LU factorization. Use this as a starting point for a literature search. The reference is available online from the library http://library.anu.edu.au/record=b2340788.
(a) Views on Gaussian Elimination with Partial Pivoting (GEPP) have changed over time. Give a half page summary of the history of the understanding of GEPP, listing the some of the key papers which have changed our understanding. Try to include at least one paper not mentioned in Higham 9.13 or 9.14. 5 pts
(b) Obtain a copy of one of the key papers you listed above, and produce a half page summary of its main arguments, in your own words. Include at least one example illustrating a key concept. 5 pts