Xref: utzoo sci.math.num-analysis:2169 comp.arch:23331
Path: utzoo!utgpu!news-server.csri.toronto.edu!bonnie.concordia.ca!thunder.mcrcim.mcgill.edu!snorkelwacker.mit.edu!think.com!samsung!olivea!uunet!validgh!dgh
From: dgh@validgh.com (David G. Hough on validgh)
Newsgroups: sci.math.num-analysis,comp.arch
Subject: Condition number of uniform random matrices; residuals
Message-ID: <397@validgh.com>
Date: 18 Jun 91 15:48:18 GMT
Followup-To: poster
Organization: validgh, PO Box 20370, San Jose, CA 95160
Lines: 69

jbs@ibm.com asserts about solving linear equation with
matrices composed of elements from a uniform random distribution:

> Problems
> of size 1000 appear to typically have condition number between 10**5
> and 10**6.

I haven't observed any evidence to support this.  If uniform random 1000x1000
matrices had such condition numbers then it would be impossible to solve
them accurately with Gaussian elimination.  I did report that 
gaussian elimination without pivoting was somewhat unstable on matrices
of that size, but that only reflects on the stability of that method,
rather than the underlying condition of the problem, independent of the
method.

For instance, the normalized residual computed in the linpack program,

        RESIDn = RESID/( N*NORMA*NORMX*EPS )

grows very slowly; typically 1.5-2 for 100x100 systems, and around 10
for 1000x1000 systems, using IEEE arithmetic and conventional partial
pivoting.  So the growth of RESID itself is only slightly faster than linear 
with N.  Of course the size of the residual b-Ax doesn't tell you too much 
about the size of the error in x, in general.   Luck is with us in this case,
however: for single precision matrices, the measured difference in the 
computed x's vs. the nominally correct x's whose elements are all exactly
one is in the same ratio, e.g. about 1e-4 for 1000x1000 vs. 1e-5 for 100x100.
Of course the actual correct answer for this problem isn't a vector of 1's
since the right hand side b was computed by rounding the result of computing
A*(1 1 1 1 1 ...) and so the computed b differs from the one corresponding
to the vector of 1's, and correspondingly the correct x corresponding to
the right hand side b actually used is somewhat different - but fortunately
that doesn't matter much for this well-conditioned problem.

But this is a good illustration of a point I made previously: in general
you don't know the error, because that would be tantamount to knowing the
answer, and most people don't spend most of their time solving problems
with known answers.  So the normalized residual is a good indication of
whether an algorithm has solved a linear equations problem about as well as
can be expected in a particular floating-point precision.  If an answer
with a small error bound is desired,
regardless of the condition of the problem, something else will have to be
done.

This reminds me of a problem that was posed to me several days in advance
of my qualifying examination:

A very eminent applied mathematician and one of his students
had submitted a paper to a prestigious
journal showing experimentally that all the theory about iterative 
improvement must be wrong.

(Iterative improvement is a method for improving the accuracy of the answer
x for linear equation problems that are moderately ill-conditioned.
When the residual can be computed in higher precision, the accuracy of x
can often be improved to nearly full working precision. )

The authors had computed solutions to problems from
several collections of test matrices of varying conditions
and to a number of random matrices and had seldom seen iterative improvement
yield more than a digit of improvement.

What was really going on here?  I couldn't figure it out before the
qualifying exam, but for better or worse I passed anyway.
-- 

David Hough

dgh@validgh.com		uunet!validgh!dgh	na.hough@na-net.ornl.gov