Computer Problem 3 – Condition numbers and error in linear systems Due in class on June 30, 1997

This Octave assignment relates to the stability of an algorithm for solving linear equations. First you have to know what stability is. When you compute the solution to A · s⃗ = b⃗ numerically using a floating point system you will probably not get the exactly correct answer s⃗, but a slightly different answer which we will call x⃗. We can imagine that x⃗ is the exact solution to (A+E)·x⃗ = b⃗, where

\frac{∥E∥}{∥A∥} ≤ CM_{ϵ}

Here C is some numerical constant and M_{ϵ} is the machine epsilon. We know that the residual r⃗ = b⃗ − A · x⃗ has the property that

\frac{∥E∥}{∥A∥} ≥ \frac{∥r⃗∥}{∥A∥∥x⃗∥} ≡ norm of the relative residual.

Thus if the norm of the relative residual is large then the solution was not computed stably. On the other hand, if the relative residual is small then the result computed stably. To see this we show that there is a matrix E such that (A + E) · x⃗ = b⃗ and that E is in some sense small. Let

E = \frac{r⃗x^{t}}{∥x⃗∥∥x⃗∥}

(here x^{t} is the transpose of x⃗ and r⃗x^{t} is an n × n matrix.)

then

b⃗ − (A + E) · x⃗ = (b⃗ − A · x⃗) − E · x⃗ = r⃗ − \frac{r⃗(x^{t} · x⃗)}{∥x⃗∥∥x⃗∥} = 0

so that (A + E) · x⃗ = b⃗. But

\frac{∥E∥}{∥A∥} = \frac{∥r⃗x^{t}∥}{∥x⃗∥∥x⃗∥∥A∥} = \frac{∥r⃗∥}{∥A∥∥x⃗∥}.

[Since

∥r⃗x^{t}∥^{2} = max \frac{(r⃗x^{t} · z⃗, r⃗x^{t} · z⃗)}{(z⃗, z⃗)} = max \frac{(z⃗, x⃗r^{t} · r⃗x^{t} · z⃗)}{(z⃗, z⃗)} = ∥r⃗∥^{2} max \frac{(z⃗, x⃗x^{t} · z⃗)}{(z⃗, z⃗)} = (∥r⃗∥∥x⃗∥)^{2}.]

The r.h.s. of the above equation, namely:

\frac{∥r⃗∥}{∥A∥∥x⃗∥}

gives a good indication of stability. A stable algorithm will yield a relative residual norm that is the order of M_{ϵ}. An unstable algorithm will yield a much larger value.

We need to be able to generate matrices with a given condition number. We can do this for 2 x 2 matrices by pre and post multiplying by a random orthogonal matrix and its transpose. First we set the random number generator to normal rather than the defualt which is uniform.

In the following example we will use Octave, although Matlab may also be used. Octave is case sensative, while Matlab is not.

rand("normal")
[Q,R] = qr(rand(2))

Q is a random orthogonal matrix.

A = Q*diag([1 1/10**2])*Q'

A is a random matrix and has condition number 10^{2}, which we check:

cond(A)

For larger than a 2 x 2 matrix use for the diagonal:

n = 4, Kappa = 10**2
[Q,R] = qr(rand(n))
A = Q*diag([1 Kappa**(-1/(n-1)) Kappa**(-2/(n-1)) Kappa**(-3/(n-1))])*Q'
cond(A)

To solve a system of equations A · x⃗ = b⃗ by Cramer’s rule you need

A = rand(2), b = rand(2,1)
x1 = (B(1)*A(2,2)-B(2)*A(1,2))/det(A)
x2 = (B(2)*A(1,1)-B(1)*A(2,1))/det(A)

To see differences in the method of solution we can use the command format long

format long
XC = [ x1  x2 ]'

To solve the same system using Gaussian elimination you would use:

XG = A\B

And to solve the system using the inverse you would use:

XI = inv(A)*B

The Octave/Matlab NORM function provides both the 2 (Euclidian) norm of a vector and the spectral norm when applied to a matrix.

Your assignement is for the set of κ = Kap = [1, 10^{4}, 10^{8}, 10^{12}, 10^{16}]

1. For the normal random matrix A with a random x⃗ and n = 4

(a) Generate the random 4 x 4 matrix A with condition number κ_{i} = kap(i)

(b) Generate a normally distributed random 4-vector x⃗:

(c) Calculate b⃗ = A · x⃗.

(d) Calculuate the solution of the system A · x⃗ = b⃗ by Gaussian elimination

(e) Calculuate the solution of the system A·x⃗ = b⃗ by inversion AND then compute and print out for i = 1 : 5 kap(i) and [ i, reg, rrg, rei, rri] where reg = relative error in the solution using Gaussian elimination rrg = relative residual norm for Gaussian elimination rei = relative error in invert and multiply solution rri = relative residual norm for invert and multiply.

2. Now for the case n = 2 using Cramer’s rule and Gaussian elimination print out for i = 1 : 5 kap(i) and [i, reg, rrg, rec, rrc] where reg = relative error in the solution using Gaussian elimination rrg = relative residual norm for Gaussian elimination rec = relative error in the solution by Cramer’s rule rrc = relative residual norm for Cramer’s rule.

3. Put a brief statement in your own words at the top of the file that you hand in containing the results as to what these results mean.

HINTS The relative residual NORM (not error) = \frac{∥r⃗∥}{∥A∥∥x⃗∥}. The relative ERROR is

\frac{∥s⃗−x⃗∥}{∥s⃗∥} where s = exact solution. You generate s⃗, then use A · s⃗ to generate b⃗.

Now solve for x⃗ which satisfies A · x⃗ = b⃗ by the various methods –Cramer’s rule, inverse and Gaussian elimination.

for j=1:4
a(j,:) = [1 1+j 2+j];
endfor
a

A is a 4 x 3 array containing the 4 rows of 3-vectors.

Special remarks for the version of Matlab running in the Lab and on the mainframe: Comments in Matlab are indicated by // while in Octave they are given by % or #. Matlab is not CASE sensative. You must be very careful when you use FOR loops in Matlab. Matlab uses .... to indicate continuation of a loop. Don’t put blank lines inside of FOR loops. Make sure that multiple line FOR loops use .. To include a comment inside a FOR loop, you must start the line with .. // which will continue the loop because of the dots and ignore the rest of the line as a comment because of the slashes. Using ”;” inhibits printout. To see intermediate steps use ”,” instead. In Matlab the following lines of code permit you to create an array of row vectors:

FOR j=1:4; ...
A(j,:) =[ 1 1+j 2+j ]; END;
A

To set rand to the normal distribution use: RAND(’NORMAL’)// To set the format to LONG use: LONG